ADI交替方向隐格式

胡坤1 ，任兰兰2

1．ADI算法的具体描述

ADI算法又称交替方向隐格式，该算法主要考虑二维热传导方程的边值问题，模型如下：

utuxxuyy,0x,yl,t0 u(x,y,0)(x,y)u(0,y,t)u(l,y,t)u(x,0,t)u(x,l,t)0

在上述模型中，取空间步长h=1，时间步长t>0，作两族平行于坐标轴的网线：

Mxxjjh,yykkh,j,k0,1,L,M,将区域0x,yl分割成M2个小矩形。

具体步骤是将第n层到第n+1层计算分为两步：

1n1（1） 第一步： 从n->n+，求uj,k2对uxx向后差分,uyy向前差分，构造出差分格式为：

2un12j,k-u2nj,k=（un12j1,k2unh12j,k2un12j1,k+unj,k12unj,kunj,k1h2）

12n1 =2(xuj,k2y2unj,k)h1n1（2） 第二步：从n+->n+1，求uj,k2对uxx向前差分,uyy向后差分，构造出差分格式

2为：

un1j,k-u2n+12j,k=（un12j1,k2unh12j,k2un12j1,k+unj,k112unj,k1unj,k11h2）

12n1 =2(xuj,k2y2unj,k1)h其中j,k1,L,M1,n0,1,2,L,其中n11表示在t=t1(n)取值

n222n12假定第n层的unj,k已求得，则由上述第一步可求出uj,k，这只需按行

1(j1,L,M1)解一些具有三对角系数矩阵的方程组；再由第二步求出unj,k，这

只需按列(k1,L,M1)解一些具有三对角系数矩阵的方程组。

2．以ADI算法分析具体实例

（1）考察例子

u1t42(uxxuyy),(x,y)G(0,1)(0,1),t0,u(0,y,t)u(1,y,t)0,0y1,t0, uy(x,0,t)uy(x,1,t)0,0x1,t0u(x,y,0)sinxcosy.上述方程精确解为：u(x,y,t)exp(（2）分析计算过程

首先设xjjh(j0,1,L,J),ykkh(k0,1,L,K),tnn(n0,1,L,N)差分

nnuu0,kJ,k0,k0,1,L,Kn解为uj,k，则边值条件为：n nnnuj,0uj,1,uj,K1uj,K,j0,1,L,J28t)sinxcosy）

初值条件为：u0j,ksinxjcosyk.

取空间步长hh1h2140，时间步长11600网比rh21，用ADI法分别计算到时间层t1.

nn从n到n+1时，根据边值条件：u0,,L,K，已经知道第0列和kuJ,k0,k0,1第K列数值全为0，故：

第一步：

1n1从n->n+，求uj,k2对uxx向后差分,uyy向前差分，构造出差分格式为：

2n12j,kn12j1,kn12j,k2n12j1,ku-u2nj,k1u=（162uuh+unj,k12unj,kunj,k1h2）

1n122n2 =(uxj,kyuj,k)216h从而得到：

11nn1n1111n11nn222ruj1,(1r)ururu(1r)uruj,k1，其中kj,kj1,kj,k1j,k321632321632j1,2,L,J1,k1,2,L,K1

即按行用追赶法求解一系列下面的三对角方程组：

111rr16321r11r1r32163211r1r3216Onn12u1,kn1f1u2,k2f12n21f3ru3,k32OO111n1r1rrfJ32321632uJ3,kf1111J2r1rrun2fJ1(J1)1321632J2,k1n211r1r(J1)1uJ1,k(J1)(J1)3216

又根据边值条件得：uj,0uj,1,uj,K1uj,K,j0,1,L,J，解出第0行unj,0和第

nnn,L,J). K行unj,K,(j0,1第二步：

1n1从n+->n+1，求uj,k2对uxx向前差分,uyy向后差分构造出差分格式为：

2un1j,k-u2n+12j,k1u=（16n12j1,k2unh12j,k2un12j1,k+unj,k112unj,k1unj,k11h2）

1n122n12 =(uuj,k)xj,ky216h从而得到：

1n1n111n11n111n1n1222ruj,k1(1r)uj,kruj,k1ruj1,k(1r)uj,kruj1,k 321632321632其中j1,2,L,J1,k1,2,L,K1

又根据边值条件得：uj,0uj,1,uj,K1uj,K,j0,1,L,J，从而得到：

nnuj,0uj,10其中(j0,1,L,J) nnuj,K1uj,K0再按列用追赶法求解一系列下面的三对角方程组：

nnnn111r11r1r321632111r1rr321632111r1rr321632OOO111r1rr321632111r1rr321632111r1rr32163211KK1unj,0n1uj,1n1uj,2n1uj,3M1unj,K3n1uj,K2n1uj,K1n1uj,Kf1f2f3 f4MfK3fK2fK1fKK1K1从而得到新的时间层的数值解.

3.MATLAB编程实现上述实例

clear clc

a = 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+1

u1(j,i) = uexact(x(i),y(j),1); end end

%数值解

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);

%---------------------------------------------------------------- 编写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+1

u(j,i) = uexact(x(i),y(j),0); end end

%------------------------------------------------------------------- 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:m

for k=2:m

d(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); end

u2(1,:)=u2(2,:); u2(m+1,:)=u2(m,:);

%--------------------------------------------------------------- %从n->n+1,u_{xx}向前差分，u_{yy}向后差分 for k=2:m

dd(1)=0;dd(m+1)=0; for j=2:m

dd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k))+u2(j,k); end

u(:,k)=zhuiganfa(aa,bb,cc,dd); end

%---------------------------------------------------------------- u2=u; end

%------------------------------------------------------------------ “追赶法”解三对角线性方程函数文件

%-------------------------------------------------------------------- %追赶法

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:n

l(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:n

y(i)=d(i)-l(i-1)*y(i-1); end %%

%求解Ux=y,\"赶\"的过程 x(n)=y(n)/u(n); for i=n-1:-1:1

x(i)=y(i)/u(i);

x(i)=(y(i)-c(i)*x(i+1))/u(i); end

%----------------------------------------------------------------- 精确解函数文件 %t时刻，u的取值；

function [ f]=uexact(x,y,t)

f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);

%-------------------------------------------------------------------

图二、精确解图像0.40.2t轴0-0.2-0.4110.500.200.4x轴0.60.8y轴 图1 以方程精确解所绘出的网格图

图一、数值解图像0.40.2t轴0-0.2-0.4110.500.200.4x轴0.60.8y轴图2 以ADI算法得出方程数值解绘出的网格图