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

1、数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 211( , ) ()()( )d22x atx atu x tx atx ata 2=0 (11.2.1)( ,0)( ) (11.2.2)( ,0)( ) ttxxtuauu xxu xx (11.2.3) , 0, 0 xta 数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计12( , )()()u x tF xatF xat2()F xat1()F xat数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计11201212120221122 ( )2( )2 ()20 ( ,) xxxxuxxxxxxxxxu

2、xxxxxx x( )0,x)(21)(21),(atxatxtxu 1x 2x 0u x ( )x 图 1 数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计-2-1.5-1-0.500.511.52-4-3-2-101234数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计01212, ( ,)( )0, ( ,)xx xxxx x, 0)(x11( , )( )d( )d22 ()()x atx atu x taaxatxat 数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计1101221020()11( )(

3、)d()()221()()2xxxxxxxxxaaxxxxa 1x 2x x ( )x 图 2 数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计-10-8-6-4-20246810-101-10-8-6-4-20246810-101-10-8-6-4-20246810-101数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计00,0 xx luu02xxttuau00( ),( )tttux ux(0, 0)xlt 11( , )( , )cossinsinnnnnnn atn atABn xu x tux tlll,.)3

4、 , 2 , 1(n数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计,cossinnnnnn xux tNtl22, arctan, nnnnnnnBn aNABAl数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计sinnn xNl120,nlllxln nn2135,2222nllllxnnnn数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计00.10.20.30.40.50.60.70.80.91-10100.10.20.30.40.50.60.70.80.91-10100.10.20.30.40.50.60.

5、70.80.91-10100.10.20.30.40.50.60.70.80.91-101数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 16对代数表达式进行可视化对代数表达式进行可视化边长为边长为b和和c的的四周固定的矩形膜振动的本征值问题为的的四周固定的矩形膜振动的本征值问题为xxyyuuu采用分离变量法可以得到本征模和本征值为采用分离变量法可以得到本征模和本征值为 22sinsinnmmnn xm yXx Yycbmnbc2ttxxyyUaUU数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 17 求前求前9个本征值个本征值 绘制前绘制前4个本征函数的图

6、形个本征函数的图形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*X./b); w21=sin(pi*Y./c).*sin(2*pi*X./b);w22=sin(pi*Y./c).*sin

7、(3*pi*X./b);对二维本征值问题进行可视化对二维本征值问题进行可视化数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 18figure(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 19figure(2)subplot(2,2,1); mesh(X,Y,w

8、11); view(0,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 20figure(3)subplot(2,2,1); contour(X,Y,w11); subplot(2,2,2); contour(X,Y,w12); subplot(2,2,3); contour(

9、X,Y,w21); subplot(2,2,4); contour(X,Y,w22); 00.511.5200.20.40.60.8100.511.5200.20.40.60.8100.511.5200.20.40.60.8100.511.5200.20.40.60.81对二维本征值问题进行可视化对二维本征值问题进行可视化数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计比较下面两副图的区别：比较下面两副图的区别：figure(4)subplot(2,2,1); C,h=contour(X,Y,w11,20); clabel(C,h,manual);subplot(2,2,2); C,h

10、=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,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,

11、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; X,Y=meshgrid(x,y); Z=zeros(51,51);p=moviein(2*3*60); %定义动画帧数定义动画

12、帧数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 endendMOVIE2AVI(p,D:A.avi) sinmnmnTtat代数表达式的代数表达式的动态动态可视化可视化二维本

13、征值问题二维本征值问题数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 23对代数表达式进行可视化对代数表达式进行可视化 二维本征值问题 动态可视化：设具有时间因子 sinmnmnTtat数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 24数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计颜色颜色色素颜色色素颜色结构色结构色数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计两种蝴蝶的翅膀及其表面微结构两种蝴蝶的翅膀及其表面微结构数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计数学物理建

14、模与计算机辅助设计数学物理建模与计算机辅助设计 天然光子晶体天然光子晶体数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计三维周期结构三维周期结构 圣甲虫壳圣甲虫壳 金龟子壳金龟子壳 吉丁虫壳吉丁虫壳数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计矿物三维周期结构矿物三维周期结构 蛋白石蛋白石数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 32数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 33数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 34光子晶体的周期性结构光子晶体的周期性结构(a) 一维光子晶体一维光子晶体(b) 二

15、二维光子晶体维光子晶体(c) 三三维光子晶体维光子晶体数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计光子晶体的电磁场本征值问题光子晶体的电磁场本征值问题Page 3521( )( )( )cE rE rr根据根据Maxwell方程方程，光子晶体中电场满足光子晶体中电场满足： 平面波展开法平面波展开法：将：将E(r)展开为一系列平面波，解上述本征方程，本征值展开为一系列平面波，解上述本征方程，本征值为为(/c)2， 为为Hermitian算符，算符，得到频率得到频率与波矢与波矢k之间的色散关系之间的色散关系。1( )rzyxzzryxzyyrxzyxxrEczyEzxEyExEEcyx

16、EzyEzExEEczxEyxEzEyE222222222222222222222222111数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计Page 36yyyrEcEzEx2222221mzamjmrer21zjknzanjEnyzyeeE22222()yyrEEzm nnmnnkac得到频谱空间中的本征方程：得到频谱空间中的本征方程：一维光子晶体一维光子晶体TE模式的本征方程：模式的本征方程：对对 和和 作傅里叶变换：作傅里叶变换：r1zjknzanjEnyzyeeE2取值范围取值范围 ，为，为z方向波矢，绘制出方向波矢，绘制出 与与 之间的色散曲线。之间的色散曲线。zkczk/

17、 ,/ aa2222()ryyzm nnEEnmnkac2222()ryyzm nnEEnmnkaczkHermitian算符算符本征值本征值一维光子晶体一维光子晶体TE模式本证值问题模式本证值问题数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计 1. 求本征方程：求本征方程：令令 ，特征值为，特征值为 ，特征矢量为，特征矢量为 。 Fourier系数：系数：2.令令 、 、 ，代入代入 中，中，利用一个循环求出利用一个循环求出Q= 再使用：再使用：omega_c=eig(Q); omega_c=sort(sqrt(omega_c);3.求出求出 与与 的的，作出程序及图像作出程序及图

18、像。Page 37rnmznmkanQ2,2yyrEmEnnmnzckan222222cyEnandSincadrnrnr1118 . 0ad25.12randSincadrnrnr111rnmznmkanQ2,2rnmznmkanQ2,2yyrEmEnnmnzckan2222czk求解求解思路思路为为 的函数的函数zkm nQ,数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计clear alld=2; %air gapa=10; %total periodicityeps_r=12.25; d_over_a=d/a;Max=50; %自定义自定义Q=zeros(2*Max+1)fre

19、q=for kz=-pi/a:pi/(10*a):pi/afor 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; endendPage 38一维光子晶体一维光子晶体TE模式本证值问题模式本证值问题数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计omega_c=eig(Q);omega_c=sort(sqrt(omega_c);freq=freq;omega_c.;