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

1、ADI法求解二维抛物方程学校：中国石油大学(华东)学院：理学院姓名：张道德时间：2013.4.271、ADI 法介绍作为模型，考虑二维热传导方程的边值问题：Ut=UxxUyy,0:二X,y:：l,t0(3.6.1)0,作两族平行于坐标轴的网线：x=xj=jh,y=yk=kh,j,k=0,1,川，M,将区域0Mx,yMl分割成 M2个小矩形。第一个 ADI 算法(交替方向隐格式)是 Peacemarff 口 Rachford(1955)提出的。方法：由第 n 层到第 n+1 层计算分为两步：rn1,求uj,k2对u。向后差分，uyy向刖差分,构造出差分格式为:1on1!(,;%k2卜2xj，k式

2、为:23 力；？+6 刈：)h2xj,kyj,k1.其中j,k=1,lll,M1,n=0,1,2,111，上表n+一表示在t=t21假定第 n 层的 ujk已求得，则由(3.6.1)1求出 uj,k2,这只需按行(1)第1从n-n+,2,n2,nuj,k-uj,k(3.6.1)12nUn1n:：；1%枭一 2%;A 十h2u；k1-2u；,ku；,kh(2)第二步:1从n+-n+12n1，求uj,k2对u。向前差分,uyy向后差分，构造出差分格n1nuj,k-uj,k(361)22n：dn1nUj4-2Uj,k2,Uj2kh2n1n-1n1uj,k2%,Uj,kh1-j=(n+一2取值n-O2

3、L(j=1,111,M-1)解一些具有三对角系数矩阵的方程组；再由(3.6.1)2求出3 3k k1.-24,1.-24,njnju1n二rUj,k1321n1n(1-16r)Uj，k_32rUj,ku；,这只需按列(k=1,(|,M1)解一些具有三对角系数矩阵的方程组,所以计算时容易实现的。2、数值例子(1)问题用 ADI 法求解二维抛物方程的初边值问题：四=!(Uxx+Uyy),(x,y)wG=(0,1)M(0,1),t0,涿4JU(0,y,t)=U(1,y,t)=0,0y0,Uy(x,0,t)=Uy(x,1,t)=0,0 x0U(X,y,0)=sin二xcos-y.2已知(精确解为：U(

4、X,y,t)=sinnxcosnyexp(-彳 t)设xj=jh(j=0,1,|lJ),yk=kh(k=0,1,|,K),tn=nf(n=0,1,|,N)差分解为u；k,un,k=u；k=0,k=0,1 川,KUj,0=Uj,1,Uj,K1=Uj,K,j=0,1JII,J初值条件为：u；,k=sin 二 xjcos 二 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,|,K,已经知道第 0 列和第 K 列数值全为 0。1,、n1一(1)从

5、 n-n+-,求 ujk2对 uxx向后差分,u 向刖差分,构造出差分格式为：2J,xxyyn41n41n41n41Uj,k2-Uj,k1Uj1,k-2uj,k2Uj42kUj,k1-2uj,kUj,k42=w(h2+h2)1cn1=(-Ujk2,切。16 卜2xj,kyj,k从而得到：则边值条件为:132132116116十1 1/V/V+ +1 1j=1,2,|,J-1,k=1,2H，K1即按行用追赶法求解一系列下面的三对角方程组:从而得到:即按列用追赶法求解一系列下面的三对角方程组:nnUj,0-Uj,1=0nn一 Uj,KJUj,K其中(j=0,1,111,J)二 011+r161r3

6、2r3211r1611n2U1,k2r321一一r32,11r321r321r161r32一r321r161-r32r321+r16(J，)(Ja又根据边值条件得:nUj,0=Uj,1,Uj,Kn=Uj,KK行 U；K,(j=0,1/li,J)。(2)第二步:1,从n+n+1式为:j-ujk从而得到:32n11n1rUj,k1(1r)Uj,kn2U2,k2n；U3,k2n-UJ2,kn4UjfkJ封(J_L)1f2f3,j=0,1J|,J,解出第 0 行 U；,o和第,n-!,，、，、，一、，,、，求Uj,k2对UXX向前差分，Uyy向后差分，构造出差分格n-n-n:Uj比k-2小卜*Uj,k

7、+Uj,k由一2Uj,k+Uj,kh2h2216h2/2n22n1(xUj,k2.,yUj,k)1nL1一一rUjk_i=一rU32j,32n2j1,k1(1r)U16n21jk-rU32n-j_2,k,其中j=1,2UJ,-k1=HI1K,-又根据边值条件得：U；,0nnUj,1,Uj,K1U；,K,j=0,1,111,J,-1-11111r1+-rr7J47QOmQOu；n2J1V,J-IIJ1d)1u#_r1+r_rj，1321632n芈u21-11j,-r1+r一一r“n平321632Uj，3-ii-111u：九-r1+r-rj,K-321632.n+Uj,K2一f11f2f3f4fK

8、_3fK_2fKfKKJ111ujKr1+rr，_K4321632,14(3)求解结果(3.1)数值解yJ1/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.145606

9、4666070102/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.00501877313137564-0.003548808514432732/42.15032106870631e-153.01985016524929e-152.11069424746919e-153/40.0035488085

10、14439730.005018773131386520.00354880851444026从而得到误差的范数为：1-范数：0.233770443573713;2-范数：0.196807761631447;8-范数：0.327253314506086(3.4)图像(3.4.1)数值解图像:(3.4.2)精确解图像:(5)主要程序(5.1)主程序%*%main_chapter 主函数%百息 10-2 张道彳惠姆号：10071223clccleara=0;b=1;%裁值范围c=0;d=1;%曲值范围tfinal=1;%R 终时刻t=1/1600;%寸间步长；h=1/40;%空间步长r=t/hA2;附

11、比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_

12、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)%*functionu=ADI(a,b,c,d,t,h,tfinal)%(a,b)x 取值范围%(c,d)y 取值范围%tfinal 最终时刻%t 时间步长；%ha 同步长r=t/hA2;附比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

13、+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);end%修正第一项与最后一项，但由于第一项与最后一项均为零，可以省略%d(1)=d(1)+u1(j,1);d(m-

14、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%*(5.(3) “追赶法”程序%*%a 赶法functionx=zhuiganfa(a,b,c,d)%寸角线下方的元素，个数比 A一个%寸角线元素%寸角线上方的元素，个数比 A一个%时方程常数项%用追赶法解三对角矩阵方程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%L 分解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%追赶法