ADI(交替方向隐格式)

1ADI算法的MATLAB编程应用实例胡坤1，任兰兰21．ADI算法的具体描述ADI算法又称交替方向隐格式，该算法主要考虑二维热传导方程的边值问题，模型如下：,0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0txxyyuuuxyltuxyxyuytulytuxtuxlt在上述模型中，取空间步长1hM=，时间步长0t，作两族平行于坐标轴的网线：,,,0,1,,,jkxxjhyykhjkM将区域0,xyl分割成2M个小矩形。具体步骤是将第n层到第n+1层计算分为两步：（1）第一步：12,12njkxxyyu从n-n+，求u对向后差分,u向前差分，构造出差分格式为：11112222,,1,,1,,1,,1221222,,2-22=21=()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）uu（2）第二步：12,12njkxxyyu从n+-n+1，求u对向前差分,u向后差分，构造出差分格式为：1111+11112222,,1,,1,,1,,12212212,,2-22=21=()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）uu其中1211,1,,1,0,1,2,,()22njkMnnn其中表示在t=t取值假定第n层的,njku已求得，则由上述第一步可求出12,njku，这只需按行(1,,1)jM解一些具有三对角系数矩阵的方程组；再由第二步求出1,njku，这只需按列(1,,1)kM解一些具有三对角系数矩阵的方程组。22．以ADI算法分析具体实例（1）考察例子21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sincos.xxyyyyuuuxyGttuytuytytuxtuxtxtuxyxy上述方程精确解为：2(,,)exp()sincos8uxyttxy）（2）分析计算过程首先设(0,1,,),(0,1,,),(0,1,,)jknxjhjJykhkKtnnN差分解为,njku，则边值条件为：0,,,0,1,1,0,0,1,,,,0,1,,nnkJknnnnjjjKjKuukKuuuujJ初值条件为：0,sincosjkjkuxy.取空间步长12140hhh，时间步长11600网比21rh，用ADI法分别计算到时间层1t.从n到n+1时，根据边值条件：0,,0,0,1,,nnkJkuukK，已经知道第0列和第K列数值全为0，故：第一步：12,12njkxxyyu从n-n+，求u对向后差分,u向前差分，构造出差分格式为：11112222,,1,,1,,1,,1221222,,2-221=1621=()16nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）uu从而得到：1112221,,1,,1,,1111111(1)(1)321632321632nnnnnnjkjkjkjkjkjkrrrrrruuuuuu，其中1,2,,1,1,2,,1jJkK3即按行用追赶法求解一系列下面的三对角方程组：121,122,123,123,122,121,(1)(1)111163211113216321111321632111132163211113216321113216nknknknJknJknJkJJrrrrrrrrrrrrrrrruuuuuu123321(1)1(1)1JJJJJffffff又根据边值条件得：,0,1,1,,,0,1,,nnnnjjjKjKuuuujJ，解出第0行,0nju和第K行,,(0,1,,)njKujJ.第二步：12,12njkxxyyu从n+-n+1，求u对向前差分构造出差,u向后差分分格式为：1111+11112222,,1,,1,,1,,12212212,,2-221=1621=()16nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）uu从而得到：111111222,1,,11,,1,111111(1)(1)321632321632nnnnnnjkjkjkjkjkjkrrrrrruuuuuu其中1,2,,1,1,2,,1jJkK又根据边值条件得：,0,1,1,,,0,1,,nnnnjjjKjKuuuujJ，从而得到：,0,1,1,00nnjjnnjKjKuuuu其中(0,1,,)jJ再按列用追赶法求解一系列下面的三对角方程组：41,01,11,21,31,31,21111113216321111321632111132163211113216321111321632111132163211njnjnjnjnjKnjKKKrrrrrrrrrrrrrrrrrruuuuuuu12343211,111,1KKnKjKKKnjKKffffffffu从而得到新的时间层的数值解.3.MATLAB编程实现上述实例clearclca=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);fori=1:m+1,forj=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)%误差分析5error=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);fori=1:m+1,forj=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;fori=1:n%---------------------------------------------------------------%从n-n+1/2,u_{xx}向后差分，u_{yy}向前差分forj=2:mfork=2:md(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1))+u(j,k);end6%修正第一项与最后一项,但由于第一项与最后一项均为零，可以省略%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}向后差分fork=2:mdd(1)=0;dd(m+1)=0;forj=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%------------------------------------------------------------------“追赶法”解三对角线性方程函数文件%--------------------------------------------------------------------%追赶法function[x]=zhuiganfa(a,b,c,d)%对角线下方的元素，个数比A少一个%%对角线元素%对角线上方的元素，个数比A少一个%d为方程常数项%用追赶法解三对角矩阵方程r=size(a);m=r(2);r=size(b);n=r(2);ifsize(a)~=size(c)|m~=n-1|size(b)~=size(d)error('变量不匹配,检查变量输入情况!');end%%%LU分解u(1)=b(1);fori=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%追赶法实现7%%%求解Ly=d，追的过程y(1)=d(1);fori=2:ny(i)=d(i)-l(i-1)*y(i-1);end%%%求解Ux=y,赶的过程x(n)=y(n)/u(n);fori=n-1:-1:1x(i)=y(i)/u(i);x(i)=(y(i)