1、数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计上机练习专题上机练习专题5：数学物理方程的计：数学物理方程的计算机求解和可视化算机求解和可视化对代数表达式进行可视化对代数表达式进行可视化 二维本征值问题二维本征值问题 矩形区域的本征模与本征振动矩形区域的本征模与本征振动边长为边长为b和和c的的四周固定的矩形膜振动的本征值问题为的的四周固定的矩形膜振动的本征值问题为xxyyuuu采用分离变量法可以得到本征模和本征值为采用分离变量法可以得到本征模和本征值为 22sinsinnmmnn xm yXx Yycbmnbc2ttxxyyUaUU2对代数表达式进行可视化对代数表达式进行可视化 二维本
2、征值问题 绘制前4个本征函数的图形%P70_1.ma=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);me
3、sh(X,Y,w12);subplot(2,2,3);mesh(X,Y,w21);subplot(2,2,4);mesh(X,Y,w22);22sinsinnmmnn xm yXx Yycbmnbc3对代数表达式进行可视化对代数表达式进行可视化 二维本征值问题 动态可视化：设具有时间因子%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);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
4、*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)sinmnmnTtat4PDE工具箱求解数学物理方程工具箱求解数学物理方程 PDE：Partial Differential Equation PDE工具箱功能工具箱功能 设定二维定解区域、边界条件以及方程的形式和系数设定二维定解
5、区域、边界条件以及方程的形式和系数 用有限元法用有限元法(FEM)生成网格、方程离散化并求出数值解生成网格、方程离散化并求出数值解 对所求解的可视化，图像和动画对所求解的可视化，图像和动画 PDE工具箱可以求解的工具箱可以求解的4种偏微分方程种偏微分方程 椭圆型方程椭圆型方程 抛物型方程抛物型方程 双曲形方程双曲形方程 特征值方程特征值方程c uauftduc uaufttduc uaufc uaudu5PDE工具箱求解数学物理方程工具箱求解数学物理方程 PDE工具箱给定的两种边界条件工具箱给定的两种边界条件 Dilichlet条件条件 Neumann条件条件 用用PDE工具箱求解二维本征值问
6、题工具箱求解二维本征值问题hurc uqugn6PDE工具箱求解数学物理方程工具箱求解数学物理方程 求解二维本征值问题7PDE工具箱求解数学物理方程工具箱求解数学物理方程1.求解定解问题0(0,0)3(0,)0,(,)sin3(,0)0,(,)sincosxxyyuuxaybyu xyu xa ybxxu x yu x ybaa取1,1,1ab。8PDE工具箱求解数学物理方程工具箱求解数学物理方程2.求矩形薄膜的横向振动。定解问题为求矩形薄膜的横向振动。定解问题为2()(0,0,0)(0,)0,(,)0(,0,)0,(,)0(,0)()sin,(,0)0ttxxyytua uuxcybtu x
7、y tu xa y tu x ytu x yb tyu x y tAx xcu x y tb取1,1,2,1abcA9差分法求解数学物理方程差分法求解数学物理方程 导数的差分公式导数的差分公式 在在x附近将函数附近将函数f(x)展开成泰勒公式展开成泰勒公式 前差公式前差公式 后差公式后差公式 二阶导数的差分公式二阶导数的差分公式(中心差分公式中心差分公式)f xhf xfx h f xhf xfxh f xf xhfxh 22fxhfxf xhf xf xhfxhh10差分法求解数学物理方程差分法求解数学物理方程 用数值方法求解微分方程的时候步骤：用数值方法求解微分方程的时候步骤：(1)变量空
8、间网格化变量空间网格化(2)微分方程变差分方程微分方程变差分方程(3)定解条件设置定解条件设置(4)整理得到递推公式整理得到递推公式11差分法求解数学物理方程差分法求解数学物理方程 一维波动问题 两端固定的弦的振动问题 设定解条件为 解的代数表达式为 20,0,0,0,0ttxxtua uutu l tu xxuxx 734sin,7700llxxxxl其它17,cossinsin143143sin 7sin 7sin 7sin 777777770,0nnnnnn atn atn xu x tABlllAnnnnnnnAB12差分法求解数学物理方程差分法求解数学物理方程 两端固定的弦的振动问题
9、的求解步骤两端固定的弦的振动问题的求解步骤(1)变量空间网格化变量空间网格化(2)微分方程变差分方程微分方程变差分方程1,2,1,2,xi xiMtj tjN ,1,121,1,2,1,11,1,2222222i ji ji jttiji jijxxi ji ji jiji jijuuuutuuuuxuuuuuuatx中心差分公式中心差分公式13差分法求解数学物理方程差分法求解数学物理方程 两端固定的弦的振动问题的求解步骤(3)定解条件设置 1,1,2,1,2,10,0,0,0,0;00jM jiiiiiiitiutu l tuuuuxuuxuuxttxxuxx ,1,2,1,2,11,0,0
10、,0,0,000;jtM jiiiiiiiiutu l tu xxuxuuxuuxutxuuxt ,1,2,1,2,11,0,0,0,0,0;00tjMiiiiiijiiuuxuuxuuxtuuttu l tu xxuxx ,2,11,2,1,10,0,0,0,000;iiiiiijMitjiutu l tu xxuuuuuuxxxxuuxtt 1,1,2,1,2,10,0,0,0,0;00jM jiiiiiiitiutu l tuuuuxuuxtuxxuxuxtx 1,1,2,1,2,10,0,0,0,0;00jM jiiiiiiitiutu l tuuuuxuuxtuxxuxuxtx思考：
11、思考：1.求时间分量的递推公式需要几个初始值？求时间分量的递推公式需要几个初始值？2.对坐标分量呢对坐标分量呢?3.两个初始值时，另一个初始值如何处理？两个初始值时，另一个初始值如何处理？14差分法求解数学物理方程差分法求解数学物理方程 两端固定的弦的振动问题的求解步骤(4)整理得到递推公式,11,1,21222,i jiji jiji ji jtcuac uuuuux ,2,1,1100jM jiiiiiuuuutuxx,21,11,1,112 12iiiiuc uuc u15差分法求解数学物理方程差分法求解数学物理方程%P164_1.mN=4001;dx=0.0024;dt=0.0005;
12、c=dt*dt/dx/dx;x=linspace(0,1,420);u(1:420,1)=0;u(181:240,1)=0.05*sin(pi*x(181:240)*7);u(2:419,2)=u(2:419,1)+c/2*(u(3:420,1)-2*u(2:419,1)+u(1:418,1);h=plot(x,u(:,1),linewidth,3);axis(0,1,-0.05,0.05);set(h,erasemode,xor,markersize,18);for k=2:N set(h,XData,x,YData,u(:,2);drawnow;u(2:419,3)=2*u(2:419,2
13、)-u(2:419,1)+c*(u(3:420,2)-2*u(2:419,2)+u(1:418,2);u(2:419,1)=u(2:419,2);u(2:419,2)=u(2:419,3);end 两端固定的弦的振动问题 程序的编制关键在于：关键在于：(1)正确的递推公式正确的递推公式(2)正确设定定解条件正确设定定解条件16本专题小结本专题小结 直接对代数表达式求解并进行可视化直接对代数表达式求解并进行可视化 利用差分法求解并进行可视化利用差分法求解并进行可视化 利用利用Matlab的的PDE工具箱求解并进行可视化工具箱求解并进行可视化17课堂练习课堂练习 1.将将P85例题例题4.1.1用偏微分方程工具箱来求解。用偏微分方程工具箱来求解。2.将将P146例题例题5.1.4关于细杆的热传导问题用偏微分关于细杆的热传导问题用偏微分方程工具箱来求解。方程工具箱来求解。18
