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

ADI法求解二维抛物方程学校：中国石油大学（华东）学院：理学院姓名：张道德时间：2013.4.271、ADI法介绍作为模型，考虑二维热传导方程的边值问题：（3.6.1）,0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0txxyyuuuxyltuxyxyuytulytuxtuxlt取空间步长1hM=，时间步长0t，作两族平行于坐标轴的网线：,,,0,1,,,jkxxjhyykhjkM将区域0,xyl分割成2M个小矩形。第一个ADI算法（交替方向隐格式）是Peaceman和Rachford（1955）提出的。方法：由第n层到第n+1层计算分为两步：（1）第一步：12,12njkxxyyu从n-n+，求u对向后差分,u向前差分，构造出差分格式为：1(3.6.1)11112222,,1,,1,,1,,1221222,,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）uu（2）第二步：12,12njkxxyyu从n+-n+1，求u对向前差分,u向后差分，构造出差分格式为：2(3.6.1)1111111222,,1,,1,,1,,12212212,,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）uu其中1211,1,,1,0,1,2,,()22njkMnnn上表表示在t=t取值。假定第n层的,njku已求得，则由1(3.6.1)求出12,njku，这只需按行(1,,1)jM解一些具有三对角系数矩阵的方程组；再由2(3.6.1)求出1,njku，这只需按列(1,,1)kM解一些具有三对角系数矩阵的方程组，所以计算时容易实现的。2、数值例子（1）问题用ADI法求解二维抛物方程的初边值问题：21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sincos.xxyyyyuuuxyGttuytuytytuxtuxtxtuxyxy已知（精确解为：2(,,)sincosexp()8uxytxyt）设(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。（2）计算过程从n到n+1时，根据边值条件：0,,0,0,1,,nnkJkuukK，已经知道第0列和第K列数值全为0。（1）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即按行用追赶法求解一系列下面的三对角方程组：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。（2）第二步：12,12njkxxyyu从n+-n+1，求u对向前差分,u向后差分，构造出差分格式为：1111111222,,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即按列用追赶法求解一系列下面的三对角方程组：1,01,11,21,31,31,21111113216321111321632111132163211113216321111321632111132163211njnjnjnjnjKnjKKKrrrrrrrrrrrrrrrrrruuuuuuu12343211,111,1KKnKjKKKnjKKffffffffu（3）求解结果(3.1)数值解yx1/42/43/41/40.1420576580925780.2008998667134840.1420576580925782/42.16292994886484e-153.03768181457584e-152.12330312762773e-153/4-0.142057658092571-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）数值解-精确解（即误差）yx1/42/43/41/4-0.00354880851443196-0.00501877313137564-0.003548808514432732/42.15032106870631e-153.01985016524929e-152.11069424746919e-153/40.003548808514439730.005018773131386520.00354880851444026从而得到误差的范数为：1-范数：0.233770443573713；2-范数：0.196807761631447；∞-范数：0.327253314506086（3.4）图像（3.4.1）数值解图像：（3.4.2）精确解图像：00.20.40.60.8100.51-0.4-0.200.20.4x轴图一、数值解图像y轴t轴00.20.40.60.8100.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/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)%误差分析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);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