ADI交替方向隱格式求解二維拋物方程含matlab程序

1、ADI法求解二維拋物方程學(xué)校：中國石油大學(xué)(華東)學(xué)院：理學(xué)院姓名：張道德時間：2013.4.271、ADI 法介紹作為模型，考慮二維熱傳導(dǎo)方程的邊值問題：Ut=UxxUyy,0:二X,y:：l,t0(3.6.1)0,作兩族平行于坐標軸的網(wǎng)線：x=xj=jh,y=yk=kh,j,k=0,1,川，M,將區(qū)域0Mx,yMl分割成 M2個小矩形。第一個 ADI 算法(交替方向隱格式)是 Peacemarff 口 Rachford(1955)提出的。方法：由第 n 層到第 n+1 層計算分為兩步：rn1,求uj,k2對u。向后差分，uyy向刖差分,構(gòu)造出差分格式為: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向后差分，構(gòu)造出差分格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)解一些具有三對角系數(shù)矩陣的方程組；再由(3.6.1)2求出3 3k k1.-24,1.-24,njnju1n二rUj,k1321n1n(1-16r)Uj，k_32rUj,ku；,這只需按列(k=1,(|,M1)解一些具有三對角系數(shù)矩陣的方程組,所以計算時容易實現(xiàn)的。2、數(shù)值例子(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)設(shè)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 網(wǎng)比r=dh2=1。用 ADI 法分別計算到時間層 t=1。(2)計算過程從n到n+1時，根據(jù)邊值條件：un,k=u：k=0,k=0,1,|,K,已經(jīng)知道第 0 列和第 K 列數(shù)值全為 0。1,、n1一(1)從

5、 n-n+-,求 ujk2對 uxx向后差分,u 向刖差分,構(gò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又根據(jù)邊值條件得: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向后差分，構(gòu)造出差分格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,-又根據(jù)邊值條件得：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)求解結(jié)果(3.1)數(shù)值解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)數(shù)值解-精確解(即誤差)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從而得到誤差的范數(shù)為：1-范數(shù)：0.233770443573713;2-范數(shù)：0.196807761631447;8-范數(shù)：0.327253314506086(3.4)圖像(3.4.1)數(shù)值解圖像:(3.4.2)精確解圖像:(5)主要程序(5.1)主程序%*%main_chapter 主函數(shù)%百息 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 函數(shù)%*%用 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)%寸角線下方的元素，個數(shù)比 A一個%寸角線元素%寸角線上方的元素，個數(shù)比 A一個%時方程常數(shù)項%用追趕法解三對角矩陣方程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%追趕法