matlab 差分解微分,基于MATLAB的偏微分方程差分解法

基于MATLAB的偏微分方程差分解法

基于 MATLAB 的偏微分方程差分解法学院：核工程与地球物理学院专业：勘查技术与工程班级：1120203学号：姓名:2014/6/111在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。近些年来，求解偏微分方程的数值方法取得进展，特别是有限差分区域分解算法，此类算法的特点是在内边界处设计不同于整体的格式， 将全局的隐式计算化为局部的分段隐式计算。使人从感觉上认为这样得到的解会比全局隐式得到的解的精度差，但大量的数值实验表明事实正好相反，用区域分解算法求得的解的精度更好。差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点(网格点)集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解(数值解)。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：(i)选取网格；(ii)对微分方程及定解条件选择差分近似，列出差分格式；(iii )求解差分格式；(iv)讨论差分格式解对于微分方程解的收敛性及误差估计。下面对偏微分方程具体例题的差分解法作一简要的介绍。§1 双曲型方程中波动方程的有限差分解法。1.1 双曲型的差分方程通过建立网格并求解中心差分方程结果为： 22,1,1,,,1()()2,31ij ijijijijuruuinL。其中为了保证上式的稳定性，必须使 /rckh。1. 2 初始值通过联立初始值及边界条件可以得到 2 2311(,)()()(1)iiiiickuxkfgffOhkh代入 ，可简化并得到一个改进的对行2的近似值差分方程：/rch22,2 1()()(2)iiiiirrfkf1.3 双曲型方程中波动方程例题的差分解法结果及程序。题： ，其中 ，边界条件为：(,)4(,)txut0 x且 t0.5(1,).(,)sin,100tutf xxg且设 0.2,.1,.hkr2解：第一步：通过联立(1)、(2)编写MATLAB程序如下：%二维双曲型偏微分方程，使用D Alembert方法function U=hyperbolic(a,b,c,n,m)% a为x的取值范围% b为t的取值范围% c为系数% n为x方向上的节点数% m为t上的节点数h=a/(n-1); % x方向上的步长k=b/(m-1); % t上的步长r=c*k/h;r2=r^2;r22=r^2/2;s1=1-r^2;s2=2-2*r^2;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); 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自带函数求取结果进行对比，发现这两种方法求得到的结果是非常理想的。§2 抛物线方程中热传导方程的有限差分解法。2.1 抛物线方程的差分方程通过建立网格并求解显示前向差分方程结果为： ,1,-1,+,(2)()(3)ijijijijuruu其中为了保证上式前向差分方程稳定性，当且仅当 满足 时。这意味着步长 必r102k须满足 。2khc2.2 抛物线方程中热传导方程例题的差分解法结果及程序。题： 其中 初始条件为 其,(,)tuxt01,0.,xt,0()12,uxfx中 边界条件为：01t12,,.tucxt且且解：第一步，分析并带入(3)并编写MATLAB求解程序如下：function U=forwdif(c1,c2,a,b,c,n,m)clch=a/(n-1);k=b/(m-1);r=(c^2*k)/(h^2);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；b=0.5；c=1；n=6；m=11；执行forwdif(0,0,1,0.1,1,11,11)；第三步：得出结果并画图入下1.等值线结果图2.三位结果图6§2.3 通过MATLAB语言提供了pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)，具体程序见附录得出的结果为：1.等值线结果图2.三维结果图2.4 结果对比通过编写MATLAB的差分方程程序求取结果和 MATLAB自带函数求取结果进行对比，发现这两种方法求得到的结果是非常相似的，差距不大证明程序编写是成功的。§3 椭圆型方程的有限差分解法。3.1 建立线型方程组 11,,,(,)21(), ()jjiinjnjimiuxyjiji在 左 边在 底 边在 右 边在 顶 边3. 2 导数边界条件Neumann边界条件确定了 边的法线的方向导数。这里使用零法线导数条件：(,)uxy7(,)0uxyN对于热传导而言，这表示边是热绝缘的而且经过边的热通量为零，从而得到： (,)(,)njxnj得出点 的Laplace差分方程为：(,)njxy1,,,1,0(4)njjnjjuu