1、基于 MATLAB 的偏微分方程差分解法学院：核工程与地球物理学院专业：勘查技术与工程班级：1120203学号：姓名:2014/6/111在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。近些年来，求解偏微分方程的数值方法取得进展，特别是有限差分区域分解算法，此类算法的特点是在内边界处设计不同于整体的格式， 将全局的隐式计算化为局部的分段隐式计算。使人从感觉上认为这样得到的解
2、会比全局隐式得到的解的精度差，但大量的数值实验表明事实正好相反，用区域分解算法求得的解的精度更好。差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分
3、方程定解问题一般需要解决以下问题：（i）选取网格；（ii）对微分方程及定解条件选择差分近似，列出差分格式；（iii ）求解差分格式；（iv）讨论差分格式解对于微分方程解的收敛性及误差估计。下面对偏微分方程具体例题的差分解法作一简要的介绍。1 双曲型方程中波动方程的有限差分解法。1.1 双曲型的差分方程通过建立网格并求解中心差分方程结果为： 22,1,1,1()()2,31ij ijijijijuruuin 。其中为了保证上式的稳定性，必须使 /rckh。1. 2 初始值通过联立初始值及边界条件可以得到 2 2311(,)()()(1)iiiiickuxkfgffOhkh代入 ，可简化并得到一个
4、改进的对行2的近似值差分方程：/rch22,2 1()()(2)iiiiirrfkf1.3 双曲型方程中波动方程例题的差分解法结果及程序。题： ，其中 ，边界条件为：(,)4(,)txut0x且 t0.5(1,).(,)sin,100tutf xxg且设 0.2,.1,.hkr2解：第一步：通过联立(1)、(2)编写MATLAB程序如下：%二维双曲型偏微分方程，使用DAlembert方法function U=hyperbolic(a,b,c,n,m)% a为x的取值范围% b为t的取值范围% c为系数% n为x方向上的节点数% m为t上的节点数h=a/(n-1); % x方向上的步长k=b/(
5、m-1); % t上的步长r=c*k/h;r2=r2;r22=r2/2;s1=1-r2;s2=2-2*r2;U=zeros(n,m);for i=2:(n-1);U(i,1)=sin(pi*h*(i-1);U(i,2)=s1.*sin(pi*h*(i-1)+k*0 .+r22.*(sin(pi*h.*(i-2)+sin(pi*h.*(i); %P116（13）end%差分方程for j=3:m;for i=2:(n-1);U(i,j)=s2*U(i,j-1)+r2*(U(i+1,j-1)+U(i-1,j-1)-U(i,j-2);%P115(7)endend U=U;%画图figure(1);
6、surf(U);figure(2);contour(U,40);第二步：输入数值并计算a=1;b=0.5;c=4n=11m=11执行hyperbolic(1,0.5,4,11,11)；第三步：得出结果并画图入下1.等值线结果图32.三位结果图1.4 通过MATLAB语言提供了pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)，具体程序见附录得出的结果为：1.等值线结果图42.三维结果图1.5 结果对比通过编写MATLAB的差分方程程序求取结果和 MATLAB自带函数求取结果进行对比，发现这两种方法求得到的
7、结果是非常理想的。2 抛物线方程中热传导方程的有限差分解法。2.1 抛物线方程的差分方程通过建立网格并求解显示前向差分方程结果为： ,1,-1,+,(2)()(3)ijijijijuruu其中为了保证上式前向差分方程稳定性，当且仅当 满足 时。这意味着步长 必r102k须满足 。2khc2.2 抛物线方程中热传导方程例题的差分解法结果及程序。题： 其中 初始条件为 其,(,)tuxt01,0.,xt,0()12,uxfx中 边界条件为：01t12,.tucxt且且解：第一步，分析并带入（3）并编写MATLAB求解程序如下：function U=forwdif(c1,c2,a,b,c,n,m)c
8、lch=a/(n-1);k=b/(m-1);r=(c2*k)/(h2);s=1-2*r;U=zeros(n,m);U(1,1:m)=c1;U(n,1:m)=c2;for i=2:n-15U(I,1)=1-abs(2*(i-1)*h-1);% U(I,1)=4*(i-1)*h-4*(i-1)*h)2;endfor j=2:mfor i=2:n-1U(I,j)=s*U(I,j-1)+r*(U(i-1,j-1)+U(i+1,j-1);endendU=U;figure(1); surf(U);figure(2);contour(U,30);第二步，代入初始条件以及边界条件：c1=0；c2=0；a=1；
9、b=0.5；c=1；n=6；m=11；执行forwdif(0,0,1,0.1,1,11,11)；第三步：得出结果并画图入下1.等值线结果图2.三位结果图62.3 通过MATLAB语言提供了pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)，具体程序见附录得出的结果为：1.等值线结果图2.三维结果图2.4 结果对比通过编写MATLAB的差分方程程序求取结果和 MATLAB自带函数求取结果进行对比，发现这两种方法求得到的结果是非常相似的，差距不大证明程序编写是成功的。3 椭圆型方程的有限差分解法。3.1 建立
10、线型方程组 11,(,)21(), ()jjiinjnjimiuxyjiji在 左 边在 底 边在 右 边在 顶 边3. 2 导数边界条件Neumann边界条件确定了 边的法线的方向导数。这里使用零法线导数条件：(,)uxy7(,)0uxyN对于热传导而言，这表示边是热绝缘的而且经过边的热通量为零，从而得到： (,)(,)njxnj得出点 的Laplace差分方程为：(,)njxy1,1,0(4)njjnjjuu利用差分方程： 1,(,)(5)2njjxniyh在（4）上式使用（5）式这个近似值，这结果为： 1,+1,-,40(4)njjnjjuu这个公式讲函数 联系起来。1,njjj， 和则
11、Neumann计算空格样板的4种情况如下所示： ,21,1,2,1,1,-40()(5-()iiiimmjjjjnnuu底 部顶 部左 部右 部3. 3 迭代方法根据（4）式和（5）式以如下形式进行迭代处理： ,(6)ijijijur式中： 1,1, 4(7)ijijijijijij ur这里 2.xnm且通过迭代公式（6）右边对所有的内部结点连续执行迭代，直到该式右边余项 “减少到,ijr零”（即， ）。可以利用逐次超松弛法（SOR）提高,|2121ijrinj且所有余项 减少到零的收敛速度。其中逐次超松弛法使用迭代公式：,ij81,1, 4(8)ijijijijijijijijijuuuu
12、r 这里参数 位于 范围内。对所有网格点应用上式（8）直到 。其中：12 ,ijr24(9)coscs11nm3.4 椭圆型方程中波动方程例题的差分解法结果及程序。题：用迭代法求解在区域 Lpalace方程的 近,:04,rxyx20似解，这里边界值为： ,2;,1804u解：第一步：通过联立差分方程及迭代方程编写MATLAB程序如下：function U=dirich(a,b,h,tol,maxl)%设DX=DY=h,且存在m，是的a=nh和b=mhn=fix(a/h)+1;m=fix(b/h)+1;ave=(a*(20+180)+b*(80+0)/(2*a+2*b); U=ave*ones
13、(n,m);U(1,1:m)=80;U(n,1:m)=0;U(1:n,1)=20;U(1:n,m)=180;U(1,1)=(U(1,2)+U(2,1)/2;U(1,m)=(U(1,m-1)+U(2,m)/2;U(n,1)=(U(n-1,1)+U(n,2)/2;U(n,m)=(U(n-1,m)+U(n,m-1)/2;w=4/(2+sqrt(4-(cos(pi/(n-1)+cos(pi/(m-1)2);err=1;cnt=0;while(errtol)for j=2:m-1for i=2:n-1relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j)
14、/4;U(i,j)=U(i,j)+relx;if(err=abs(relx)err=abs(relx);9endendendcnt=cnt+1;relx;endU=flipud(U);figure(1); surf(U);figure(2);contour(U,30);第二步：输入数值并计算a=4; b=4; h=0.5; tol=0.1; maxl=20;执行dirich(4,4,0.5,0.1,20)；第三步：得出结果并画图入下1.等值线结果图2.三位结果图103 附录双曲型及抛物线方程MATLAB语言提供的pdepe()函数程序。%写主函数clcx=linspace(0,1,50);t=
15、linspace(0,0.1,50);m=0;sol=pdepe(m,pdefun,pdeic,pdebc,x,t);u=sol(:,:,1);figure(numbertitle,off,name,PDE Demoby Matlabsky)surf(x,t,u)title(The Solution of u_1)xlabel(X)ylabel(T)zlabel(U)figure(2);contour(u,40);% 目标PDE函数function c,f,s=pdefun (x,t,u,du)%p141-4c=1;f=du;s=0;% %p139-4% c=4;% f=du;% s=0;% 边界条件函数function pa,qa,pb,qb=pdebc(xa,ua,xb,ub,t)%a表示下边界，b表示上边界%p141-4pa=ua;qa=0;pb=ub;qb=0;% %p139-4% pa=ua;% qa=0;% pb=ub;% qb=0;% 初值条件函数function u0=pdeic(x)%p141-4u0=1-abs(2*x-1);11% %p139-4% u0=sin(pi*x);
