下载文档原格式
本科毕业论文（设计）题目二维双曲线型方程的交替方向隐格式解法院（系）数学系专业数学及应用数学学生姓名周玲玲学号 10022156指导教师陈淼超职称讲师论文字数 9500完成日期: 2014年6月8日巢湖学院本科毕业论文(设计)诚信承诺书本人郑重声明：所呈交的本科毕业论文(设计)，是本人在导师的指导下，独立进行研究工作所取得的成果。
除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。
对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。
本人完全意识到本声明的法律结果由本人承担。
本人签名：日期：巢湖学院本科毕业论文 (设计)使用授权说明本人完全了解巢湖学院有关收集、保留和使用毕业论文 (设计)的规定，即：本科生在校期间进行毕业论文(设计)工作的知识产权单位属巢湖学院。
学校根据需要，有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许毕业论文 (设计)被查阅和借阅；学校可以将毕业论文(设计)的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编毕业，并且本人电子文档和纸质论文的内容相一致。
保密的毕业论文(设计)在解密后遵守此规定。
本人签名：日期：导师签名：日期：二维双曲型方程的交替方向隐格式解法摘要在解决偏微分方程中二维双曲线型方程的问题求解初边值问题时，显示格式的稳定性是有条件的，并且多维的稳定性条件更严，为得到稳定性好的格式，隐式格式受到重视。
用隐式格式求解二维问题得到的线性方程组其系数矩阵为宽带状，因此求解不甚便利，采用交替方向隐式（ADI ）格式可以避免此问题。
本文以基本概念和基本方法为主，同时结合算例实现算法。
第一部分介绍二维双曲线型方程的交替方向隐格式解法的基本概念，引入本文的研究对象——二维波动方程： ，x R ∈，],0(T t ∈第二部分介绍上述方程的二维双曲线型方程的交替方向隐格式及这种格式的存在性、收敛性及稳定性。
GAGGAGAGGAFFFFAFAF偏微分數值解法實驗報告實驗名稱：六點對稱格式，ADI 法，預校法，LOD 法解二維拋物線方程的初值問題實驗成員： 吳興 楊敏 姚榮華 于瀟龍 余凡 鄭永亮實驗日期：2013年5月17日 指導老師：張建松一、 實驗內容用六點對稱格式，ADI 法，預校法和LOD 法求解二維拋物線方程的初值問題：21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sin cos .xx yy y y u u u x y G t t u y t u y t y t u x t u x t x t u x y x y ππ∂⎧=+∈=⨯>⎪∂⎪⎪==<<>⎨⎪==<<>⎪=⎪⎩已知（精確解為：2(,,)sin cos exp()8u x y t x y t πππ=-）設(0,1,,),(0,1,,),(0,1,,)j k n x jh j J y kh k K t n n N τ======差分解為GAGGAGAGGAFFFFAFAF,nj k u ，則邊值條件為：0,,,0,1,1,0,0,1,,,,0,1,,n n k J k nn n nj j j K j K u u k K u u u u j J-⎧===⎪⎨===⎪⎩初值條件為：0,sin cos j k j k u x y ππ=取空間步長12140h h h ===，時間步長11600τ=網比。
1: ADI 法：由第n 层到第n+1层计算分成两步：先先第n 層到n+1/2層，對uxx 用向后差分逼近，對uyy 用向前差分逼近，對uyy 用向后差分逼近，于是得到了如下格式：11112222,,1,,1,,1,,1221222,,2-22=21()n n n n n n n nj kj kj kj k j kj k j k j k n nx j k y j k hh hτδδ+++++-+-+-+-+=+uu uuuu u u （+） （1）u u1111111222,,1,,1,,1,,12212212,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hh hτδδ++++++++-+-++-+-+=+u u uuuu u u （+） （2)u u其中j,k=1,2,…,M-1,n=0,1,2,…,上标n+1/2表示在t=tn+1/2+(n+1/2)τ取值。
二维抛物型方程的交替方向隐格式二维抛物型方程的交替方向隐格式是一种数值解法，用于求解二维抛物型方程的数值解。
这种解法将二维问题分解为两个一维问题，并采用隐式差分方法来求解。
具体来说，二维抛物型方程可以表示为：u_t = a^2 u_{xx} + b^2 u_{yy} + f(x,y,u)其中，u是待求解的函数，t是时间，a和b是两个不同的参数，f是右侧的非线性函数。
为了求解这个问题，我们可以采用交替方向隐式差分方法，将问题分解为两个一维问题：1、在x方向上，从左到右扫描每一行数据，更新每个点的u值。
这个过程可以使用隐式差分方法来实现：u^[i,j]^(n+1) = u^[i,j]^(n) + dt a^2 (u^[i+1,j]^(n) - 2u^[i,j]^(n) + u^[i-1,j]^(n)) / h^2 + dt f^[i,j]^(n) 其中，u^[i,j]^(n)表示第n个时间步中，位置为(x[i],y[j])的点的u值，h是x方向上的步长，dt是时间步长，f^[i,j]^(n)表示第n个时间步中，位置为(x[i],y[j])的点上的非线性函数f的值。
2. 在y方向上，从上到下扫描每一列数据，更新每个点的u值。
这个过程也可以使用隐式差分方法来实现：u^[i,j]^(n+1) = u^[i,j]^(n) + dt b^2 (u^[i,j+1]^(n) - 2u^[i,j]^(n) + u^[i,j-1]^(n)) / h^2 + dt f^[i,j]^(n) 其中，u^[i,j]^(n)表示第n个时间步中，位置为(x[i],y[j])的点的u值，h是y方向上的步长，dt是时间步长，f^[i,j]^(n)表示第n个时间步中，位置为(x[i],y[j])的点上的非线性函数f的值。
通过交替更新每个点的u值，我们可以逐步逼近方程的数值解。
这种解法具有二阶精度和稳定性，可以应用于多种二维抛物型方程的问题。
关于解二维椭圆型高精度差分方程的交替方向迭代法交替方向迭代法（alternating direction iterative method，简称ADI法）是一种常用于解二维椭圆型高精度差分方程的数值方法。
它的基本思想是将二维问题拆分为两个一维问题，然后依次迭代求解，从而达到解二维问题的目的。
一、AD方法的原理二维椭圆型方程可以表示为：>(A1u(i+1,j)+B1u(i,j+1)-(2A1+B1+C1)u(i,j)+B1u(i,j-1)+A1u(i-1,j))/Δx^2+(D1Δy^2)u(i,j)=f(i,j)，式中i,j分别为空间坐标，f(i,j)为已知的外力场，D1为系数。
对于前一半时间步，我们考虑j方向上的一维问题，可将上述方程改写为：>(A1+B1)u(i+1,j)+(C1-2A1-B1)u(i,j)+A1u(i-1,j)=f1(i,j)，式中f1(i,j)为已知。
对于后一半时间步，我们考虑i方向上的一维问题，可将上述方程改写为：>(A2+B2)u(i,j+1)+(C2-2A2-B2)u(i,j)+A2u(i,j-1)=f2(i,j)，式中f2(i,j)为已知。
通过交替迭代这两个方程，逐步逼近真实解。
二、AD方法的步骤1.选取初始解u0(i,j)，以及参数Δx，Δy和迭代终止准则ε。
2.开始迭代，设定初始误差e0=ε+1，n=0。
3. 使用前一半时间步的方程求解每个i所对应的j方向的一维问题，得到一个临时解ut(i,j)。
4.使用后一半时间步的方程求解每个j所对应的i方向的一维问题，得到一个新的解u(n+1)(i,j)。
5. 计算误差en=，u(n+1)-u(n)，2，n=n+16. 如果en ≤ ε，迭代结束，得到近似解u(n+1)；否则返回第3步。
7.算法结束。
三、AD方法的优势1.AD方法是一种分离的迭代方法，当网格较大时其计算量相对较小。
2.AD方法是一种隐式迭代方法，对稳定性和收敛性有较好的保证。
ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种）理论部分参看孙志忠：偏微分方程数值解法注意：1.最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。
2.中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。
3.下面有三个程序4.具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。
Matlab程序：1.function [u u0 p e x y t]=ADI1(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u0(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))]; u0([1 m2+1],j,k)=[exp(0.5*x(j)-t0(k)) exp(0.5*(1+x(j))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r2*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+...h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r2*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend2.function [u p e x y t]=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(D'Yakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u'(i,j,k)(引入的过渡层),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定义t0是为了f(x,y,t)~=0的情况% for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));%编程时-t0(k)写成了+t0(k),导致错误;endendend%初始条件for i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endend%边界条件for k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))];endendr=h2/(h1*h1);r4=1+r;r5=r/2;for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=r4*u(i,[1 m1+1],k+1)-r5*(u(i-1,[1 m1+1],... k+1)+u(i+1,[1 m1+1],k+1));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=2*r4*ones(1,m1-1);d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+...u(i,3,k))+r2*u(i,2,k)+r3*(u(i-1,1,k)+...u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k))+2*h2*f(i,2,k);for l=2:m1-2d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+...u(i,l+2,k))+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+...u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k))+2*h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+...u(i+1,m1,k)+u(i,m1+1,k))+r2*u(i,m1,k)+...r3*(u(i-1,m1-1,k)+...u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k))+2*h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);end%回代过程%u0(i,m1,k)=d(m1-1)/b(m1-1);for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=2*r4*ones(1,m2-1);d(1)=r*u(1,j,k+1)+2*u0(2,j,k);for l=2:m2-2d(l)=2*u0(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend3.function [u u0 p e x y t]=ADI5(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，未截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u1(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=u1(i,[1 m1+1],k)-r2*(u(i-1,[1 m1+1],k+1)-...2*u(i,[1 m1+1],k+1)+u(i+1,[1 m1+1],k+1)-u(i-1,[1 m1+1],k)+...2*u(i,[1 m1+1],k)-u(i+1,[1 m1+1],k));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendfor k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组%for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r3*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+... h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r3*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend[up e x y t]=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001)) t=1的误差曲面下面是三种方法的误差比较：1.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2）代替）（t=1时的误差）2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) surf(x,y,e(:,:,11))(表示t=1时的误差）下面是相关数据：1: [u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00040947 0.00025182 0.00019077 0.00017112 0.000176040 0.00057359 0.00042971 0.00035402 0.00032565 0.000336280 0.00066236 0.00054689 0.00047408 0.00044596 0.000462670 0.00072152 0.00062001 0.00055081 0.00052442 0.000545530 0.00076164 0.0006576 0.00058522 0.00055732 0.000579840 0.00078336 0.00065993 0.00057557 0.00054161 0.000562090 0.00078161 0.00061872 0.00051646 0.00047429 0.000489640 0.00073621 0.0005148 0.00039979 0.00035439 0.000363130 0.00056964 0.00031688 0.00022051 0.0001884 0.000191920 0 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00027006 0.00016305 0.00012104 0.0001071 0.000109950 0.00037754 0.00027817 0.0002253 0.00020483 0.000211160 0.00043539 0.00035386 0.00030207 0.00028124 0.00029140 0.00047398 0.00040104 0.00035113 0.00033111 0.000344050 0.0005003 0.00042535 0.00037309 0.0003519 0.000365710 0.00051479 0.00042699 0.00036681 0.00034164 0.00035410 0.00051415 0.00040056 0.00032887 0.0002985 0.000307640 0.00048504 0.0003335 0.00025411 0.0002221 0.000227060 0.00037609 0.00020532 0.00013956 0.00011718 0.000119020 0 0 0 0 03.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 8.6469e-006 1.4412e-005 1.8364e-005 2.091e-005 2.2174e-0050 1.4412e-005 2.4777e-005 3.2047e-005 3.6716e-005 3.8961e-0050 1.8364e-005 3.2047e-005 4.1789e-005 4.8054e-005 5.1008e-0050 2.091e-005 3.6716e-005 4.8054e-005 5.5353e-005 5.8764e-0050 2.2174e-005 3.8961e-005 5.1008e-005 5.8764e-005 6.2389e-0050 2.2118e-005 3.8698e-005 5.0523e-005 5.8126e-005 6.171e-0050 2.055e-005 3.5581e-005 4.6157e-005 5.2942e-005 5.6197e-0050 1.707e-005 2.8951e-005 3.7128e-005 4.2365e-005 4.4952e-0050 1.0851e-005 1.7698e-005 2.2265e-005 2.5203e-005 2.672e-0050 0 0 0 0 01.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）Columns 7 through 110 0 0 0 00.00020348 0.00026228 0.00038338 0.00066008 00.00038607 0.00048321 0.00064717 0.00091668 00.00052635 0.00064203 0.00081637 0.0010517 00.0006174 0.00074272 0.00092111 0.0011417 00.00065651 0.00078964 0.00097724 0.0012051 00.00064051 0.00078116 0.00098594 0.0012433 00.00056474 0.00070822 0.00093332 0.0012478 00.00042547 0.00055526 0.00078616 0.0011844 00.00022735 0.00030946 0.00049004 0.00092402 00 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)Columns 7 through 110 0 0 0 00.00012826 0.00016798 0.00024986 0.00043637 00.00024444 0.00031023 0.00042173 0.00060513 00.00033401 0.00041257 0.00053179 0.00069358 00.00039216 0.00047742 0.00059986 0.00075263 00.00041704 0.00050761 0.00063642 0.00079439 00.00040657 0.0005021 0.00064226 0.00081984 00.00035784 0.000455 0.00060828 0.00082334 00.00026866 0.00035628 0.00051263 0.0007824 00.00014262 0.00019789 0.00031956 0.0006113 00 0 0 0 0 3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) Columns 7 through 110 0 0 0 02.2118e-005 2.055e-005 1.707e-005 1.0851e-005 03.8698e-005 3.5581e-005 2.8951e-005 1.7698e-005 05.0523e-005 4.6157e-005 3.7128e-005 2.2265e-005 05.8126e-005 5.2942e-005 4.2365e-005 2.5203e-005 06.171e-005 5.6197e-005 4.4952e-005 2.672e-005 06.1116e-005 5.5803e-005 4.4817e-005 2.6785e-005 05.5803e-005 5.1239e-005 4.1529e-005 2.5153e-005 04.4817e-005 4.1529e-005 3.42e-005 2.126e-005 02.6785e-005 2.5153e-005 2.126e-005 1.3869e-005 00 0 0 0 0。
二维抛物型方程的交替方向隐格式在数学领域中，二维抛物型方程是一类重要的偏微分方程，它们在众多实际问题的数学建模中起着关键作用。
对于这类方程，交替方向隐格式是一种常用且有效的数值求解方法。
本文将详细介绍二维抛物型方程及其交替方向隐格式的原理和应用，希望能为读者提供一份生动、全面且有指导意义的参考材料。
首先，我们来了解什么是二维抛物型方程。
通常，二维抛物型方程可以表示为以下形式：∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²) + βu + f(x, y, t)其中，u是未知函数，t是时间变量，x和y是空间变量，α和β是常数，f(x, y, t)是已知函数。
二维抛物型方程广泛应用于物理、工程、生物等领域的问题求解，比如热传导、扩散、扩散反应等。
为了求解这类方程，数学家们开发了各种求解方法，交替方向隐格式是其中一种。
交替方向隐格式是一种时间和空间交替迭代的求解方法，它通过将二维抛物型方程离散化为一组代数方程，然后通过迭代求解这些代数方程得到数值解。
具体来说，交替方向隐格式先将时间方向离散化，将时间变量t划分为一系列离散时间步长。
然后，对于每个时间步长，交替方向隐格式将二维抛物型方程中的时间导数∂u/∂t进行近似。
最常用的近似方法是向后差分格式，即用u(n+1) - u(n)来近似∂u/∂t，其中u(n)表示第n个时间步长的数值解，u(n+1)表示第n+1个时间步长的数值解。
这样，二维抛物型方程可以离散为一组代数方程。
接下来，交替方向隐格式将空间方向离散化，将空间变量x和y划分为一系列离散网格点。
然后，在空间离散化的基础上，通过引入交替方向（例如，先按x方向更新，再按y方向更新，或者反之）和隐格式（例如，使用向后差分格式近似二阶导数项），将二维抛物型方程中的空间导数进行近似。
通过交替迭代求解这组离散代数方程，我们可以得到二维抛物型方程的数值解。
当离散网格点的数量足够多时，数值解将趋近于方程的解。
二维抛物线方程数值解法方法二维抛物线方程是数学中常见的偏微分方程之一，表示一个二维区域中的传热、扩散、波动等现象。
求解二维抛物线方程通常需要使用数值解法，其中ADI（Alternating Direction Implicit）隐式交替法是一种有效的数值解法之一、本文将详细介绍ADI隐式交替法的原理和步骤。
首先，我们来看一下二维抛物线方程的一般形式：∂u/∂t=α(∂²u/∂x²+∂²u/∂y²)+f(x,y,t)其中u(x,y,t)是未知函数，表示二维区域中的一些物理量（如温度），α是扩散系数，f(x,y,t)是源（即外力项）。
ADI隐式交替法的基本思想是将偏微分方程中的时间导数项进行分离，然后通过交替方向的半隐式差分逼近来求解。
具体步骤如下：1.将时间导数项进行分离，得到两个互相独立的方程：∂u/∂t=α(∂²u/∂x²)+α(∂²u/∂y²)+f(x,y,t)等价于：∂u/∂t=α(∂²u/∂x²)+f(x,y,t)，考虑到y方向上的导数已经被分离出来。
∂u/∂t=α(∂²u/∂y²)+f(x,y,t)，考虑到x方向上的导数已经被分离出来。
2. 对第一个方程应用隐式差分（比如Crank-Nicolson差分格式），得到一个关于未知函数u的线性方程组Ax=b。
3.利用ADI思想，对第二个方程同样进行隐式差分，得到另一个关于未知函数u的线性方程组Ay=c。
4.对线性方程组Ax=b和Ay=c进行求解，得到u的一个时间层的近似解。
5.重复步骤2至4，直到得到整个时间区间内的数值解。
ADI隐式交替法的关键在于递归地求解两个互相独立的线性方程组。
在每个时间步长中，先求解一个方向（比如x方向）上的方程组，然后用该方程组的解代入另一个方向（比如y方向）的方程组中求解。
通过交替进行，可以逼近二维抛物线方程的解。
解二阶抛物型方程含参数高精度两层差分格式引言在数值计算中，抛物型偏微分方程是一类重要的方程，涉及到热传导、扩散、反应扩散等许多实际问题的数值模拟。
而解二阶抛物型方程含参数的问题则更具有挑战性，建立高精度的差分格式成为研究的重点之一。
本文将详细探讨解二阶抛物型方程含参数的高精度两层差分格式。
二阶抛物型方程含参数先考虑一个边界平滑、参数依赖的二阶抛物型方程：其中，u(x,t)是未知函数，a(x)是参数函数，f(x,t)是外力项。
本文的目标是构建一个高精度的差分格式来求解这个方程。
高精度差分格式的结构由于我们要构建高精度的差分格式，自然而然地思路是采用两层差分格式的结构。
两层差分格式具有更高的精度和更好的稳定性，对于求解复杂的偏微分方程尤为有效。
一维差分格式首先，我们先介绍一维的差分格式。
对于一维问题，我们可以将区域进行离散化，将连续的问题转换为离散的问题。
一维差分格式可以表示为：这个差分格式包含三个部分，分别是时间项、空间项和外力项。
其中，时间项计算了时间导数，空间项使用中心差分方法计算了空间导数，外力项贡献了外力。
多维差分格式对于二维问题，我们可以将其推广为多维差分格式。
多维差分格式的基本思想是利用乘积形式构造更高维的差分格式。
在这里，我们只考虑二维情况：这个差分格式在时间、x方向和y方向都有相应的差分项。
同样，时间项计算了时间导数，空间项使用了中心差分方法计算了空间导数，外力项贡献了外力。
参数依赖的高精度两层差分格式在介绍了一维和多维差分格式后，我们接下来考虑含参数的高精度两层差分格式。
对于含参数的问题，我们需要在空间项中考虑该参数的影响。
高精度的两层差分格式可以表示为：其中，αi,j n是参数依赖的系数。
高精度差分格式的稳定性和收敛性分析在构建高精度差分格式时，我们不仅需要关注其精度，还需要考虑其稳定性和收敛性。
在这一部分，我们将对所构建的高精度差分格式进行稳定性和收敛性的分析。
稳定性稳定性是差分格式是否收敛的重要条件之一。
解二阶抛物型方程含参数高精度两层差分格式
解二阶抛物型偏微分方程是许多领域中非常重要的问题，例如热传导、扩散、化学反应等。
然而，在实际计算中，由于参数的存在，往往需要使用高精度的差分格式来保证数值计算的准确性。
以下是一种含参数的高精度两层差分格式，可以用于解决二阶抛物型偏微分方程：
1. 先验估计
为了使用高精度两层差分格式，我们首先需要进行先验估计，以确定合适的时间步长和空间步长。
具体来说，我们可以使用稳定性分析来确定时间步长和空间步长的上限值。
2. 差分格式
在确定了时间步长和空间步长之后，我们可以开始使用高精度两层差分格式来求解二阶抛物型偏微分方程。
该差分格式通常包括以下几个步骤：
（1）先用向前差分公式求解第一层，得到一个中间解。
（2）再采用Crank-Nicolson格式对第二层进行差分，同时使用前一步得到的中间解进行修正。
（3）最后，将得到的数值解反推回到未知函数的值域中，得到方程的
数值解。
需要注意的是，在使用这种高精度差分格式进行计算时，我们需要使
用高精度的算法来保证计算的准确性。
3. 参数调节
由于实际问题中经常存在参数不确定性的情况，因此，在进行数值计
算时，我们需要对参数进行调节和优化。
具体来说，我们可以通过多
次求解不同的二阶抛物型偏微分方程，来不断调节参数并逐步优化计
算结果。
以上是一个含参数的高精度两层差分格式，可以用于解决二阶抛物型
偏微分方程的计算问题。
该方法能够保证数值计算的高精度和准确性，同时也能够应对实际问题中的参数不确定性。
李荣华二维抛物方程的初边值问题在数学和工程领域中扮演着重要角色。
通过使用Matlab，我们能够更深入地理解和解决这一问题。
在本文中，我们将从基础概念开始，逐步深入讨论李荣华二维抛物方程的初边值问题，并结合Matlab进行实际分析和解决。
一、初边值问题的概念和涵义1. 初边值问题的定义初边值问题是指在偏微分方程中，除了部分边界条件外，还需给出一些初始条件。
李荣华二维抛物方程的初边值问题即是在方程中给定初值条件和边界条件的问题。
2. 李荣华二维抛物方程简介李荣华二维抛物方程是描述热传导、扩散等现象的数学模型，在物理学和工程领域有着广泛的应用。
它的形式通常为一个关于未知函数和其在空间和时间上的导数的方程。
3. Matlab在初边值问题中的应用Matlab作为一种强大的数学建模和仿真工具，能够帮助我们求解各种类型的偏微分方程，包括初边值问题。
通过Matlab，我们可以对李荣华二维抛物方程进行数值求解和模拟，从而更清晰地理解问题本质。
二、李荣华二维抛物方程的初边值问题解析1. 方程的建立我们需要建立李荣华二维抛物方程的数学模型，并给定相应的初值条件和边界条件。
在Matlab中，我们可以通过定义矩阵和方程的形式来实现这一过程。
2. 数值求解我们利用Matlab提供的数值求解方法，如有限差分法或有限元法，对初边值问题进行求解。
通过程序的编写和运行，我们可以得到方程的数值解，并对物理过程有更直观的认识。
3. 结果分析在求得数值解之后，我们可以对结果进行分析和可视化。
通过Matlab 的绘图功能，我们能够直观地观察李荣华二维抛物方程的演化过程，以及初值条件和边界条件对解的影响。
三、个人观点和总结在本文中，我们深入探讨了李荣华二维抛物方程的初边值问题，并结合Matlab进行了实际分析和解决。
通过对问题的全面评估和数值求解，我们对方程的特性和行为有了更深入的理解。
个人观点上，我认为Matlab是一个非常强大的工具，能够帮助我们更直观、高效地处理数学问题。
基于matlab解抛物型方程的交替隐方向p-r差分格式的实现1. 引言1.1 概述本文旨在利用MATLAB中的抛物型方程解析方法，具体实现交替隐方向p-r差分格式。
抛物型方程是一类常见的偏微分方程，在科学计算和工程领域中有着广泛的应用。
该类方程描述了许多自然界和社会系统中的动态过程，如热传导、扩散、弹性形变等。
而交替隐方向p-r差分格式则是一种高效解法，适用于求解抛物型方程。
1.2 文章结构本文将按以下结构展开详细论述：- 第2节将简要介绍抛物型方程及其解析方法概述，并特别关注MATLAB在此过程中的应用。
- 第3节将深入探讨交替隐方向p-r差分格式的原理，并对其稳定性和精确度进行分析。
- 第4节将重点阐述基于MATLAB实现交替隐方向p-r差分格式的步骤，包括空间离散化方法选择与实现、时间离散化方法选择与实现、以及迭代求解过程描述与收敛性分析。
- 最后，第5节将呈现数值实验设置，并展示数值结果，同时对结果进行讨论。
1.3 目的本文的目的在于通过MATLAB解析抛物型方程，并实现交替隐方向p-r差分格式，从而提供一种高效、稳定、精确的数值计算方法。
此研究对于处理抛物型方程相关问题具有实际应用意义，为科学计算和工程领域中的相关研究提供了指导和借鉴。
我们期望该研究能够拓展数值计算方面的知识，促进在实践中解决复杂系统动态过程模拟与分析的能力。
2. 抛物型方程解析2.1 抛物型方程简介抛物型方程是一类常见的偏微分方程，它描述了许多自然现象和数学模型中的动态行为。
一般而言，抛物型方程包括一个时间变量和多个空间变量，并且通常具有二阶时间导数和二阶或更高阶的空间导数。
典型的抛物型方程包括热传导方程、扩散方程和波动方程等。
2.2 解析方法概述解析方法是指通过使用数学分析和解析推导来求解偏微分方程的方法。
在抛物型方程的解析研究中，常用的方法包括变量分离法、相似变量法、格林函数法等。
这些方法基于物理建模和数学推导，可以得到精确的解或者近似解。
一类二维抛物型方程的ADI格式【摘要】本文针对一类二维抛物型方程，建立了一个在空间和时间方向上均具有二阶精度的ADI格式，并分析其稳定性. 比较以往算法，此格式具有精度较高，无条件稳定等优点，同时，该方法通过求解两个线性代数方程得到原问题的解，避免了非线性迭代运算，提高了计算效率.【关键词】二维抛物型方程；ADI格式；稳定性；截断误差1.引言抛物型偏微分方程在研究热传导过程、一些扩散现象及电磁场传播等许多问题中都有广泛的应用，对这一类方程数值解法的研究一直是科研工作者关注的热点问题之一，其中高精度的有限差分方法更是受到了越来越多的重视. 考虑如下的初边值问题[1]：其中，为常数.在文献[1]中对问题（1）建立了差分格式，格式的截断误差阶为.本文将在文献[1]的基础上进一步研究问题（1）的高效差分格式，建立了一个高精度的交替方向隐式差分格式（即ADI格式），提高了时间方向上的精度，并给出相应的稳定性分析。
2.差分格式的建立为了得到问题（1）的有限差分格式，首先将求解区域进行网格剖分，结点为. 选取正整数L和N，并令为方向上的网格步长，为方向上的网格步长，记假定第层的已知，则由第（Ⅰ）步，通过解一个三对角线性代数方程组求出，再由第（Ⅱ）步，再解一个三对角线性代数方程组即可求出. 所以，只要利用追赶法求解两个三对角线性代数方程组即可，此时计算量与储存量都较小. 另外，在处理边界条件时，为了提高精度，采取中心差分，这样会出现虚拟值，此时，只要再与格式中的方程联立，即可消去虚拟值[2].3. 稳定性分析下面采用von Newmann方法[3]对上述D格式进行稳定性分析. 一般地，低阶项不影响差分格式的稳定性，所以不妨略去项，并对（3）、（5）式消去中间变量得：利用Taylor展开式求误差，可知此处建立的D格式的截断误差阶为.参考文献：[1]管秋琴.一类二维抛物型方程的有限差分格式[J]. 赤峰学院学报（自然科学版）. 2010，26（1）：7.[2]Richtmyer R D ，Morton KW. Difference methods for initial - value problems （2nd edition）[J ] . John wiley &sons. 1967.[3]戴嘉尊，邱建贤. 微分方程数值解法[M]. 南京：东南大学出版社.2002.；安庆师范学院青年科研基金项目（KJ201020）。
二维变系数抛物型方程的一个高阶ADI差分格式马小霞;颜晓琳;陈汝栋【摘要】A high accuracy alternation direction implicit scheme (ADI) for solving the two-dimensional parabolic equations is presented, and the scheme is absolutely stable and the truncation error is O(τ2+h4). The experiments show the scheme is effective and advantage, and the theory is right by a numerical example.%针对二维变系数抛物型方程，构造出了一个高精度、恒稳定的交替方向隐式（ADI）差分格式，格式的截断误差阶达O（τ2+h4）。
通过数值实验，验证了理论分析的正确性和差分格式的精确性与有效性。
【期刊名称】《天津工业大学学报》【年(卷),期】2014(000)001【总页数】5页(P77-80,84)【关键词】抛物型方程;ADI格式;截断误差;恒稳定【作者】马小霞;颜晓琳;陈汝栋【作者单位】焦作大学基础部，河南焦作 454003;天津工业大学理学院，天津300387;天津工业大学理学院，天津 300387【正文语种】中文【中图分类】O241.82抛物型方程在处理废料污染、渗透、驱动、海水入侵以及半导体等工程实际问题中有着广泛的应用，因此研究其高精度、高稳定和计算量较小的数值解法具有重要的意义.用有限差分方法研究这类问题的数值方法目前已做了许多工作[1-5].但这些工作大多是对常系数而言的.文献[4]中对二维变系数抛物型方程数值方法仅对系数依赖于一个变量的情况进行了研究，本文的研究是对系数依赖于两个变量的情形进行的.应用Taylor展开、算子方法[6]以及粘结系数法[7]得到了一个高精度（截断误差阶达O（τ2+h4））、恒稳定的ADI格式.格式的建立和稳定性分析都比文献[4]简单得多，文末的数值实验证明了本文理论分析的正确性和所得格式的精确性与有效性.考虑如下的二维变系数非齐次抛物型方程初边值问题其中：0＜c1≤a（x，y）≤c2；Γ为Ω的边界.设τ=Δt=T/N为时间步长，h=Δx=Δy=1/M为空间步长，N、M均为正整数.表示在节点（jh， kh，nτ）处的网函数值，微分方程问题（1）—（3）的解函数为u（x，y，t），并记u（jh，kh，nτ）=u（j，k，n），由Taylor展开式在（j，k，n+1）处考虑方程（1），根据粘结系数法[7]有其中：分别为关于x、y的二阶中心差分算子；η1、η2、η3均为待定参数.适当选取这些参数使得式（9）—式（11）成立.由于即式（9）—式（11）的3个式子均成立，由式（8）、（14）、（16）和（17）得：将所得的等式左端加上右端加上，经过简单的计算可知，所加项的误差阶均不超过O（τ2+h4），则可得在（18）中略去小量项，并令，则可建立如下差分格式：易知差分格式（19）—（21）的截断误差为对格式（19）—（21）引进中间过渡量，可以得到问题（1）—（3）的ADI格式：当第n层上的值已知时，对于固定的k，（22）是关于的三对角线性方程组，可用追赶法求解.在解方程组（22）时，需要边值和（0≤k≤M-1），这可从方程（23）得到，在方程（23）中，分别令j=0和j=M，有当过渡层已求得时，对于固定的j，式（23）是关于的三对角线性方程组，也可用追赶法求解.解方程组（23）所需要的边界值和（0≤k≤M-1）可直接从已知的边界值条件得到，即：为对格式（22）—（23）进行稳定性分析，可将过渡层变量消去，即得格式（19），故两格式有相同的截断误差和稳定性质.因为差分格式关于初值稳定一定可以推出关于右端稳定，故只对齐次格式进行讨论，根据稳定性分析的Fourier方法[21]，令将上式代入格式（19）式的齐次形式中，经计算整理，并利用关系式这里，可得格式（19）的传播因子为由s1的取值范围可知，所以|G（s1，s2）|根据Von neumann条件知格式（19）也即格式（22）—（23）恒稳定.综合上节与本节的论述，并根据Lax的稳定性与收敛性等价原理可得：定理本文的ADI格式恒稳定且以o（τ2+h4）的收敛阶收敛.在区域D：{0≤x，y≤1，t≥0}上对初值边问题用本文格式求数值解并与精确解u（x，y，t）=e-2tsin（x+y），相比较，取h=1/10，τ=rh2=r/100，r=1/2计算到n=200时的结果如表1所示.由表1结果看出，对所取的r，本文格式解与精确解均有很好的吻合，这与理论分析完全一致.本文构造了二维变系数抛物型方程的ADI格式，数值例子表明，它具有易于计算（能三对角追赶法求解）、精度高（较现有格式）、无条件稳定等特点.本文构造出的格式可以推广到任意维变系数抛物型方程，对当前关于变系数抛物型方程差分格式的构造与研究有重要的理论和实践意义.【相关文献】[1]DAI Weizhong，NASSAR Raja.A Second-order ADI scheme for three-dimensional parabolic differential equations[J].Numer Math，1964（6）：428-453.[2]DAI Weizhong，NASSAR pact ADI method for solving parabolic differential equations[J].Numer Methods Partial Differential Eq，2002（18）：129-142.[3] 孙志忠，李雪玲.反应扩散方程的紧交替方向差分格式[J].计算数学，2005，19（2）：209-224.[4] 孙志忠，李雪玲.二维变系数反应扩散方程的紧交替方向差分格式[J].高等学校计算数学学报，2006，28（1）：83-95.[5] 马菊意，杨辉.二维抛物型方程的一个高精度PC格式[J].工程数学学报，2008，25（2）：373-376.[6]余德浩，汤华中.微分方程数值解法[M].北京：科学出版社，2003：102-174.[7] 戴嘉尊，邱建贤.微分方程数值解法[M].南京：东南大学出版社，2003：91-93[8]李立康，於崇华，朱政华.微分方程数值方法[M].上海：复旦大学出版社，1999：208-213.。
ADI 法求解二维抛物方程学校：中国石油大学（华东） 学院：理学院 姓名：张道德 时间：2013.4.271、ADI 法介绍作为模型，考虑二维热传导方程的边值问题：（3.6.1）,0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0t xx yy u u u x y l t u x y x y u y t u l y t u x t u x l t ϕ=+<<>⎧⎪=⎨⎪====⎩取空间步长1hM，时间步长0，作两族平行于坐标轴的网线：,,,0,1,,,j k x x jh y y kh j k M =====将区域0,x y l ≤≤分割成2M 个小矩形。
第一个ADI 算法（交替方向隐格式）是Peaceman 和Rachford （1955）提出的。
方法：由第n 层到第n+1层计算分为两步：（1） 第一步： 12,12n j k xx yy u +从n->n+，求u 对向后差分,u 向前差分，构造出差分格式为：1(3.6.1)11112222,,1,,1,,1,,1221222,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ+++++-+-+-+-+=+uu uuuu u u （+）u u（2） 第二步：12,12n j k xx yy u +从n+->n+1，求u 对向前差分,u 向后差分，构造出差分格式为：2(3.6.1)1111111222,,1,,1,,1,,12212212,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hh hτδδ++++++++-+-++-+-+=+uu uuuu u u （+）u u其中1211,1,,1,0,1,2,,()22n j k M n n n τ+=-=+=+上表表示在t=t 取值。
如何用matlab求解二阶微分方程，以及程序实例(2010-03-25 12:56:36)转载matlab程序实例二阶微分方程dsolve文化微分方程的求解主要用到一个dsolve函数，如下面的“ Uc仁dsolve('D2y+1000*Dy+10A6*y=0','y(0)=10','Dy(0)=0','t'); ”，可以看出，函数的第一部分是所要求解的微分方程，其次是初始条件，最后是对自变量的说明。
下面给出的程序实例是用于分析一个最简单零输入的二阶电路。
其中C=1uf,L=1H。
R是不确定的，他的值的选取将会直接影响到方程解的形式以及最后画出的曲线形状，在此我取R的值分别为1000,2000,3000欧姆。
R=1000;while (R<=3000)if R<2000 Uc1=dsolve('D2y+1000*Dy+10A6*y=0','y(0)=10','Dy(0)=0','t');It1=-1*diff(Uc1)*(1e-6);Ul1=diff(It1);elseif R==2000 Uc2=dsolve('D2y+2000*Dy+10A6*y=0','y(0)=10','Dy(0)=0','t');It2=-1*diff(Uc2)*(1e-6);Ul2=diff(It2);elseUc3=dsolve('D2y+3000*Dy+10A6*y=0','y(0)=10','Dy(0)=0','t');lt3=-1*diff(Uc3)*(1e-6);Ul3=diff(It3);endR=R+1000;end%while i<=3figure (1)xlabel('t') ylabel('Uc(t)') hold on;。
ADI 法求解二维抛物方程学校：中国石油大学（华东）学院：理学院：道德时间：2013.4.271、ADI 法介绍作为模型，考虑二维热传导方程的边值问题：（3.6.1）,0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0t xx yy u u u x y l t u x y x y u y t u l y t u x t u x l t ϕ=+<<>⎧⎪=⎨⎪====⎩取空间步长1hM，时间步长0，作两族平行于坐标轴的网线：,,,0,1,,,j k x x jh y y kh j k M =====将区域0,x y l ≤≤分割成2M 个小矩形。
第一个ADI 算法（交替方向隐格式）是Peaceman 和Rachford （1955）提出的。
方法：由第n 层到第n+1层计算分为两步：（1） 第一步： 12,12n j k xx yy u +从n->n+，求u 对向后差分,u 向前差分，构造出差分格式为：1(3.6.1)11112222,,1,,1,,1,,1221222,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ+++++-+-+-+-+=+uu uuuu u u （+）u u（2） 第二步：12,12n j k xx yy u +从n+->n+1，求u 对向前差分,u 向后差分，构造出差分格式为：2(3.6.1)1111111222,,1,,1,,1,,12212212,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hh hτδδ++++++++-+-++-+-+=+uu uuuu u u （+）u u其中1211,1,,1,0,1,2,,()22n j k M n n n τ+=-=+=+上表表示在t=t 取值。
假定第n 层的,n j ku已求得，则由1(3.6.1)求出12,n j k u +，这只需按行(1,,1)j M =-解一些具有三对角系数矩阵的方程组；再由2(3.6.1)求出1,n j k u +，这只需按列(1,,1)k M =-解一些具有三对角系数矩阵的方程组，所以计算时容易实现的。
2、数值例子（1）问题用ADI 法求解二维抛物方程的初边值问题：21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sin cos .xx yy y y u u u x y G t t u y t u y t y t u x t u x t x t u x y x y ππ∂⎧=+∈=⨯>⎪∂⎪⎪==<<>⎨⎪==<<>⎪=⎪⎩ 已知（精确解为：2(,,)sin cos exp()8u x y t x y t πππ=-）设(0,1,,),(0,1,,),(0,1,,)j k n x jh j J y kh k K t n n N τ======差分解为,n j k u ，则边值条件为：0,,,0,1,1,0,0,1,,,,0,1,,n nk J k nn n nj j j K j K u u k Ku u u u j J-⎧===⎪⎨===⎪⎩初值条件为：0,sin cos j k j k u x y ππ=取空间步长12140h h h ===，时间步长11600τ=网比21r h τ==。
用ADI 法分别计算到时间层1t =。
（2）计算过程从n 到n+1时，根据边值条件：0,,0,0,1,,n n k J k u u k K ===，已经知道第0列和第K 列数值全为0。
（1）12,12n j k xx yy u +从n->n+，求u 对向后差分,u 向前差分，构造出差分格式为：11112222,,1,,1,,1,,1221222,,2-221=1621()16n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ+++++-+-+-+-+=+u u u u uu u u （+）u u从而得到：1112221,,1,,1,,1111111(1)(1)321632321632n n n n n n j k j k j k j k j k j k r r r r r r ++++-+--++-=+--u u u u u u ，其中1,2,,1,1,2,,1j J k K =-=-即按行用追赶法求解一系列下面的三对角方程组：121,122,123,123,122,121,(1)(1)111163211113216321111321632111132163211113216321113216n k n k n k n J k n J k n J k J J r r r r r r rr r r r r r r r r ++++-+-+--⨯-⎡⎤⎡+-⎢⎥⎢⎢⎥⎢⎢⎥-+-⎢⎢⎥⎢⎢⎥⎢⎢⎥-+-⎢⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+⎣⎢⎥⎣⎦u u u u u u 123321(1)1(1)1J J J J J f f f f f f ----⨯-⨯⎤⎥⎥⎡⎤⎥⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥⎦ 又根据边值条件得：,0,1,1,,,0,1,,nnnnj j j K j K u u u u j J -===，解出第0行,0nj u 和第K 行,,(0,1,,)n j K u j J =。
（2）第二步：12,12n j k xx yy u +从n+->n+1，求u 对向前差分,u 向后差分，构造出差分格式为：1111111222,,1,,1,,1,,12212212,,2-221=1621()16n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ++++++++-+-++-+-+=+uu u uuu u u （+）u u从而得到：111111222,1,,11,,1,111111(1)(1)321632321632n n n n n n j k j k j k j k j k j k r r r r r r +++++++-+--++-=+--u u u u u u ，其中1,2,,1,1,2,,1j J k K =-=-又根据边值条件得：,0,1,1,,,0,1,,nnnnj j j K j K u u u u j J -===，从而得到：,0,1,1,0n nj j n nj K j K u u u u -⎧-=⎪⎨-+=⎪⎩其中(0,1,,)j J = 即按列用追赶法求解一系列下面的三对角方程组：1,01,11,21,31,31,21111113216321111321632111132163211113216321111321632111132163211n j n j n j n j n j K n j K K Kr r r r r r r r r r r rr r r r r r +++++-+-⨯-⎡⎤⎢⎥⎢⎥-+-⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥-⎢⎥⎣⎦u u u u u u u 12343211,111,1K K n K j K K K n j K K f f f f f f f f --+--⨯+⨯⎡⎤⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦u （3） 求解结果从而得到误差的数为：1- 数：0.3713； 2-数：0.1447； ∞-数：0.6086（3.4）图像（3.4.1）数值解图像：（3.4.2）精确解图像：（5）主要程序（5.1）主程序%*************************************************************%main_chapter主函数%信息10-2——道德%学号：10071223clccleara = 0; b=1; %x取值围c=0; d=1; %y取值围tfinal = 1; %最终时刻t=1/1600;%时间步长；h=1/40;%空间步长r=t/h^2;%网比x=a:h:b;y=c:h:d;%**************************************************************%精确解m=40;u1=zeros(m+1,m+1);for i=1:m+1,for j=1:m+1u1(j,i) = uexact(x(i),y(j),1);endend%数值解u=ADI(a,b,c,d,t,h,tfinal);%***************************************************************** %绘制图像figure(1) ;mesh(x,y,u1)figure(2); mesh(x,y,u)%误差分析error=u-u1;norm1=norm(error,1);norm2=norm(error,2);norm00=norm(error,inf);%*****************************************************************（5.2）ADI函数%****************************************************************% 用ADI法求解二维抛物方程的初边值问题% u_t = 1/16（u_{xx} + u_{yy}）（0,1）*（0,1）% 精确解: u(t,x,y) = sin(pi*x) sin(pi*y)exp(-pi*pi*t/8) %****************************************************************function [u]=ADI(a,b,c,d,t,h,tfinal )%(a , b) x取值围%(c, d) y取值围%tfinal最终时刻%t时间步长；%h空间步长r=t/h^2;%网比m=(b-a)/h;%n=tfinal/t; %x=a:h:b;y=c:h:d;%****************************************************************** %初始条件u=zeros(m+1,m+1);for i=1:m+1,for j=1:m+1u(j,i) = uexact(x(i),y(j),0);endend%****************************************************************** u2=zeros(m+1,m+1);a=-1/32*r*ones(1,m-2);b=(1+r/16)*ones(1,m-1);aa=-1/32*r*ones(1,m);cc=aa;aa(m)=-1;cc(1)=-1;bb=(1+r/16)*ones(1,m+1);bb(1)=1;bb(m+1)=1;for i=1:n%**************************************************************%从n->n+1/2,u_{xx}向后差分，u_{yy}向前差分for j=2:mfor k=2:md(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1))+u(j,k);end% 修正第一项与最后一项,但由于第一项与最后一项均为零，可以省略%d(1)=d(1)+u1(j,1);d(m-1)=d(m-1)+u1(j,m+1);u2(j,2:m)=zhuiganfa(a,b,a,d);endu2(1,:)=u2(2,:);u2(m+1,:)=u2(m,:);%**************************************************************%从n->n+1,u_{xx}向前差分，u_{yy}向后差分for k=2:mdd(1)=0;dd(m+1)=0;for j=2:mdd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k))+u2(j,k); endu(:,k)=zhuiganfa(aa,bb,cc,dd);end%****************************************************************u2=u;end%********************************************************************（5.3）“追赶法”程序%******************************************************************** %追赶法function [x]=zhuiganfa(a,b,c,d)%对角线下方的元素，个数比A少一个% %对角线元素%对角线上方的元素，个数比A少一个%d为方程常数项%用追赶法解三对角矩阵方程r=size(a);m=r(2);r=size(b);n=r(2);if size(a)~=size(c)|m~=n-1|size(b)~=size(d)error('变量不匹配,检查变量输入情况!');end%%%LU分解u(1)=b(1);for i=2:nl(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);v(i-1)=(b(i)-u(i))/l(i-1);end%追赶法实现%%%求解Ly=d，"追"的过程y(1)=d(1);for i=2:ny(i)=d(i)-l(i-1)*y(i-1);end%%%求解Ux=y,"赶"的过程x(n)=y(n)/u(n);for i=n-1:-1:1x(i)=y(i)/u(i);x(i)=(y(i)-c(i)*x(i+1))/u(i);end%********************************************************************（5.4）精确解函数%t时刻，u的取值；function [ f]=uexact(x,y,t)f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);%********************************************************************。
