上机练习专题6：数学物理方程的计算机求解和可视化.ppt

《上机练习专题6：数学物理方程的计算机求解和可视化.ppt》由会员分享，可在线阅读，更多相关《上机练习专题6：数学物理方程的计算机求解和可视化.ppt（41页珍藏版）》请在装配图网上搜索。

1、数学物理建模与计算机辅助设计,上机练习专题6：数学物理方程的计算机求解和可视化,Page 2,一维无界弦自由振动（无强迫力）齐次定解问题为：,直接利用达朗贝尔公式求解：,（ ）,利用行波法求解数学物理定解问题,达朗贝尔解的物理意义,：正向波（或右行波）,：反向波（或左行波）,传播速度为 a,实例1：一维无界弦振动定解问题:,根据达朗贝尔公式，位移为 ：,clear u0=2; u(1:560)=0; x1=-1;x2=1; x=linspace(-2,2,560); u(141:280)=2*u0*(x(141:280)-x1)/(x2-x1); u(281:420)=2*u0*(x2-x(2

2、81:420)/(x2-x1); uu=u; h=plot(x,u,linewidth,3); axis(2*x1,2*x2,-4,4); set(h,EraseMode,xor) for at=2:140 lu(1:560)=0; ru(1:560)=0; lx=141:420-at; rx=141:420+at; lu(lx)=0.5*uu(141:420); ru(rx)=0.5*uu(141:420); u=lu+ru; set(h,Xdata,x,YData,u); drawnow; pause(0.1) end,某时刻无界弦振动图型,u(x,t),x,实例2：一维无界弦振动定解问题

3、:,由达朗贝尔公式得:,其中:,clear x1=0;x2=1;t=0:0.005:8;x=-10:0.1:10;a=1; X,T=meshgrid(x,t); xpat=X+a*T; xpat(find(xpat=x2)=1; xmat=X-a*T; xmat(find(xmat=x2)=1; jf1=1/2/a*(xpat-x1);jf2=1/2/a*(xmat-x1);jf3=1/2/a*(xpat-x1)-(xmat-x1); subplot(3,1,1) h1=plot(x,jf1(1,:),linewidth,3); set(h1,erasemode,xor); axis(-10

4、10 -1 1) hold on subplot(3,1,2) h2=plot(x,jf2(1,:),linewidth,3); set(h2,erasemode,xor); axis(-10 10 -1 1) hold on subplot(3,1,3) h3=plot(x,jf3(1,:),linewidth,3); set(h3,erasemode,xor); axis(-10 10 -1 1) hold on for j=2:length(t) pause(0.01) set(h1,ydata,jf1(j,:); set(h2,ydata,jf2(j,:); set(h3,ydata,

5、jf3(j,:); drawnow; end,左行波:,右行波:,叠加后的波形:,利用分离变量法求解数学物理定解问题,长为 l ，两端固定的均匀弦的自由振动定解问题：,通解：,特解改写为：,驻波叠加,角频率,初位相,振幅:,波节:,波腹:,function zb t=0:0.005:2.0; x=0:0.001:1; ww1=wfun(1,0);ww2=wfun(2,0); ww3=wfun(3,0);ww4=wfun(4,0); ymax1=max(abs(ww1); figure subplot(4,1,1) h1=plot(x,ww1,r,linewidth,5); axis(0,1,-

6、ymax1,ymax1) subplot(4,1,2) ymax2=max(abs(ww2); h2=plot(x,ww2,r,linewidth,5); axis(0,1,-ymax2,ymax2) subplot(4,1,3) ymax3=max(abs(ww3); h3=plot(x,ww3,r,linewidth,5); axis(0,1,-ymax3,ymax3) subplot(4,1,4) ymax4=max(abs(ww4); h4=plot(x,ww4,r,linewidth,5); axis(0,1,-ymax4,ymax4),for n=2:length(t) ww1=w

7、fun(1,t(n); set(h1,ydata,ww1); ww2=wfun(2,t(n); set(h2,ydata,ww2); drawnow; ww3=wfun(3,t(n); set(h3,ydata,ww3); ww4=wfun(4,t(n); set(h4,ydata,ww4) drawnow; end,function wtx=wfun(k,t) x=0:0.001:1; a=1; wtx=cos(k*pi*a*t).*sin(k*pi*x);,n=1,n=2,n=3,n=4,Page 16,对代数表达式进行可视化,二维本征值问题 矩形区域的本征模与本征振动,边长为b和c的的四

8、周固定的矩形膜振动的本征值问题为,采用分离变量法可以得到本征模和本征值为,Page 17,求前9个本征值 绘制前4个本征函数的图形,b=2; c=1; m,n=meshgrid(1:3); L=(m*pi./c).2+(n*pi./b).2) L = 12.3370 41.9458 91.2938 19.7392 49.3480 98.6960 32.0762 61.6850 111.0330 x=0:0.01:b; y=0:0.01:c; X,Y=meshgrid(x,y); w11=sin(pi*Y./c).*sin(pi*X./b); w12=sin(2*pi*Y./c).*sin(pi

9、*X./b); w21=sin(pi*Y./c).*sin(2*pi*X./b);w22=sin(pi*Y./c).*sin(3*pi*X./b);,对二维本征值问题进行可视化,Page 18,figure(1) subplot(2,2,1); mesh(X,Y,w11); subplot(2,2,2); mesh(X,Y,w12); subplot(2,2,3); mesh(X,Y,w21); subplot(2,2,4); mesh(X,Y,w22);,对二维本征值问题进行可视化,Page 19,figure(2) subplot(2,2,1); mesh(X,Y,w11); view(0

10、,90) subplot(2,2,2); mesh(X,Y,w12); view(0,90) subplot(2,2,3); mesh(X,Y,w21); view(0,90) subplot(2,2,4); mesh(X,Y,w22); view(0,90),对二维本征值问题进行可视化,Page 20,figure(3) subplot(2,2,1); contour(X,Y,w11); subplot(2,2,2); contour(X,Y,w12); subplot(2,2,3); contour(X,Y,w21); subplot(2,2,4); contour(X,Y,w22);,对

11、二维本征值问题进行可视化,比较下面两副图的区别： figure(4) subplot(2,2,1); C,h=contour(X,Y,w11,20); clabel(C,h,manual); subplot(2,2,2); C,h=contour(X,Y,w12,20); clabel(C,h,manual); subplot(2,2,3); C,h=contour(X,Y,w21,20); clabel(C,h,manual); subplot(2,2,4); C,h=contour(X,Y,w22,20); clabel(C,h,manual); figure(5) subplot(2,2

12、,1); C,h=contour(X,Y,w11,20); clabel(C,manual); subplot(2,2,2); C,h=contour(X,Y,w12,20); clabel(C,manual); subplot(2,2,3); C,h=contour(X,Y,w21,20); clabel(C,manual); subplot(2,2,4); C,h=contour(X,Y,w22,20); clabel(C,manual);,Page 21,对二维本征值问题进行可视化,Page 22,动态可视化：设具有时间因子,b=2; c=1; x=0:0.02:b;y=0:0.02:c

13、; X,Y=meshgrid(x,y); Z=zeros(51,51); p=moviein(2*3*60); %定义动画帧数 for m=1:2 for n=1:3 for i=1:60 a=sqrt(m*pi/c).2+(n*pi/b).2); Z=sin(a*i*.02*pi)*sin(m*pi*Y./c).*sin(n*pi*X./b); mesh(X,Y,Z); t=本征振动：,m=,int2str(m), n=,int2str(n); title(t); axis(0 b 0 c -1 1); p(:,(m-1)*3+(n-1)*60+i)=getframe; end end en

14、d MOVIE2AVI(p,D:A.avi),代数表达式的动态可视化二维本征值问题,Page 23,对代数表达式进行可视化,二维本征值问题 动态可视化：设具有时间因子,Page 24,光子晶体专题,光子晶体中的本征值问题 平面波展开法求解并进行可视化,颜色,色素颜色,结构色,结构色的起因,结构色（ structural colour），又称物理色(physical colour)，是一种由光的波长引发的光泽。 物质表面周期性的微观结构对不同波长的光产生色散现象，从而从宏观上观察到不同的颜色形态。,两种蝴蝶的翅膀及其表面微结构,一维周期结构 鲍鱼壳,二维周期结构 蝴蝶翅膀,天然光子晶体,二维周期

15、结构孔雀羽毛,三维周期结构 圣甲虫壳 金龟子壳 吉丁虫壳,矿物三维周期结构 蛋白石,Page 32,Page 33,光子晶体光纤,Page 34,光子晶体的周期性结构,(a) 一维光子晶体,(b) 二维光子晶体,(c) 三维光子晶体,光子晶体具有周期性结构的人工合成（或天然）材料。,光子晶体的电磁场本征值问题,Page 35,根据Maxwell方程，光子晶体中电场满足：,平面波展开法：将E(r)展开为一系列平面波，解上述本征方程，本征值 为(/c)2， 为Hermitian算符，得到频率与波矢k之间的色散关系。,Page 36,得到频谱空间中的本征方程：,一维光子晶体TE模式的本征方程：,对

16、和 作傅里叶变换：,取值范围 ，为z方向波矢，绘制出 与 之间的色散曲线。,Hermitian算符,本征值,一维光子晶体TE模式本证值问题,1. 求本征方程： 令 ，特征值为 ，特征矢量为 。 Fourier系数： 2.令 、 、 ，代入 中， 利用一个循环求出Q= 再使用：omega_c=eig(Q); omega_c=sort(sqrt(omega_c); 3.求出 与 的，作出程序及图像。,Page 37,求解思路,为 的函数,clear all d=2; %air gap a=10; %total periodicity eps_r=12.25; d_over_a=d/a; Max=5

17、0; %自定义 Q=zeros(2*Max+1) freq= for kz=-pi/a:pi/(10*a):pi/a for m=1:2*Max+1 for n=1:2*Max+1 M=m-Max; N=n-Max; kn=(1-1/eps_r)*d_over_a.*sinc(pi*(M-N)*d_over_a)+(M-N)=0)*1/eps_r; Q(m,n)=(2*pi*N/a+kz).2.*kn; end end,Page 38,一维光子晶体TE模式本证值问题,omega_c=eig(Q); omega_c=sort(sqrt(omega_c); freq=freq;omega_c.;

18、end kz=-pi/a:pi/(10*a):pi/a; plot(kz,freq(:,1:10),.-) xlabel(Kz) ylabel(omega/c) title(PBG of 1D PC with d=2,a=10,epsilon_r=12.25),Page 39,一维光子晶体TE模式本证值问题,背景为空气，d/a=0.2，周期性介质层的介电常数为12.25。,Page 40,一维光子晶体TE模式能带曲线,考核内容及注意事项,注：要求交纸质版答卷和电子版答卷。纸质版答卷包括每题的解题思路和程序，12月23日（星期五）上午3、4节课到逸夫楼406教研室交纸质版答卷。电子版答卷于12月23日前发到邮箱397379519。,Page 41,考核内容和比重： 完成课后作业 40% 上机情况 10% 期末考试试题 50%,