1、ADI 法求解二维抛物方程学校：中国石油大学（华东） 学院：理学院 姓名：张道德 时间：2013.4.271、ADI 法介绍作为模型，考虑二维热传导方程的边值问题：（3.6.1）,0,0(,)(),(,)(,)0txyultltuxtlt取空间步长 ，时间步长 ，作两族平行于坐标轴的网线：1hM=0t将区域 分割成 个小矩形。,1,jkxyj ,xyl2M第一个 ADI算法（交替方向隐格式）是 Peaceman和 Rachford（1955）提出的。方法：由第 n层到第 n+1层计算分为两步：（1） 第一步： ，构造出差分格12,njkxyu从 -+， 求 对 向 后 差 分 ,向 前 差 分
2、式为： 1(3.6)111222,1,12122,-=()nnnnnjkjjkjkjkjjkjnnxjkyjkhhuu（ ） u（2） 第二步： ，构造出差分12,2njkxy从 n+-， 求 对 向 前 差 分 ,u向 后 差 分格式为： 2(3.61)111122,21221,-=()nnn nnjkjjkjkjkjjkjnnxjkyjkhh u u（ ） u其中 。12,1,0, ()njkM 上 表 表 示 在 t=取 值假定第 n层的 已求得，则由 求出 ，这只需按行,njku1(3.6)12,njku解一些具有三对角系数矩阵的方程组；再由 求出(1,)jM 2(3.61)，这只需按
3、列 解一些具有三对角系数矩阵的方程组，,njku(1,)所以计算时容易实现的。2、数值例子（1）问题用 ADI法求解二维抛物方程的初边值问题： 2(),(0,1),0,40,)1(,),sinco.xyyyuuGtttytt xuxx已知（精确解为： ）2(,)sincoep()8uytxyt设 差分解为 ，(0,1,(0,1),0,1)j k nxhJhKtN ,njku则边值条件为： ,0,1,nnnjjKjuuJ 初值条件为： ,sicojkjkxy取空间步长 ，时间步长 网比 。用 ADI法分124h16021rh别计算到时间层 。t（2）计算过程从 n到 +时 ，根据边值条件： ，已
4、经知道第 0列和第 K列数值全为0,0,1nkJuK0。（1） ，构造出差分格式为：12,njkxy从 -+， 求 对 向 后 差 分 ,u向 前 差 分111222,1,12122,-=6()nnnnnjkjjkjkjkjjkjnnxjkyjkhhuuuu（ +） 从而得到：，其中111222, , ,1,1() ()363632nnnnnnjkjkjkjkjkjkrrrrrr uuuu,JK 即按行用追赶法求解一系列下面的三对角方程组： 12,123,123,12,(1)16321632116326nknknJknJkJrrrrr uu123321()()Jff又根据边值条件得： ，解出第
5、 0行 和,0,1,0,nnnjjjKjuu ,nju第 行 。K,()njuJ（2）第二步： ，构造出差分格12,2njkxyu从 +-， 求 对 向 前 差 分 ,向 后 差 分式为： 111 1122,21221,-=62()nnn nnjkjjkjkjkjjkjnnxjkyjkhh u u（ ） u从而得到：，其中111111222, , , , ,() ()3263263nnnnnnjkjkjkjkjkjkrrrrrr u u1,2,12,jJkK 又根据边值条件得： ，,0,1,0,1nnnjjjjKuuJ从而得到：其中,0,1,njjnjKju(,)j即按列用追赶法求解一系列下面
6、的三对角方程组： 1,0,213,12113263216321632163njjnjKnKrrrrrrr uu1234321,j Kfff（3） 求解结果(3.1)数值解y x 1/4 2/4 3/41/4 0.142057658092578 0.200899866713484 0.1420576580925782/4 2.16292994886484e-15 3.03768181457584e-15 2.12330312762773e-153/4 -0.142057658092571 -0.200899866713473 -0.142057658092570（3.2）精确解y x 1/4 2
7、/4 3/41/4 0.145606466607010 0.205918639844859 0.1456064666070102/4 1.26088801585392e-17 1.78316493265431e-17 1.26088801585392e-173/4 -0.145606466607010 -0.205918639844859 -0.145606466607010（3.3）数值解-精确解（即误差）y x 1/4 2/4 3/41/4 -0.00354880851443196 -0.00501877313137564 -0.003548808514432732/4 2.1503210
8、6870631e-15 3.01985016524929e-15 2.11069424746919e-153/4 0.00354880851443973 0.00501877313138652 0.00354880851444026从而得到误差的范数为：1- 范数：0.233770443573713； 2-范数：0.196807761631447；-范数：0.327253314506086（3.4）图像（3.4.1）数值解图像：0 0.20.4 0.60.8 100.51-0.4-0.200.20.4x、 、 、 、 、 、 、 、y、t、（3.4.2） 精确解图像：0 0.20.4 0.60
9、.8 100.51-0.4-0.200.20.4x、 、 、 、 、 、 、 、y、t、（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/h2;%网比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%
10、数值解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=A
11、DI(a,b,c,d,t,h,tfinal )%(a , b) x取值范围%(c, d) y取值范围%tfinal最终时刻%t时间步长；%h空间步长r=t/h2;%网比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
12、(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,:)=u
13、2(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=
14、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);%*
