六點對稱法-ADI法-預校法-和LOD法解二維拋物線方程

1、精選優(yōu)質(zhì)文檔-傾情為你奉上 偏微分數(shù)值解法實驗報告 實驗名稱：六點對稱格式，ADI法，預校法，LOD法解二維拋物線方程的初值問題實驗成員： 吳興 楊敏 姚榮華 于瀟龍 余凡 鄭永亮實驗日期：2013年5月17日指導老師：張建松一、 實驗內(nèi)容用六點對稱格式，ADI法，預校法和LOD法求解二維拋物線方程的初值問題：已知（精確解為：）設差分解為，則邊值條件為：初值條件為：取空間步長，時間步長網(wǎng)比。1: ADI法：由第n層到第n+1層計算分成兩步：先先第n層到n+1/2層，對uxx用向后差分逼近，對uyy用向前差分逼近，對uyy用向后差分逼近，于是得到了如下格式：其中j,k=1,2,M-1,n=0,1

2、,2,上標n+1/2表示在t=tn+1/2+(n+1/2)取值。假定第n層的已求得，則由（1）求出，再由（2）求出。2:預-校法差分格式： 先通過U的n層求解U的n+1/4層，在通過U的n+1/4層求U的n的n+1/2層，最后通過U的n+1/2層求解U的n+1層，下為計算的預算格式：3:LOD算法：由第n層到第n+1層計算分為兩步：（1） 第一步： ，構(gòu)造出差分格式為：（2） 第二步：，構(gòu)造出差分格式為：其中。假定第n層的已求得，則由求出，這只需按行解一些具有三對角系數(shù)矩陣的方程組；再由求出，這只需按列解一些具有三對角系數(shù)矩陣的方程組，所以計算時容易實現(xiàn)的。4:六點對稱格式： 將向前差分格式和

3、向后差分格式做算術(shù)平均，即可以得到得六點對稱格式：二、程序代碼1: ADI法：%交替方向差分格式 ADIclcx_a=0; x_b=1;%x的區(qū)間端點y_a=0; y_b=1;%y的區(qū)間端點N=40;%控制空間區(qū)域劃分h=1/N;%空間步長x=x_a:h:x_b;y=y_a:h:y_b;T=1600;tao=1/T;%時間步長r=tao/(h2);%網(wǎng)比a=1/16;U=ones(N+1,N+1);%迭代矩陣%按題意將邊界點的值取為0for j=1:N+1 U(1,j)=0; U(N+1,j)=0;end%初值條件for i=2:N for j=1:N+1 U(i,j)=sin(pi*x(i)

4、*cos(pi*y(j); endend%差分格式方程組的系數(shù)矩陣diag_0=(1+r*a)*ones(N-1,1);diag_1=(-r*a/2)*ones(N-2,1)'A=diag(diag_0)+diag(diag_1,1)+diag(diag_1,-1);%組裝系數(shù)矩陣A2=zeros(N+1);A2(2:N,2:N)=A;A2(1,1)=1;A2(N+1,N+1)=1;A2(1,2)=-1;A2(N+1,N)=-1;A2(2,1)=-r*a/2;A2(N,N+1)=-r*a/2;f=zeros(N-1,1);f2=zeros(N+1,1);for n=1:T %計算到時間

5、層t=1%x方向的迭代 for k=2:N for j=1:N-1 %邊界值為0，不必特殊處理j=1和N-1的情況 f(j)=r*a/2*(U(j,k)+U(j+2,k)+(1-r*a)*U(j+1,k); end U(2:N,k)=Af; end%y方向的迭代for j=2:N for k=2:N f2(k)=r*a/2*(U(j,k+1)+U(j,k-1)+(1-r*a)*U(j,k); end U(j,:)=(A2f2)'endend%構(gòu)造t=1時精確解網(wǎng)格函數(shù)jingquejie=zeros(N+1,N+1);for i=1:N+1 for j=1:N+1 jingquejie

6、(i,j)=sin(pi*x(i)*cos(pi*y(j)*exp(-pi2/8); endenddeta=abs(U-jingquejie);%絕對誤差deta_max=max(max(deta);fprintf('最大誤差%fn',deta_max)figure(1);x_l,y_l=meshgrid(x);%生成網(wǎng)格采樣點mesh(x_l,y_l,deta);title('誤差網(wǎng)格分布');figure(2);mesh(x_l,y_l,jingquejie');%精確值的網(wǎng)格函數(shù)值title('精確解');figure(3);mes

7、h(x_l,y_l,U');%數(shù)值解的網(wǎng)格函數(shù)title('數(shù)值解');U;2:預-校法差分格式：%用預-校法解拋物型方程clearclcformat longJ = 40;%x,y方向上的劃分個數(shù)N = 1600;%t方向上的劃分個數(shù)，這里只求到t=1h=1/J;%x和y方向上的步長t=1/N;%t方向上的步長r=1;%網(wǎng)格比a=1/16; %方程中的系數(shù)U=zeros(J+1,J+1,N+1);%使用預-校法計算值U1=zeros(J+1,J+1,N+1);%真值%計算真值for n=1:N+1 for i=1:J+1 for j=1:J+1 U1(i,j,n)=s

8、in(pi*(i-1)*h)*cos(pi*(j-1)*h)*exp(-pi2*(n-1)*t/8); end endend%邊值條件U在t=0層有U=sin（pi*x(i)）cos(pi*y(k)for j=1:J+1 for k=1:J+1 U(j,k,1)=sin(pi*(j-1)*h)*cos(pi*(k-1)*h); endend% U1(:,:,1)-U(:,:,1)%驗證初值條件%追趕法l=ones(1,J+1);l=l*(-a*r/2);v=l;u=ones(1,J+1);for i=2:J u(1,i)=1+a*r;endb = zeros(1,J+1);b1 = zeros

9、(1,J+1);y = zeros(1,J+1);x = zeros(1,J+1);y1 = zeros(1,J+1);x1 = zeros(1,J+1);u(1,1) = u(1,1);for i = 2 : J+1 l(1,i) = l(1,i)/u(1,i - 1); u(1,i) = u(1,i) - l(1,i)*v(1,i - 1);end%求解U,求解到t=1，按層for n=2:N+1 u1=zeros(J+1)*(J+1),1);%構(gòu)造u的1/4層 for k=1:J+1 %按行求 for i=1:J+1 b(1,i)=U(i,k,n-1); end y(1,1) = b(1

10、,1); for i = 2 : J+1 y(1,i) = b(1,i) - l(1,i)*y(1,i - 1); end x(1,J+1) = y(1,J+1)/u(1,J+1); u1(J+1)*k,1)=x(1,J+1); for i = J : -1 : 1 x(1,i)=(y(1,i) - v(1,i)*x(1,i + 1)/u(1,i); u1(J+1)*k-(J+1-i),1)=x(1,i); end end u2=zeros(J+1)*(J+1),1);%g構(gòu)造u的1/2層 for k=1:J+1 %按列求 for i=1:J+1 b1(1,i)=u1(k+(J+1)*(i-1

11、),1); end y1(1,1) = b1(1,1); for i = 2 : J+1 y1(1,i) = b1(1,i) - l(1,i)*y1(1,i - 1); end x1(1,J+1) = y1(1,J+1)/u(1,J+1); u2(J+1)*k,1)=x1(1,J+1); for i = J : -1 : 1 x1(1,i)=(y1(1,i) - v(1,i)*x1(1,i + 1)/u(1,i); u2(J+1)*k-(J+1-i),1)=x1(1,i); end end %求解U for i=1:J+1 %邊值條件：u（0,k,n）=u(J,k,n)=0,k=0,.,K U

12、(1,i,n)=0; U(J+1,i,n)=0; end for j=2:J %按列求解 for k=2:J %按行 U(j,k,n)=U(j,k,n-1)+r*a*(u2(k+(J+1)*j,1)+u2(k+(J+1)*(j-2),1)+u2(k+1+(J+1)*(j-1),1)+u2(k-1+(J+1)*(j-1),1)-4*u2(k+(J+1)*(j-1),1); end U(j,1,n)=U(j,2,n); %邊值條件u（j,0,n）=u(j,1,n),j=0,.,J U(j,J+1,n)=U(j,J,n); %邊值條件u（j,K-1,n）=u(j,K,n),j=0,.,J enden

13、d%在節(jié)點（xi，yj）=（i/4,j/4),j,k=1 2 3的計算結(jié)果for i=1:3 for j=1:3 UTRUE(i,j)=Ut(i*10+1,j*10+1);%精確解 PrU(i,j)=UU(i*10+1,j*10+1);%lod差分解 endendErrors=PrU-UTRUE;%誤差format longUTRUE'PrU'format shortErrors3:LOD算法%主程序%求解方程ut=(4（-2）)*(uxx +uyy)%x軸的邊值條件u(0,y,t)=u(1,y,t)=0%y軸的邊值條件uy(x,0,t)=uy(x,1,t)=0%初值條件u(x

14、,y,0)=sin(pi*x)*cos(pi*y) %LOD法主函數(shù)function=LOD()clearclcA=4(-2);%方程右邊系數(shù)ax=0;bx=1;%(ax，bx) x取值范圍ay=0;by=1;%(ay, by) y取值范圍t0=1;% (0,t0) 時間范圍h=1/40;%h 空間步長tao=1/1600;%t 時間步長LOD_chafen(A,ax,bx,ay,by,t0,h,tao)end%真實解函數(shù)function fT=True(x,y,t)fT=sin(pi*x)*cos(pi*y)*exp(-pi2*t/8);end%LOD差分函數(shù)%function=LOD_ch

15、afen(A,ax,bx,ay,by,t0,h,tao)ticNX=(bx-ax)/h;%x方向剖分份數(shù)NY=(by-ay)/h;%x方向剖分份數(shù)N=NX+1;Node=N2;%結(jié)點個數(shù)r=A*tao/(h2);%網(wǎng)比coefM=sparse(eye(Node);%系數(shù)矩陣R=sparse(zeros(Node,1);%不考慮邊界條件時的系數(shù)矩陣for j=2:N-1 for i=2:N-1 k=i+(j-1)*N; coefM(k,k-1)=-r/2;%對角線左下 coefM(k,k)=1+r;%對角線 coefM(k,k+1)=-r/2;%對角線右上 endend%初始條件Mat= spa

16、rse(zeros(Node,1);for i=1:N for j=1:N Mat(i-1)*N+j)=sin(pi*(i-1)*h)*cos(pi*(j-1)*h); endend%循環(huán)求解for m=1:10 %第n層到n+1/2層 for i=1:N coefM(i,i+N)=-1; end for i=Node-N+1:Node coefM(i,i-N)=-1; end %每一列開頭和結(jié)尾為0 for i=2:N-1 coefM(1+(i-1)*N,2+(i-1)*N)=0; coefM(i*N,i*N-1)=0; end for j=2:N-1 for i=2:N-1 R(i+(j-

17、1)*N)=r/2*Mat(j+(i-2)*N)+r/2*Mat(j+i*N)+(1-r)*Mat(j+(i-1)*N); end end Mat,z=bicgstab(coefM,R,1e-6,100); %n+1/2層到n+1層 %右端項 for i=2:N-1 for j=2:N-1 R(j+(i-1)*N)=r/2*Mat(j-2)*N+i)+(1-r)*Mat(i+(j-1)*N)+r/2*Mat(i+j*N); end end %將上次賦值取消 for i=1:N coefM(i,i+N)=0; end for i=Node-N+1:Node coefM(i,i-N)=0; end

18、 %考慮邊界條件 for i=2:N-1 coefM(1+(i-1)*N,2+(i-1)*N)=-1; coefM(i*N,i*N-1)=-1; end Mat,z=bicgstab(coefM,R,1e-6,100);end%一維轉(zhuǎn)化為二維并求真實解lod=sparse(zeros(N,N);true=sparse(zeros(N,N);for i=1:N for j=1:N lod(i,j)=Mat(j+(i-1)*N); true(i,j)=True(i-1)*h,(j-1)*h,tao*10); endend%在節(jié)點（xi，yj）=（i/4,j/4),j,k=1 2 3的計算結(jié)果for

19、 i=1:3 for j=1:3 TRUE(i,j)=true(i*NX/4,j*NX/4);%精確解 LOD(i,j)=lod(i*NX/4,j*NX/4);%差分解 endenderror=LOD-TRUE;%誤差 % 作圖 x=ax:h:bx; y=ay:h:by; xx,yy=meshgrid(x,y); %畫出精確解圖像 figure(1) subplot(1,2,1); surf(xx,yy,true') title('精確解') axis(ax bx ay by -1 1) %畫出差分解圖像 subplot(1,2,2); surf(xx,yy,lod&#

20、39;) title('LOD差分解') axis(ax bx ay by -1 1) %畫出精確解與差分解之間的誤差圖 figure(2) surf(xx,yy,(true-lod)') title('精確解與差分解的誤差') tocend%精確解function true=TRUE(J,K,T)for n=1:T+1 for j=1:J+1 for k=1:K+1 true(j,k,n)=sin(pi*(j-1)*h)*cos(pi*(k-1)*h)*exp(-pi2/8)*(n-1)*te); end endend4:六點對稱格式%主程序clccl

21、ear%用六點對稱格式求解二維拋物方程的初邊值問題ut=(4-2)*(uxx +uyy)%x軸的邊值條件u(0,y,t)=u(1,y,t)=0%y軸的邊值條件uy(x,0,t)=uy(x,1,t)=0%初值條件u(x,y,0)=sin(pi*x)*cos(pi*y)a1=0;b1=1; %x的取值范圍0<x<1;%y的取值范圍0<y<1m=40;n=40;th=1600;tao=1/th;a=4;n1=input('輸入要計算的時間層n1=');u=six(a,a1,b1,m,n,th,n1);%計算精確解h=(b1-a1)/m; for j=1:m+1

22、 for k=1:n+1 uu(j,k)=uexact(j-1)*h,(k-1)*h,n1*tao); endend%在節(jié)點（xj，yk）=（j/4,k/4),j,k=1 2 3的計算結(jié)果for j=1:3 for k=1:3 Exact(j,k)=uu(j*m/4,k*n/4);%精確解 Differencesolution(j,k)=u(j*m/4,k*n/4);%差分解 endendExactDifferencesolution % 作圖 x=a1:h:b1; y=a1:h:b1; xx,yy=meshgrid(x,y); %畫出精確解圖像 figure(1) surf(xx,yy,uu

23、') title('精確解') %畫出差分解圖像 figure(2) surf(xx,yy,u') title('差分解') %畫出精確解與差分解之間的誤差圖 figure(3) surf(xx,yy,(uu-u)') title('精確解與差分解的誤差')%數(shù)值差分function u=six(a,a1,b1,m,n,th,n1)%用六點對稱格式求解二維拋物方程的初邊值問題ut=(4-2)*(uxx +uyy)%x軸的邊值條件u(0,y,t)=u(1,y,t)=0%y軸的邊值條件uy(x,0,t)=uy(x,1,t)=0

24、%初值條件u(x,y,0)=sin(pi*x)*cos(pi*y)m=40;% 對x軸分割的份數(shù)n=40;% 對y軸分割的份數(shù)tao=1/th;%時間步長knots=(m+1)*(n+1);% 結(jié)點個數(shù)h=1/m; %空間步長r=tao/(a*a*h*h);%其中拋物方程中的au0=;%用于存儲初值即第0層的值for j=1:m+1 for k=1:n+1 u0(j,k)=uexact(j-1)*h,(k-1)*h,0); endend% 形成系數(shù)矩陣和右端項XS=sparse(eye(knots);Rhs=sparse(zeros(knots,1);solution=sparse(zeros(knots,1);%邊界x=0 for j=1:n+1 l=j; XS(l,l)=1; Rhs(l)=0; end %邊界x=1 for j=1:n+1 l=m*(n+1)+j; XS(l,l)=1; Rhs(l)=0; end %邊界y=0for j=1:m+1 l=(j-1)*(n+1)+1; XS(l,l)=1;