MATLAB计算四类数学物理方程的举例求解

数学物理建模与计算机辅助设计第5章四类数学物理方程的求解举例数学物理建模与计算机辅助设计本章内容§5.1求解本征值型数学物理方程§5.2求解稳定型数学物理方程§5.3求解热传导型数学物理方程§5.4求解波动型数学物理方程数学物理建模与计算机辅助设计本征值问题简介•用分离变量法解数学物理方程的时候，最后都会归结到求解本征值问题上去。•在利用本征函数系展开法求解数学物理方程的时候，需要对所用的本征函数系有较好的理解•在本部分中讲介绍一维和二维本征值问题和相应的本征函数系§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题•一维波动方程令则分离变量可得20ttxxuauxl,uxtXxTt2''0;0TaTXX§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题•一维热传导方程令分离变量可得20txxuauxl,uxtXxTt2''0;'0TaTXX§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题•二维拉普拉斯方程令分离变量可得00,0xxyyuuxayb,uxyXxYy12''0;''0XXYY§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题•特点：–分离变量后的方程都是一元一阶或二阶常微分方程–一阶常微分方程表明变量是衰减的，具有阻尼，即解中具有e的负指数因子–二阶常微分方程表明变量是振荡的，将具有一次或多次倍频的正弦或余弦的形式，即为本征函数。这儿的频率既可以指空间频率，也可指时间频率。–下面讲解二阶常微分方程获得的本征函数系,expuxtXxat§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题的本征函数系•本征函数系求解本征值问题注意到这里的两个边界上都是第一类齐次边界条件。其本征值和本征函数分别是000XXXXlsinnxl222,sin,1,2,3,nnxnll§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题的本征函数系•本征函数系求解本征值问题注意到这里的两个边界上都是第二类齐次边界条件。其本征值和本征函数分别是0'0'0XXXXlcosnxl222,cos,0,1,2,nnxnll§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题的本征函数系•本征函数系–求解本征值问题–注意到这里的两个边界中左边界是第一类边界条件，右边界上是第二类齐次边界条件。其本征值和本征函数分别是00'0XXXXl12sinnxl2221212,sin,0,1,2,nnxnll§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计一维本征值问题的本征函数系•本征函数系–求解本征值问题–注意到这里的两个边界中左边界是第二类边界条件，右边界上是第一类齐次边界条件。其本征值和本征函数分别是0'00XXXXl12cosnxl2221212,cos,0,1,2,nnxnll§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计本征函数系图像x=0:0.001:1;(p65)A=sin(pi*[1:4]'*x);B=cos(pi*(0:3)'*x);C=sin(pi*(1/2:7/2)'*x);D=cos(pi*(1/2:7/2)'*x);subplot(4,1,1);plot(x,A);subplot(4,1,2);plot(x,B);subplot(4,1,3);plot(x,C);subplot(4,1,4);plot(x,D);采用了矢量化技术§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计本征函数系运动图像(振动问题)functionzb%存储于文件zb.m中(p65)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));figuresubplot(2,1,1);h1=plot(x,ww1,'r','linewidth',5);holdon;h3=plot(x,ww3,'g','linewidth',5);axis([01–ymax1ymax1]);ymax4=max(abs(ww4));subplot(2,1,2);h2=plot(x,ww2,'r','linewidth',5);holdon;h4=plot(x,ww4,'g','linewidth',5);axis([01–ymax4ymax4]);csoinsnnxaltl§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计本征函数系运动图像(振动问题)forn=2:length(t)ww1=wfun(1,t(n));set(h1,'ydata',ww1);ww3=wfun(3,t(n));set(h3,'ydata',ww3);drawno=wfun(2,t(n));set(h2,'ydata',ww2);ww4=wfun(4,t(n));set(h4,'ydata',ww4);drawnowendfunctionwtx=wfun(k,t)%调用函数存储于wfun.m中x=0:0.001:1;a=1;wtx=cos(k*pi*a*t).*sin(k*pi*x);§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计本征函数系运动图像(热传导问题)functionsfb%存放于sfb.m中(p67)x=0:0.01:1;t=1e-5:0.0001:0.005;ww1=sfbf(1,t(1));ww2=sfbf(2,t(1));ww3=sfbf(3,t(1));ww4=sfbf(4,t(1));ymax3=max(abs(ww3));figuresubplot(2,1,1);h1=plot(x,ww1,'r','linewidth',3);holdon;h3=plot(x,ww3,'g','linewidth',3);axis([01–ymax3ymax3]);ymax4=max(abs(ww4));subplot(2,1,2);h2=plot(x,ww2,'r','linewidth',3);holdon;h4=plot(x,ww4,'g','linewidth',3);axis([01–ymax4ymax4]);2exsinpatnxl§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计本征函数系运动图像(热传导问题)forn=2:length(t)ww1=sfbf(1,t(n));set(h1,'ydata',ww1);ww3=sfbf(3,t(n));set(h3,'ydata',ww3);drawno=sfbf(2,t(n));set(h2,'ydata',ww2);ww4=sfbf(4,t(n));set(h4,'ydata',ww4);drawnowendfunctioncht=sfbf(k,t)%调用函数存储于wfun.m中x=0:0.01:1;cht=sin(k*pi*x).*exp(-k^2*pi^2*2^2*t);§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况•矩形区域的本征模与本征振动边长为b和c的的四周固定的矩形本征模的本征值问题为采用分离变量法可以得到本征模和本征值为xxyyuuu22sinsinnmmnnxmyXxYycbmnbc§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况•下面绘制前4个本征函数的图形(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).*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);§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况•下面绘制前4个本征函数的运动图形：波动图形波动问题中的时间因子的形式为这里的不妨取其中的正弦项加以展示cossinmnmnmnmnmnTtCatDatsinmnmnTtatfunctionp71_1%p71_1.mb=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);form=1:2forn=1:3fori=1:60a=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([0b0c-11]);p(:,((m-1)*3+(n-1))*60+i)=getframe;endendendMOVIE2AVI(p,'D:\A.avi')§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：矩形区域情况§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：圆形区域情况半径为周边固定的圆形膜的本征问题是本征模是其中下面讨论第一类其次边界条件下的本征值情形,,uusin,cosmmJmJm0,0,1,2,,0,1,2,mnmnx0§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：圆形区域情况m=0，n=0,1,2,3的贝塞尔函数的前4个本征函数图形r=0:0.01:1;%(p71_1.m)y=inline('besselj(0,x)','x');x(1)=fzero(y,1);plot(r,y(x(1)*r))holdonfork=1:3x(k+1)=fzero(y,x(k)+3);plot(r,y(x(k+1)*r));end§5.1求解本征值型的数学物理方程数学物理建模与计算机辅助设计二维本征值问题：圆形区域