求解數(shù)學物理方程的計算機方法課件

1、數(shù)學物理建模與計算機輔助設計第4章 求解數(shù)學物理方程的計算機方法 1本章內(nèi)容4.1 特殊函數(shù)的繪制Gamma函數(shù)的繪制連帶勒讓德函數(shù)的繪制球函數(shù)的圖形的繪制幾種貝塞爾函數(shù)圖形的繪制4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程求解利用PDE工具箱求解數(shù)學物理方程利用差分法求解數(shù)學物理方程Matlab的特殊函數(shù)help matlabspecfun airy- Airy functions. besselj- Bessel function of the first kind. bessely- Bessel function of the second kind. besselh- Bes

2、sel functions of the third kind (Hankel function). besseli- Modified Bessel function of the first kind. besselk- Modified Bessel function of the second kind. beta- Beta function. betainc- Incomplete beta function. betaln- Logarithm of beta function. ellipj- Jacobi elliptic functions. ellipke- Comple

3、te elliptic integral.4.1 特殊函數(shù)的繪制Matlab的特殊函數(shù) erf- Error function. erfc- Complementary error function. erfcx- Scaled complementary error function. erfinv- Inverse error function. expint- Exponential integral function. gamma- Gamma function. gammainc- Incomplete gamma function. gammaln- Logarithm of ga

4、mma function. psi- Psi (polygamma) function. legendre- Associated Legendre function. cross- Vector cross product. dot- Vector dot product.4.1 特殊函數(shù)的繪制函數(shù)(Gamma函數(shù))定義基本性質(zhì)拓展后函數(shù)在除z=0，-1，-2，之外各點解析4.1 特殊函數(shù)的繪制勒讓德(Legendre)函數(shù)問題來由：分離變量法求解拉普拉斯方程分離變量可得歐拉型常微分方程和球諧函數(shù)方程4.1 特殊函數(shù)的繪制勒讓德(Legendre)函數(shù)進一步對球諧函數(shù)分離變量可得到關于函數(shù)的

5、常微分方程此方程稱作l階連帶勒讓德方程。特別地：如果球諧函數(shù)具有旋轉(zhuǎn)對稱性，則解與方位角無關，則m=0，此時方程稱作l階勒讓德方程。4.1 特殊函數(shù)的繪制勒讓德(Legendre)函數(shù)l階勒讓德方程的解為l階連帶勒讓德方程的解為4.1 特殊函數(shù)的繪制勒讓德(Legendre)函數(shù)Matlab計算勒讓德函數(shù)指令 legendre(N,x)結(jié)果為所有N階連帶勒讓德函數(shù)的值例：legendre(2,0.0:0.1:0.2) -0.5000 -0.4850 -0.4400 0 -0.2985 -0.5879 3.0000 2.9700 2.88004.1 特殊函數(shù)的繪制勒讓德(Legendre)函數(shù)繪

6、制前6個勒讓德多項式的圖像(p20_1.m)x=0:0.01:1;y1=legendre(1,x);y2=legendre(2,x);y3=legendre(3,x);y4=legendre(4,x);y5=legendre(5,x);y6=legendre(6,x);plot(x,y1(1,:), x,y2(1,:), x,y3(1,:), x,y4(1,:), x,y5(1,:), x,y6(1,:)title(勒讓德多項式)4.1 特殊函數(shù)的繪制球函數(shù)問題來由：用分離變量法求解亥姆霍茲方程其中令同樣可以分離變量得到兩個方程4.1 特殊函數(shù)的繪制球函數(shù)與球諧函數(shù)方程再分離變量可以得到兩組本

7、征值問題4.1 特殊函數(shù)的繪制球函數(shù)前幾個球函數(shù)4.1 特殊函數(shù)的繪制球函數(shù)球函數(shù)的圖形球函數(shù)是在球面上的二元函數(shù)，它是從球函數(shù)方程的本征問題中得到的本征函數(shù)系。球函數(shù)的圖形是空間圖形，為了畫出其圖形，必須指定球的半徑對復數(shù)形式的球函數(shù)，必須對其實部和虛部分別作圖可以用角變量作為平面上的x、y軸的變量畫出球函數(shù)圖4.1 特殊函數(shù)的繪制球函數(shù)繪制球函數(shù)的程序(p81_1.m)l=3; m=2; R=4; A=3; delta=pi/40;theta0=0:delta:pi; phi=0:2*delta:2*pi;phi,theta=meshgrid(phi, theta0);Ymn=legend

8、re(l,cos(theta0); Ymn=Ymn(m+1,:);L=size(theta,1); yy=repmat(Ymn,1,L);Reyy=yy.*cos(m*phi); Imyy=yy.*sin(m*phi);ReM=max(max(abs(Reyy);Rerho=R+A*Reyy/ReM;Rer=Rerho.*sin(theta); Rex=Rer.*cos(phi);Rey=Rer.*sin(phi); Rez=Rerho.*cos(theta);4.1 特殊函數(shù)的繪制球函數(shù)subplot(1,2,1); surf(Rex,Rey,Rez); light; lighting ph

9、ong;axis(square); axis(-5 5 -5 5 -5 5); axis (off); view(40,30)title(實球諧函數(shù));ImM=max(max(abs(Imyy);Imrho=R+A*Imyy/(ImM+eps*(ImM=0);Imr=Imrho.*sin(theta); Imx=Imr.*cos(phi);Imy=Imr.*sin(phi); Imz=Imrho.*cos(theta);subplot(1,2,2); surf(Imx,Imy,Imz); light; lighting phong;axis(square); axis(-5 5 -5 5 -5

10、 5); axis (off); view(40,30)title(虛球諧函數(shù));4.1 特殊函數(shù)的繪制球函數(shù)figuresubplot(3,1,1); surf(theta,phi,Reyy);xlabel(theta); ylabel(phi);subplot(3,1,2); surf(theta,phi,Imyy);xlabel(theta); ylabel(phi);subplot(3,1,3); surf(theta,phi,(Reyy.2+Imyy.2);xlabel(theta); ylabel(phi);4.1 特殊函數(shù)的繪制球函數(shù)4.1 特殊函數(shù)的繪制球函數(shù)4.1 特殊函數(shù)的

11、繪制球函數(shù)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)問題來由：求解固定邊界的圓膜振動作變量分離方程變成4.1 特殊函數(shù)的繪制貝塞爾函數(shù)分解后的這兩個問題中第一個是對時間變化規(guī)律求解第二個是對空間變化規(guī)律求解對空間方程繼續(xù)分離變量令 ，方程化為v階貝塞爾方程4.1 特殊函數(shù)的繪制貝塞爾函數(shù)貝塞爾方程的解v不為整數(shù)，貝塞爾方程通解為v取任意值，貝塞爾方程通解為v取任意值，由兩類貝塞爾函數(shù)可以構成第三類貝塞爾(柱)函數(shù)，又稱漢克爾函數(shù)第一類貝塞爾(柱)函數(shù)（簡稱貝塞爾函數(shù)）第二類貝塞爾(柱)函數(shù)（又稱諾依曼函數(shù)）4.1 特殊函數(shù)的繪制貝塞爾函數(shù)貝塞爾函數(shù)的表達式特殊貝塞爾函數(shù)：虛宗量貝塞爾函數(shù)4.1 特殊函數(shù)

12、的繪制貝塞爾函數(shù)貝塞爾函數(shù)貝塞爾函數(shù)的計算Besselj第一類貝塞爾函數(shù)，簡稱貝塞爾函數(shù)Bessely第二類貝塞爾函數(shù)，又稱諾依曼函數(shù)Besselh第三類貝塞爾函數(shù)，又稱漢克爾函數(shù)Besseli第一類虛宗量貝塞爾函數(shù)虛宗量貝塞爾函數(shù)Besselk第二類虛宗量貝塞爾函數(shù)虛宗量漢克爾函數(shù)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)例1：繪制貝塞爾函數(shù)曲線y=besselj(0:3,(0:0.2:10);figure(1)plot(0:0.2:10),y)legend(J0,J1,J2,J3)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)補充：貝塞爾函數(shù)的零點的尋找例：尋找0,50之間的零點方法I(插值法) x=0:0.05

13、:50; y=besselj(0,x); LD=; for k=1:length(y)-1 if y(k)*y(k+1)LDLD = Columns 1 through 11 2.4049 5.5201 8.6537 11.7915 14.9309 18.0711 21.2116 24.3525 27.4935 30.6346 33.7758 Columns 12 through 16 36.9171 40.0584 43.1998 46.3412 49.48264.1 特殊函數(shù)的繪制貝塞爾函數(shù)方法II(利用Matlab內(nèi)建函數(shù))，尋找20個零點 j=1; a= ; y=inline(bess

14、elj(0,x),x) for k=1:40 while (besselj(0,j)*besselj(0,j+1)0) j=j+1; end q=fzero(y,j); j=j+1; a=a,q; k=k+1; end4.1 特殊函數(shù)的繪制貝塞爾函數(shù)例2：繪制諾依曼函數(shù)的圖形y=bessely(0:1,(0:0.02:10);plot(0:0.02:10),y)legend(N0,N1)axis (0 10 -3.5 1)grid on4.1 特殊函數(shù)的繪制貝塞爾函數(shù)例3：虛宗量貝塞爾函數(shù)I=besseli(0:1,(0.1:0.3:3);plot(0.1:0.3:3),I)legend(I0

15、,I1)axis (0 3 0 5)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)例4：虛宗量漢克爾函數(shù)K=besselk(0:1,(0.1:0.1:3);plot(0.1:0.1:3),K)legend(K0,K1)axis (0 3 0 10)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)例5：球貝塞爾函數(shù)，球諾依曼函數(shù)，球漢克爾函數(shù)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)(1)球貝塞爾函數(shù) x=eps:0.2:15; y1=sqrt(pi/2./x).*besselj(1/2,x); y2=sqrt(pi/2./x).*besselj(3/2,x); y3=sqrt(pi/2./x).*besselj(5/2,x); y4

16、=sqrt(pi/2./x).*besselj(7/2,x); plot(x,y1,x,y2,x,y3,x,y4) legend(j0,j1, j3,j4)4.1 特殊函數(shù)的繪制貝塞爾函數(shù)(2)球諾依曼函數(shù) x=0.8:0.2:15; y1=sqrt(pi/2./x).*bessely(1/2,x); y2=sqrt(pi/2./x).*bessely(3/2,x); y3=sqrt(pi/2./x).*bessely(5/2,x); y4=sqrt(pi/2./x).*bessely(7/2,x); plot(x,y1,x,y2,x,y3,x,y4) axis (1 10 -0.5 0.4)

17、 legend(N0,N1, N3,N4) grid on4.1 特殊函數(shù)的繪制數(shù)學物理方程的求解方法數(shù)學物理定解問題的求解方法1.行波法；2.分離變量法；3.冪級數(shù)解法；4.格林函數(shù)法； 5.積分變換法；6.保角變換法； 7.變分法；8.計算機仿真解法；9.數(shù)值計算法4.2 數(shù)學物理方程的計算機求解方法利用數(shù)學工具軟件(MATLAB、Mathematic、Mathcad等)和常用的計算機語言(Visual C+等)實現(xiàn)數(shù)學物理方程的求解利用有限差分法、蒙特卡洛方法等方法求解數(shù)學物理方法的數(shù)值解對解析解進行編程例1：二維本征值問題矩形區(qū)域的本征模與本征振動邊長為b和c的的四周固定的矩形本征模的

18、本征值問題為采用分離變量法可以得到本征模和本征值為4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程然后對解析解用編程的方法進行仿真下面繪制前4個本征函數(shù)的圖形 (p70_1.m)a=2; b=1;m,n=meshgrid(1:3);L=(m*pi./b).2+(n*pi./b).2);x=0:0.01:a; y=0:0.01:b;X,Y=meshgrid(x,y);w11=sin(pi*Y./b).*sin(pi*X./a); w12=sin(2*pi*Y./b).*sin(pi*X./a); w21=sin(pi*Y./b).*sin(2*pi*X./a);w22=sin(pi*Y./b

19、).*sin(3*pi*X./a);figuresubplot(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);4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程繼續(xù)可以對波動進行動態(tài)（波動方程圖形）仿真波動問題中的時間因子的形式為這里的不妨取其中的正弦項加以展示4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程function p71_1 (p71_1.m)b

20、=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);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)=getf

21、rame; end endendMOVIE2AVI(p,D:A.avi)4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程例2：三維拉普拉斯方程環(huán)形電流磁感應強度圓線圈的半徑為a，通有電流I，求空間任意一點的磁感應強度B4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程程序p104_1.m4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程例3：泊松方程和格林函數(shù)球域的格林函數(shù) 在半徑為a的導體球內(nèi)(或?qū)w球外)，距球心r0處放置電量為40q的點電荷，求它形成的靜電場定解問題是解析解是4

22、.2 數(shù)學物理方程的計算機求解方法對解析解進行編程(1) 球外有一點電荷(p133_1.m)4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程(2) 球內(nèi)有一點電荷(p135_1.m)4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程例4：泊松方程和格林函數(shù)圓域的格林函數(shù) 在半徑為a的圓外，距圓心r0處放置電量為40q的點電荷，求它形成的靜電場定解問題是問題的解析解為4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程(1) 圓外區(qū)域電荷(p136_1.m)4.2 數(shù)學物理方程的計算機求解方法對解析解進行編程(2) 圓內(nèi)區(qū)域電荷(p137_1.m)4.2 數(shù)學物理方程的計算機求解方法利用P

23、DE工具箱求解數(shù)學物理方程偏微分方程工具箱(PDEToolbox)功能(1)設置PDE(偏微分方程)定解問題 設定二維定解區(qū)域、邊界條件以及方程的形式和系數(shù)(2)用有限元法(FEM)求解PDE 網(wǎng)格的生成、方程的離散以及求出數(shù)值解(3)解的可視化4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程PDE工具箱對偏微分方程的描述：(1)方程類型的描述1)橢圓型方程： 2)拋物型方程： 3)雙曲形方程： 4)特征值方程：4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程PDE工具箱對偏微分方程的描述：(2)邊界條件的描述1) Dilichlet條件2) Neuma

24、nn條件4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程例1：用PDE工具箱求解二維本征值問題(P71)4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程

25、4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程例2：二維波動問題圓形膜的振動(p186)定解問題的解析解為4.2 數(shù)學物理方程的計算機求解方法利用PDE工具箱求解數(shù)學物理方程用PDE工具箱進行仿真4.2 數(shù)學物理方程的計算機求解方法利用差分法求解數(shù)學物理方程導數(shù)的差分公式在x附近將函數(shù)f(x)展開成泰勒公式可以得到進而可以得到前差公式也可得到后差公式4.2 數(shù)學物理方程的計算機求解方法利用差分法求解數(shù)學物理方程由前后差公式可以得到二階導數(shù)的差分公式(中心差分公式)4.2 數(shù)學物理方程的計算機求解方法利用差分法求解數(shù)學物理方程用數(shù)值方法求解微分方程的時候步驟：(1)變量空間網(wǎng)格化(2)微分方程變差分方程(3)定解條件設置(4)整理得到遞推公式4.2 數(shù)學物理方程的計算機求解方法利用差分法求解數(shù)學物理方程例1：一維波動問題兩端固定的弦的振動問題(p163)設邊界條件為4.2 數(shù)學物理方程的計算機求解方法利用差分法求解數(shù)學物理方程解析解結(jié)果4.2 數(shù)學物理方程