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

1、ADI法求解二维抛物方程学校：中国石油大学(华东)学院：理学院 姓名：张道德 时间：2013.4.271、ADI法介绍作为模型，考虑二维热传导方程的边值问题：Ut =Uxx Uyy,0 :二 X, y ； l,t 0(3.6.1 ) <u(x,y,0) =?(x,y)u(0, y,t) =u(l, y,t) =u(x,0,t) =u(x,l,t) =0取空间步长h=1M ,时间步长t>0,作两族平行于坐标轴的网线： x=xj =jh,y = yk=kh,j,k=0,1,川，M,将区域 0 M x, y M l 分割成 M 2 个小矩形。第一个ADI算法(交替方向隐格式) 是Peac

2、emarff口 Rachford (1955)提出的。方法：由第n层到第n+1层计算分为两步：. r 一1 n 1(1)第一步： 从n->n+；求Uj,k2对u。向后差分，uyy向前差分,构造出差分格式为:,n 2 , ,nUj,k -Uj,k(361)12nJ=(5n '1 n:：；1-2%二h2可上书-2un,k+un,fh1o n 1-J2G2Ujk2 卜2 x j，k(2)第二步:1从 n+->n+12n 1，求 Uj,k2对 u打向前差分,uyy向后差分，构造出差分格式为:n 1 nUj,k-Uj,k(3.6.1)22n：dn 1n 1-Uj .2,k -2Uj,

3、k2 , Uj2kh2n 1n -1n 1Uj,k 书2Uj,k +uj,kh223力；？+62%：) h 2 x j,k y j,k1 - -j =(n +一2取值 n -O2L1 .其中 j,k=1,lll,M 1,n =0,1,2,111，上表 n+一表示在 t=t2n 1假定第n层的u；,k已求得，则由(3.6.1)1求出Uj,k2 ,这只需按行(j =1,111,M -1)解一些具有三对角系数矩阵的方程组；再由 (3.6.1)2求出u；,这只需按列(k=1,|,M -1)解一些具有三对角系数矩阵的方程组,所以计算时容易实现的。2、数值例子(1)问题用ADI法求解二维抛物方程的初边值问

4、题：四=! (Uxx + Uyy),( x,y) w G = (0,1) M (0,1),t > 0, 涿 4u(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.2已知(精确解为：u(x, y,t) =sin nxcosnyexp(-彳t)设 xj = jh(j =0,1/11 J), yk =kh(k =0,1,|,K),tn =nf(n =0,1,| ,N)差分解为 /水 ,则边值条件为:un,k =un,

5、k =0,k=0,1川,Kuj,0 =Uj,1,Uj,K=Uj,K,j =0,1JII,J初值条件为：u0,k = sin 二 xj cos二 yk取空间步长h =h1 =h2 =1/40,时间步长7=1/1600网比r =dh2 =1。用ADI法分别计算到时间层t=1。(2)计算过程从n到n+1时，根据边值条件：un,k=u：k=0,k=0,1,|i,K ,已经知道第0列和第K列数值全为01n 1_(1)从n->n+ ,求ujk2对uxx向后差分，u yy向前差分,构造出差分格式为：2j ,xxyyn ：；-n ：!1n：!|ln：!1Uj,k2-un,k1 / Uj/,k 2Uj,k

6、2+Uj_2k un,k书2un,k+un,k二、3P=w(h2+h2)1c n 1=G2Ujk2 , Wk)16卜2x j ,ky j,k从而得到：132116 十1/V+132uk1.-24,n ju1；1、n 1 n=3215+。一石厂皿水一而四,其中j =1,2J|,J -1,k =1,2"H，K1即按行用追赶法求解一系列下面的三对角方程组:11 +r161 r 321r3211 r161一一r321一一 r32,11 r161r321r321 r161r321一r321 r161- r321 r 32n +Uj或n+Uj IkngjUJ_Lk 1(jnK一 f1 1f2f3

7、又根据边值条件得:Uj' = U；1,UjnK=un,K , J =0,1,III,J ,解出第。行 Un,0 和第K 行 u；K,(j =0,1/礼J)。(2)第二步:1从 n+ >n+12n 1求Uj,k2对UXX向前差分，Uyy向后差分，构造出差分格式为:unUn -n - n -u 九 - 2u -2，u 5 uni-j 1,k UJ j,kJ j,k + j,k 1h22成；u；,;h21, 2 n 22 n 1、南( xj2' yUj,k)从而得到:1 n 11、n 1-ruj，k1 (1 r)Uj'k 32I61n 11一一ruj k_i = 一 r

8、u32 j，32n 2 j 1,k1(1 r)u16n 21j k - 一 ruj， 32n -jj,k,其中j =1, 2U J,-七 HH , k,-又根据边值条件得：un,0 =un,1,un,K =un,K, j =0,1,III,J ,从而得到:'n nUj,0 - Uj,1 = 0nnUj,KJuj,K其中(j =0,1,111, J)二0即按列用追赶法求解一系列下面的三对角方程组:- 1-11111r 1 +- r r71n 4 7QOUaQOU i c2 J1 V,J£-II J 1/1U：才_r 1 +r_rj，1321632n 芈Uj21-11-r1 +r

9、 一一r.：4321632Uj，3-ii-111u：九-r 1 +r- rj,K-321632.n +Uj,K 21一 11一-r1+r-r.n+321632Uj，K。n-1一 f1 1f2f3f4fK_3fK_2fKfKKJ111U j K-r 1+r 一一r，_K4321632 , 14(3)求解结果(3.1)数值解y J.x1/42/43/41/40.1420576580925780.2008998667134840.1420576580925782/42.16292994886484e-153.03768181457584e-152.12330312762773e-153/4-0.14

10、2057658092571-0.200899866713473-0.142057658092570(3.2)精确解yx1/42/43/41/40.1456064666070100.2059186398448590.1456064666070102/41.26088801585392e-171.78316493265431e-171.26088801585392e-173/4-0.145606466607010-0.205918639844859-0.145606466607010(3.3)数值解-精确解(即误差)y_.x1/42/43/41/4-0.00354880851443196-0.00

11、501877313137564-0.003548808514432732/42.15032106870631e-153.01985016524929e-152.11069424746919e-153/40.003548808514439730.005018773131386520.00354880851444026从而得到误差的范数为：1-范数：0.233770443573713; 2-范数：0.196807761631447;8-范数：0.327253314506086(3.4)图像(3.4.1 )数值解图像:(3.4.2 ) 精确解图像:( 5)主要程序( 5.1 )主程序%*%main_

12、chapter 主函数%信息10-2张道德%学号：10071223clccleara = 0; b=1;%裁值范围c=o； d=i；%曲值范围tfinal = 1;%最终时刻t=1/1600; %时间步长；h=1/40; %空间步长r=t/hA2;附比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);

13、 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 最终时刻

14、%t时间步长；%曲间步长r=t/hA2;附比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;fo

15、r i=1:n%* %An->n+1/2,u_xx向后差分，u_yy向前差分 for j=2:m for 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,:);%* %An->n+1,u_xx向前差分，u_yy向后差分 for k=2:mdd(

16、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一个%时方程常数项%用追赶法解三对角矩阵方程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%L分解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%追赶法实现%求解