基于Matlab语言数理方程解可视化01

1、基于Matlab语言的数理方程解的可视化 马 列 （玉溪师范学院理学院物理系08级物理1班 云南玉溪 653100） 指导老师：倪永勤 摘要：介绍Matlab语言在数学物理方程中的运用，并将数理方程的解用Matlab编程语言的形式表达出来，在采用Matlab软件中的绘图工具展示出来，使原本枯燥难懂的数学物理方程的解变得形象直观，有助于物理系学生的学习和教师的教学。关键词：Matlab 数理方程 可视化1、 引言数学物理方法是物理系及其他理工可学生的基础课程之一，其课程内容多而难，题目繁而杂，是一门公认的较难的课程。在教学过程中，主要是采用公式推导，大量的讲解和演算，所得结果往往是一个复杂的积分

2、或级数，甚至还有一些比较特殊的函数。尽管数学物理方法的习题通常都有明确的物理意义，可是怎样才能使学生从眼花缭乱的数学表达式中看出其中所表达的物理图像？这不仅是物理系学生的困惑，也是老师们的棘手问题。在解决这一问题的过程中，不同的教师或学生有不同的方法，有的可能会利用C语言，有的也会使用FORTRAN语言，但应用最广泛的还是Matlab语言。Matlab是一个集数学运算、图形处理、程序设计和系统建模的著名语言软件，具有功能强大、使用简单和编程简短易调试的优点。主要介绍Matlab在数理方程中的应用，并将数理方程的解可视化，让老师易教省事，学生易学易懂。 介绍了运用MATLAB求解数理方程中的几类

3、简单的方程：复变函数、波动方程、热传导方程，并用动画的形式演示了方程的解析解。2、 复变函数图形2.1 在Matlab指令窗口键入复变函数 下面举几个例子：2.1.1 2.1.2 2.1.3 2.1.4 在指令窗中，按下面的格式输入有关的指令，就能得出结果，在此，>>为指令提示符，ans后的内容是计算机输出的结果。>> sin(2+3i)ans = 9.1545 - 4.1689i>> log(-1) %log(x)函数求x的自然对数ans = 0 + 3.1416i>> z=2+3iz = 2.0000 + 3.0000i >> (

4、cosh(z)2-(sinh(z)2 %sinh、cosh函数分别实现双曲正弦和双曲余弦ans = 1.0000 - 0.0000i >> exp(2*z*i-3*sin(z)*i) %abs(x)函数求x的绝对值 ans = -9.0185e-010 +9.1311e-009i 2.2 Matlab画复变函数的图形为了能形象地表示复变函数的特性，可以画出它的图像，Matlab表现四位数据的方法是用三维空间坐标再加颜色，类似于地球仪用颜色表示海洋和高山一样，具体画法是以xy平面表示自变量所在的复平面，以z轴表示复变函数的实部，而用颜色来表示复变函数的虚部。再用Matlab语言画复变

5、函数的图形时，需要了解两个重要的函数，即CPLXGRID 构建一个极坐标的复数数据网格Z=CPLXMAP(m) 这是一个(m+1)*(2*m+1)的复数的极坐标下的数据网格，以下是几个画复变函数图形的例子2.2.1画的图形，所用程序如下>> z=cplxgrid(5);>> cplxmap(z,exp(2*z*i-3*sin(z)*i);>> xlabel('X');>> ylabel('Y');>> zlabel('exp(2*z*i-3*sin(z)*i)');>> co

6、lorbar('vert'); 图（2.2） >> title('z=exp(2*z*i-3*sin(z)*i)');结果如图（2.1） 2.2.2 画的图形，所用程序如下>> z=cplxgrid(5);>> w=log(z);>> for k=0:3 w=w+i*2*pi; surf(real(z),imag(z),imag(w),real(w); hold on 图（2.2） xlabel('X');ylabel('Y');zlabel('Z'); title(

7、'lnz')；colorbar('vert');end>> view(-75,30) 结果如图（2.2）从图（2.1）可以看出，自变量z的取值在水平面的单位圆内。X轴是实轴，Y轴是虚轴，变化范围都为-1+1。画函数时，是以坐标系的z轴表示实部，它所形成的曲面随自变量x、y和z的变化而变化,z有一个最大值2.5和一个最小值0，有两个低谷。图形右边的颜色条就是的坐标，不同的颜色标有不同的值，在曲面图中不难看出，随着X轴的变化，曲面颜色有绿色变为蓝色，其值域为-1+1。图（2.2）中，随着x和y的变化，z的值在单位圆内形成了螺旋上升的四层，即z有四个不同的

8、值，而不同的z值的自然对数的函数值却是相同的，即都为-1.60。3、 第一类数学物理方程：波动方程3.1 无限长的弦的自由振动泛定方程为： 因为无限长的线没有边界，即无边界条件，只有初始条件 ，问题的解是达朗贝尔公式： 当，而是，设初速度为 上述初始条件表明初速度的值为1，只分布于01之间。这是一个简单的可积分函数，由达朗贝尔公式得 由初始条件得 以上两式代表将函数 表示的波形向左和向右以的速度移动。利用上式，定义，通过以下程序作出解的模拟动画。在编写程序之前要先给定变量，在的范围内成立，再给定变量，在的范围内成立，超出这些区域的函数值按照上面的表达式由程序中接下来的四个语句来规定。该程序结构

9、如下：t=0:0.005:8; x=-10:0.1:10; a=1; X,T=meshgrid(x,t); u1=X+a*T; u1(find(u1<=0)=0; u1(find(u1>=1)=1;u2=X-a*T; u2(find(u2<=0)=0; u2(find(u2>=1)=1; jf=1/2/a*(u1-u2); h=plot(x,jf(1,:),'Linewidth',2) set(h,'Erasemode','xor');axis(-10 10 -1 1)hold onfor j=2:length(t) pa

10、use(0.01) set(h,'ydata',jf(j,:); drawnow;end 图（3.1）3.2、两端固定弦的自由振动1、两端固定弦的自由振动之一（初位移不为零初速度为零） 两端固定的均匀弦的自由振动定解问题是（参考书【数理方法】P180）(3.1)(3.2)(3.3)从上式可以看出，弦是有限长的，即为，且有两个固定端点，振动在这两个端点间来回传播，在遇到端点后又被反射回来，以此往复，从而形成驻波。各点的振幅X随地点x的变化而变化，即X是x的函数，将泛定方程分离变量，得： （3.4） 把（3.4）式代入弦振动方程（3.1）和边界条件（3.2），得 （3.5）在（3.

11、4）式中，自变量x只出现在X中，自变量t只出现在T中，把驻波的一般表示式代入弦振动方程（3.1）和边界条件（3.2）、(3.3）式，得： （3.6） （3.7）条件（3.6）的意义是：无论任何时刻,弦的两端点是不动的，即振幅为零，得：， （3.8）由（3.5）式，可以将时间函数和位移函数写在等号两边，既得： （3.9）两边相等是不可能的，而只能等于一个常数，即本征值。将（3.8）式分解为两个常微分方程： （3.10） （3.11） 求解（3.9）和（3.10）式可得泛定方程（3.1）的解： (n=1,2,3,) (3.12)其中 （n=1,2,3,） n 为正整数，即两端固定弦上的波数，对于每

12、一个不同的n 值，都有一种驻波，在用MATLAB编程时，取N=10，对应解的动画如图（3）所示，如果还想了解不同的驻波，可以再MATLAB的M文件中修改N的取值，这里只做N=10时的介绍。对于有限长的弦，如果将时间范围缩小来讨论，边界的影响是可以忽略的，即可以参照无限长弦的方法求解，为了研究方便，将其分为两种情况来考虑：（一）、，；（二）、，下面的程序介绍了第一种情况，设 取，得 将系数代入解的表达式中，级数中的每一项都是一个驻波。级数解有无穷项，为了画出级数解的图形，需要给定n值，以下程序画出了N=10的动画，使用了子程序wfun.m来计算级数中不同k的求和项，然后在主程序jxj中将他们加起

13、来。function jxjN=10t=0:0.005:2.0;x=0:0.001:1;ww=wfun(N,0);ymax=max(abs(ww);h=plot(x,ww,'linewidth',2);axis(0,1,-ymax,ymax)sy= ;for n=2:length(t) ww=wfun(N,t(n); set(h,'YData',ww); drawnow; pause(0.05) sy=sy,sum(ww);endfunction wtx=wfun(N,t)x=0:.001:1; a=1; wtx=0;for I=1:N if I=7 wtx=w

14、tx+0.05*(sin(pi*(7-I)*4/7)-sin(pi*(7-I)*3/7)/(7-I)/pi. -(sin(pi*(7+I)*4/7)-sin(pi*(7+I)*3/7)/(7+I)/pi)*cos(I*pi*a*t). .*sin(I*pi*x); else wtx=wtx+0.05/7*cos(I*pi*a*t).*sin(I*pi*x); endend 两端固定弦的振动之一（） 图（3.2）3.3、两端固定的弦振动问题之二（初位移为零初速度不为零） 与问题二比较，可以设初速度为： 上式的解析解是它的系数为：为了方便起见，不妨设, ,时间t与位置x的取值上一问题一致，同样取N

15、=50，使用子函数psi来计算级数中的各项，然后在主函数中调用，运用下面的程序可以作出解析解的动画，图（4）是动画的几个画面。function psiN=50;t=0:0.005:2.0;x=0:0.001:1;ww=psi1fun1(N,0);h=plot(x,ww,'Linewidth',3,'color','r');axis(0,1,-0.08,0.08)sy= ;for n=2:length(t) ww=psi1fun1(N,t(n); set(h,'ydata',ww); pause(0.05) title('两

16、端固定弦振动之二') %标题命名 drawnow; sy=sy,sum(ww);endfunction wtx=psi1fun1(N,t)x=0:0.001:1; a=1; wtx=0;for I=1:N BI=2./(I.*I.*pi*pi).*(cos(3.*I*pi/7)-cos(4*I*pi/7); wtx=wtx+BI*sin(I*pi*t)*sin(I*pi*x);end 图（3.3） 3.4、两端固定弦的振动问题三（有阻尼）弦振动如果遇到阻尼，振幅会不断减小，假定阻尼与速度成正比，定解问题是 首先将方程写成差分形式整理后得 在使用上式表示初始条件时，初速度为零可表示为在程

17、序中，初位移由表示，由初速度为零和初位移的表达式可求出，然后利用上式可以求出全部数值解。程序如下，取clearN=800;dx=0.01;dt=0.0005;c=400*dt*dt/dx/dx;g=5;u(1:100,1)=0;x=linspace(0,1,100)'u(2:66,1)=0.05/66*(2:66)'u(66:99,1)=0.05/34.*(100-(66:99)');u(2:99,2)=2*u(2:99,1)+c*(u(3:100,1)-2*u(2:99,1)+u(1:98,1);u(2:99,2)=1/(2+g*dt)*(u(2:99,2)+g*dt

18、.*u(2:99,1);h=plot(x,u(:,1),'Linewidth',3);set(h,'EraseMode','xor')axis(0,1,-0.05,0.05);for k=2:N set(h,'XData',x,'YData',u(:,2); title('有阻尼的两端固定弦振动') drawnow; u(2:99,3)=2*u(2:99,2)-u(2:99,1)+c*(u(3:100,2). -2*u(2:99,2)+u(1:98,2); u(2:99,3)=1/(1+g*dt)*

19、(u(2:99,3)+g*dt.*u(2:99,3); u(:,1)=u(:,2); u(:,2)=u(:,3); pause(0.01);end 图（3.4）3.5、偶极声源发射的稳恒声振动的速度势三维动画半径为r0的球面径向速度分布为 （3.5.1）试求解这球面所发射的稳恒声振动的速度势,设远小于声波的波长。由于（4.1）式中的即为，故本例称为偶极声源该例题没有初始条件，用球坐标，极点取在球心，定解问题是 （3.5.2） 边界条件的即。写成了，这要求约定计算的最后结果也应取其实部，查表可知：(1)问题与无关，（2）对的依赖关系为即，（3）边界条件中的时间因子，查得 （3.5.3）且，即，本

20、例只有这个唯一的值，所以无需叠加， （3.5.4） 为了确定系数，把（4.5）代入边界条件（4.3）， (3.5.5)即 （3.5.6）因，即很小，所以上式括号中第一项的绝对值远远大于其余两项，上式可简化为 由此，.于是得出答案 在远场区域即大的地方，运用渐近公式， 其实部为 这是振幅按减小的球面波。其对空间中方向的依赖也由描写，因此是偶极声源的声场。(参考书369)clearr0=0.2; v0=2; k=60; a=2;theta=linspace(0,2*pi,16);rho=0.1:0.1:2;Th,Rh=meshgrid(theta,rho);X,Y=pol2cart(Th,Rh);

21、rh=sqrt(X.2+Y.2);th=atan(Y./X);for t=0:0.001:0.03 u=real(v0/2*r03*(-1./rh.2+i*k./rh).*. cos(th).*exp(k*(rh+2*t)*i); surf(X,Y,u) view(-32,28) title('偶极声源发射的稳恒声振动的速度势') pause(0.5)end 图（3.5）将程序中的surf(X,Y,u) view(-32,28)改为contour(X,Y,u); axis(-4 4 -4 4); axis square; 偶极声源发射的稳恒声振动的速度势的等值线动画 图（3.6

22、）3.6、四极声源发射的稳恒声振动的速度势三维动画 半径为的球面径向速度分布为，是求解这球面所发射的稳恒声振动的速度势，设远小于声波的波长。运用球坐标，极点取在球心，定解问题如下 球面的边界条件是即为，在上式中写成，要求在计算结果中也只取实部。问题的解析解为 取其实部即为所求。在远场可以取渐进公式后再取实部，得 用下列程序画出解析解的动画，编写程序时分别取，；clearr0=2; v0=3;k=60; a=2;theta=linspace(0,2*pi,16);rho=0.1:0.1:2;Th,Rh=meshgrid(theta,rho);X,Y=pol2cart(Th,Rh);rh=sqrt

23、(X.2+Y.2);th=atan(Y./X);P2=1/4.*(3.*cos(2.*th)+1);for t=0:0.002:0.05 u=real(i*r04*( -1./rh.3+i*60./rh.2+. 3600/3./rh).*P2.*exp(i*k*(rh+2*t); surf(X,Y,u) view(-45,45) title('四极声源反射的稳恒声振动的速度势') pause(0.5)end 图（3.7）图（3.7）是动画中的几幅画面，表示由球面向外传播的球面波，有四个波源，故而朝四个方向传播，从远场近似解的表达式可以看出声波的传播具有明显的四极方向性，我们还可

24、以从图（3.8）的等值线中清楚地看到。 图（3.8）4、 第二类数学物理方程：输运方程4.1、有限长细杆的导热问题（齐次方程）有限长细杆在第一类边界条件下的热传导问题，其定界问题如下： ， （4.1） （4.2） （4.3） 其中取，即，取且 （4.4）所得解析解为 （4.5）在程序中用子函数w=rcdf(N,t(n)来计算各个分波，然后在主函数function jxj中调用，这里画出前100个分波的合成图形，程序如下：function jxjN=100;t=1e-5:0.00001:0.005; x=0:0.21:20;w=rcdf(N,t(1);h=plot(x,w,'linewi

25、dth',4);axis(0,20 ,0,1.5)for n=2:length(t) w=rcdf(N,t(n); set(h,'YData',w); drawnow; title('第一类边界条件下热传导问题') pause(0.01);endfunction u=rcdf(N,t)x=0:0.21:20; u=0;for k=1:2*N cht=2/k/pi*(cos(k*pi*10/20)-cos(k*pi*11/20)*sin(k*pi*x./20); u=u+cht*exp(-(k2*pi2*102/400*t); %a=10,l=20end图

26、（4.1）是动画中的两幅画面， 图（4.1）从图中可以看到，在开始时，温度分布是在附近的一个脉冲状分布，随着时间的变化，热量向两边传播，形成一个平缓的波包，如果时间足够长，最终杆上的温度会全部为零。4.2、球体内的热传导问题 匀质球，半径为，初始时刻球体温度均匀为。把它放入温度为的烘箱，是球面温度保持为。求解球内各处温度的变化情况。取球坐标系，极点在球心，定解问题是 首先把边界条件化为齐次，为此，移动温标零点，即，则定解问题为 （4.2.1）由此，上式的解析解为 （4.2.2）在解中只有一个空间变量，并将表示成，的函数，画成一个随时间变化的二维图形，结果见图（4.2）。等温线是一些同心圆，随着

27、时间的增加，这些同心圆由内到外的半径不断减小，函数值也不断增大，其变化范围为，在作图过程中分别取，在级数求和中去了前10项。U0=2; u0=0; a2=2; N=100;r=eps:0.05:1; theta=linspace(0,2*pi,16);t=0.1:0.001:0.2;RR,TT=meshgrid(r,t);R,TH=meshgrid(theta,r);X,Y=pol2cart(R,TH);for tt=1:100 un=0; for k=1:N unn=2*(U0-u0)*(-1).k.*sin(k.*pi.*(X.2+Y.2).(1/2).*. exp(-k2.*pi.2.*

28、a2.*t(tt)./(pi.*k.*(X.2+Y.2).(1/2); un=un+unn; endmesh(X,Y,un); title('匀质球内的热传导') xlabel('r') ylabel('t') axis(-1 1 -1 1 -0.4 0); pause(0.1)end 图（4.2） 5、 结束语 论文对数学物理方程的两种类型的方程的解进行了简单的可视化，通过长时间查阅大量的关于MATLAB的书籍，进一步了解了MATLAB 7.11.0的操作界面及其工作坏境，掌握了有关方面的编程技巧及绘图函数的运用，但MATLAB是一款博大精深的数学工程软件，要全面深入地了解或真正掌握，绝非一朝一夕之功。再加上数学物理方程繁多而复杂，要将其进行可视化，还需要深厚的数