ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)

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);%********************************************************************。