基于MATLAB语言求偏微分方程

1、MATLAB语言课程论文基于MATLAB语言求偏微分方程 姓名：马兰 学号： 专业：通信工程 班级：2010级通信班 指导老师：汤全武 学院：物理电气信息学院 完成日期：2011年12月24日基于MATLAB语言求偏微分方程（ 马兰 2010级通信班）摘要数学物理方法中解偏微分方程难度非常大，因内容抽象，数学推导繁琐，平常学习起来感觉非常吃力。通过对MATLAB的学习，发现MATLAB是高性能的数值计算型数学类科技应用软件，具有优秀的数值计算功能和强大的数据可视化能力。应用MATLAB语言求解偏微分方程，一方面可以提高解题的速度，另一方面可将抽象的解和一些特殊函数以图形的形式显示出来，直观明了

2、，物理意义明确，为我们学习带来极大地方便。关键词偏微分方程 MATLAB语言 图形绘制1、 问题的提出 MATLAB是近几年传播速度最快，影响最大的数学类软件。应用MATLAB求解数学物理方法中的偏微分方程，使原来繁琐的手工计算变得简便，而且可以使一些解和特殊函数可以而且可以使一些解和特殊函数可以用图形的形式显示出来，形象，直观，便于理解。MATLAB强大的科学运算，灵活的程序设计，便捷的与其它程序语言接口的功能，显示其很强的优越性。 数学物理方法中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程。偏微分方程的求解是一个非常复杂的问题，除了几种极少数情况外，要求出它们精确解是很难

3、的。随着科学技术的发展，近几十年来其近似解法在理论上和方法上都有很大的进展，而且在各个领域内的应用也愈来愈广泛。下面我们就用MATLAB求解偏微分方程数值解作简单介绍。有限元方法可用来求解常微分方程(组)边值问题和偏微分方程（组），当求解的偏微分方程中空间变量的几何条件复杂时，有限元方法是解决问题的首选方法。而通过MATLAB解决偏微分显得比较容易。二、求解偏微分方程类型特征1、 椭圆型偏微分方程椭圆型偏微分方程的一般形式为 （1）其中：若，为的梯度，则其定义为 （2）散度的定义为 （3）这样，可以更明确地表示为 （4）若为常数，则进一步化简为 （5）其中，又称为Laplace算子。这样椭圆型

4、偏微分方程可以简单地写为 （6）2、抛物型偏微分方程抛物型偏微分方程的一般形式为 （7）根据上面叙述，若为常数，则该方程可以更简单地写为 （8）3、双曲型偏微分方程双曲型偏微分方程的一般形式为 （9）若为常数，则可以将该方程简化为 （10）三类方程的直接的区别在于对的导数的阶次。若对没有求导，可以理解为其值为常数，故称为椭圆型的。若取对时间的一阶导数，则与对的二阶导数直接构成了抛物线关系，故称为抛物型偏微分方程。若取对时间的二阶导数，称其为双曲型偏微分方程。4、特征值型偏微分方程特征值型偏微分方程为 （11）对常数该方程可以化简为 （12）该方程是椭圆型偏微分方程的特例。三、求解偏微分方程例证

5、 通过上面的了解学习，我知道了偏微分方程分类和利用MATLAB求解偏微分方程的方法，下面我就针对典型例题利用MATLAB求解，以实现检验与实践的目的。1、热传导方程的差分公式 热传导方程可以写成差分形式(右边取t时刻的值计算) （13）令,上式可写成显式差分公式稳定条件为 ，截断误差为 。 问题1：细杆传热问题定解问题： （14）其中且 （15）根据上面的公式,编写MATLAB程序如下clear all %清除x=0:20; %x从0到20取值t=0:0.01:1; %t从0到1取值a2=10; % 给a2赋值r=a2*0.01; %r的表达式u=zeros(21,101); %产生零矩阵uu

6、(10:11,1)=1; %对u中的数替换for j=1:100 %j从1到100循环u(2:20,j+1)=(1-2*r)*u(2:20,j)+r*( u(1:19,j)+ u(3:21,j);plot(u(:,j); %绘图函数调用axis(0 21 0 1); %设置坐标轴pause(0.1) %便于观察设置效果endsurf(u) %绘制三维表面图运行得到细杆传热问题结果如图1图1 细杆传热问题运行结果热传导方程还可以写成如下差分形式(右边取时刻的值计算) （16）引入与上面相同的足标,且令 （17）则得到隐式公式为 （18）则显式公式写成 （19）这个公式对任何步长都是恒稳定的。 含

7、时间的一维薛定谔方程 一维运动的粒子的含时间的薛定格方程可以化成如下形式的抛物型问题(设方程右边的系数为1), （20） 在使用MATLAB编程时,先计算出每一时刻的对应的程序中是把的系数表示为, 然后利用MATLAB的矩阵方程求解的功能求，再用它去求下一时刻的，程序在一个220点的格子上，定义解析形式的位势(方形势阱或势垒，Gauss势阱或势垒，位势台阶或抛物线势阱)，建立一个初始的Gauss波包或Lorentz波包，它具有规定的平均位置、动量和空间宽度，并在保持格子的端点上为零的边界条件下演化。用动画在每一时间步长上都显示概率密度和位势，以及波包在一指定点的左边和右边两部分每部分的总概率和

8、平均位置。 问题2 选位势为一个方形位垒，高度是0.18, 位于格程序使用的数据是，选位势为一个方形位垒，高度是0.18, 位于格点编号为105和115 的位置。初始的波包是一高斯波包，其平均波数和宽度为;中心位置x0在编号为40的格点上。时间步长为1，演化时间为130。动画演示了波包向势垒行进，在势垒上发生透射和反射，最后分为透射波和反射波二部分向不同的方向运动。编写MATLAB程序如下：clear all %清除NPTS=220; %设定格点范围sigma=10; %设置奇异值频域图k0=0.6; %设置平均波数x0=40; %设置中心位置time=130; %设定时间v(NPTS)=0;

9、 %设定初值v(105:115)=+0.18; %替换A=diag(-2+2i+v)+diag(ones(NPTS-1,1),1)+diag(ones(NPTS-1,1),-1); %建立对角阵运算PHI0(1)=0;PHI0(NPTS)=0; %设定初值for x=2:NPTS-1; %循环PHI0(x)=exp(0.5i*x)*exp(-(x-x0)2*log10(2)/sigma2); End %结束PHI(:,1)=PHI0.; %for y=2:timeCHI(:,y)=4i*(APHI(:,y-1);PHI(:,y)=CHI(:,y)-PHI(:,y-1);endmo=moviei

10、n(time);for j=1:timeplot(1:NPTS,abs(PHI(:,j).2/norm(PHI(:,j).2,1:NPTS,v/2);%输出mo(:,j)=getframe; endmovie(mo) %播放图NPTS=220; %设定格点范围通过上述程序得到运行结果如图2所示 图2 某一时刻透射波与反射波运动图像2、差分法与松弛法解椭圆型方程利用二阶导数的中心差分公式得到二维拉普拉斯方程 的差分公式为 （21） 记得显式 （22）稳定条件是. 问题3： 差分法解平面温度场,设正方形的一边温度为10度,其余各边的温度为零度,求稳定的温度场。按照上面的差分公式，任意设定内部各点的

11、温度为零，并假定当两次迭代计算的值误差小于0.01时，就停止计算。编写MATLAB程序如下：clear all %清除u=zeros(100,100); %产生100行100列的全零矩阵 u(100,:)=10; %取十个数uold=u+1; %定义表达式 unew=u; %把u赋给unewfor k=1:500 %循环语句if max(max(abs(u-uold)=0.01 %循环语句unew(2:99,2:99)=0.25*(u(3:100,2:99)+u(1:98,2:99)+.u(2:99,3:100)+u(2:99,1:98); %叠加语句uold=u; %赋值 u=unew; %

12、赋值End %结束endsurf(u) %绘制三维图形运行上述程序，结果如图3所示图3 差分法解平面温度运行结果 如果在差分公式中，随时将上一步算得的格点上的值替代旧值，并且每次算出的新值也替换成新值与旧值的“组合”，则得到下列松弛法的计算公式,其中。这个公式可以用变分原理去证明。 （23）问题4： 松弛法计算平面温度场,已知定解问题如下 (24)为了进行数值计算，取。编写MATLAB程序如下:clear all %清除omega=1.5; %定义常量值x=linspace(0,3,30); y=linspace(0,2,20); %调用linspace函数产生线性坐标phi(:,30)=si

13、n(3*pi/2*y); %定义表达式phi(20,:)=(sin(pi*x).*cos(pi/3*x); %运算函数for N=1:100for i=2:19for j=2:29ph=(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1); %运算表达式phi(i,j)=(1-omega)*phi(i,j)+0.25*omega*(ph); %定义表达式endendendcolormap(0.5,0.5,0.5); %设定颜色值figure(2), surfc(phi) %建立图形运行上述程序所的结果结果如图4 图4松弛法计算平面温度场结果图问题5： 单位圆

14、上的Poisson方程边值问题： （25）这一问题的精确解为 （26）编写MATLAB程序如下： clear all %清除 p,e,t=initmesh(circleg,hmax,1); %初始化网格 error=;err=1; %错误 while err0.001; %循环语句 p,e,t=refinemesh(circleg,p,e,t); %加密网格 u=assempde(circleg,p,e,t,1,0,1); %求解 exatc=(1-p(1,:).2-p(2,:).2)/4; %定义精度 err=norm(u-exact,inf); %错误 error=error err; E

15、nd %结束 pdemesh(p,e,t) %绘制网格图 pdesurf(p,t,u) %绘制求解图形 pdesurf(p,t,u-exact) %绘制误差图在MATLAB窗口输入 p 运行结果如下p = Columns 1 through 3 -1.0000 0.0000 1.0000 -0.0000 -1.0000 0 Columns 4 through 6 0.0000 -0.7071 0.7071 1.0000 -0.7071 -0.7071 Columns 142 through 144 0.4932 -0.1362 0.4255 -0.0428 -0.5976 0.5949输入 e

16、 运行结果如下e = Columns 1 through 3 1.0000 9.0000 10.0000 9.0000 10.0000 11.0000 0 0.1250 0.2500 0.1250 0.2500 0.3750 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0 0 0 Columns 31 through 32 31.0000 32.0000 32.0000 1.0000 0.7500 0.8750 0.8750 1.0000 4.0000 4.0000 1.0000 1.0000 0 0输入 t 运行结果如下t = Columns 1 th

17、rough 6 32 14 20 26 29 17 1 2 3 4 8 6 89 97 81 98 84 92 1 1 1 1 1 1 Columns 253 through 254 61 43 102 140 144 144 1 1输入u 运行结果如下u = 0 0 0 0.1712 0.0783 0.0783 0.1693 0.1706 0 0 0 0 0.1005 0.0945 0.1313最终解决单位圆上的Poisson方程边值问题所得的运行结果如图5，6，7所示 图5 运算网格图图6 运算误差图图7 运算结果图通过以上例解决单位圆上的Poisson方程边值问题，我们可以清晰地看到用m

18、atlab解决该类偏微分问题是多么的简单，大量节省了计算绘图所带来的困扰与麻烦。四、结论上述所解答的实例问题，是我们平时在偏微分方程问题的时候所遇到非常的麻烦几类，但通过运用MATLAB解决此类问题，就显得非常的简单了。从以上利用MATLAB语言对偏微分方程求解，可见偏微分方程是一门实用性很强的学科。对于理论研究和实际应用中提出的偏微分方程问题，由于其边界和边界条件复杂等原因，寻求解析解是非常困难甚至不可能的，利用计算机研究相应问题的数值解法是十分必要的。而借助于MATLAB这个工具，我们可以从繁琐，共性的求解步骤中解脱出来而专注于问题的核心即问题的描述，定义及简化，边界条件的确定，求解方法和精度控制的选择等，大大提高了求解效率。在应用MATLAB解决偏微分问题是要特别注意边界的处理，及相关表达式的格式定义，还要注意每个函数所代表的含义，当我们不清楚函数定义时，可采用在MATLAB窗口输入