ADI(交替方向隱格式)求解二維拋物方程

1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)ADI法求解二維拋物方程學校：中國石油大學（華東） 學院：理學院 姓名：張道德 時間：2013.4.271、ADI法介紹作為模型，考慮二維熱傳導方程的邊值問題：（3.6.1）取空間步長，時間步長，作兩族平行于坐標軸的網(wǎng)線：將區(qū)域分割成個小矩形。第一個ADI算法（交替方向隱格式）是Peaceman和Rachford（1955）提出的。方法：由第n層到第n+1層計算分為兩步：第一步： ，構(gòu)造出差分格式為：第二步：，構(gòu)造出差分格式為：其中。假定第n層的已求得，則由求出，這只需

2、按行解一些具有三對角系數(shù)矩陣的方程組；再由求出，這只需按列解一些具有三對角系數(shù)矩陣的方程組，所以計算時容易實現(xiàn)的。2、數(shù)值例子（1）問題用ADI法求解二維拋物方程的初邊值問題：已知（精確解為：）設差分解為，則邊值條件為：初值條件為：取空間步長，時間步長網(wǎng)比。用ADI法分別計算到時間層。（2）計算過程根據(jù)邊值條件：，已經(jīng)知道第0列和第K列數(shù)值全為0。（1），構(gòu)造出差分格式為：從而得到：，其中即按行用追趕法求解一系列下面的三對角方程組：又根據(jù)邊值條件得：，解出第0行和第行。（2）第二步：，構(gòu)造出差分格式為：從而得到：，其中又根據(jù)邊值條件得：，從而得到：其中即按列用追趕法求解一系列下面的三對角方程組

3、：求解結(jié)果(3.1)數(shù)值解yx1/42/43/41/40.25780.34840.25782/42.484e-153.584e-152.773e-153/4-0.2571-0.3473-0.2570（3.2）精確解yx1/42/43/41/40.70100.48590.70102/41.392e-171.431e-171.392e-173/4-0.7010-0.4859-0.7010（3.3）數(shù)值解-精確解（即誤差）yx1/42/43/41/4-0.-0.-0.2/42.631e-153.929e-152.919e-153/40.0.0.從而得到誤差的范數(shù)為：范數(shù)：0.3713； 2-范數(shù)：0

4、.1447；-范數(shù)：0.6086（3.4）圖像（3.4.1）數(shù)值解圖像：精確解圖像：（5）主要程序（5.1）主程序%*%main_chapter主函數(shù)%信息10-2張道德%學號：clccleara = 0; b=1; %x取值范圍c=0; d=1; %y取值范圍tfinal = 1; %最終時刻t=1/1600;%時間步長；h=1/40;%空間步長r=t/h2;%網(wǎng)比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+1 u1(j,i) = uexact(x(i),y(j),1); endend%數(shù)值解u=AD

5、I(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_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,

6、c,d,t,h,tfinal ) %(a , b) x取值范圍 %(c, d) y取值范圍%tfinal最終時刻%t時間步長；%h空間步長r=t/h2;%網(wǎng)比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+1 u(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(

7、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:m for k=2:m d(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); end u2(1,:)=u2(2,:);

8、u2(m+1,:)=u2(m,:); %* %從n-n+1,u_xx向前差分，u_yy向后差分 for k=2:m dd(1)=0;dd(m+1)=0; for j=2:m dd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k)+u2(j,k); end u(:,k)=zhuiganfa(aa,bb,cc,dd); end %* u2=u;end%*（5.3）“追趕法”程序%*%追趕法function x=zhuiganfa(a,b,c,d)%對角線下方的元素，個數(shù)比A少一個% %對角線元素%對角線上方的元素，個數(shù)比A少一個%d為方程常數(shù)項%用追趕法解三對角矩陣

9、方程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%LU分解u(1)=b(1);for i=2:n l(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%追趕法實現(xiàn)%求解Ly=d，追的過程y(1)=d(1);for i=2:n y(i)=d(i)-l(i-1)*y(i-1);end%求解Ux=y,趕的過程x(n)=y(n)/u(n);for i=