六点对称法,ADI法,预校法,和LOD法解二维抛物线方程

1、(1)偏微 分数值 解法实 验报告实验名称： 六点对称格式，ADI法，预校法，LOD法解二维抛物线方 程的初值问题实验成员：吴兴 杨敏 姚荣华 于潇龙 余凡 郑永亮实验日期：2013年5月17日指导老师：张建松实验内容用六点对称格式，ADI法，预校法和LOD法求解二维抛物线方程的初值问题:Uyy),(X,y)(0,1) (0,1),t0,u(0, y,t) u(1,y,t)0,0 y 1,t0,Uy(x,0, t)Uy(X,1,t) 0,0 X 1,t0u(x, y,0) sin xcos y.已知（精确解为：2u(x,y,t) sin xcos yexp(t)8设Xjjh(j 0,1,L ,

2、 J), ykkh(k 0,1 ,L ,K),tnn (n 0,1,L ,N)差分解为 u；,k，则边值条件为:u0,ku：,k0,k 0,1JL ,Knn nnUj,0 Uj,1,Uj,K 1 Uj,K,J 0,1 ,L ,J初值条件为：0 Uj,k sin Xj cos yk取空间步长h h1 h2 1 40，时间步长11600网比1: ADI 法:由第n层到第n+1层计算分成两步：先先第n层到n+1/2层，对Uxx用向后差分逼近, 得到了如下格式:1nU j ,k2 -U j ,k对uyy用向前差分逼近，对uyy用向后差分逼近,于是1n _(Uj 12k1 n _ 2Uj,k21n 空n

3、U j 1,k + U j ,k 12U：,kU j,k 1)1x2U；:U：,k)n 1 nu j,k -u j,k21 n _ Uj 12k丄 n ,k2Ujh2n2Ujn 1+ Uj,kn 12Uj,kn 1Uj,k 1)(2)1:u：F2 n 1、yUj,k)其中 j,k=1,2，,M-1,n=0,1,2,;上标 n+1/2 表示在 t=tn+1/2+（n+12）取值。假定第1n层的U；,k已求得，则由（1）求出u；/，再由（2）求出u；。2:预-校法差分格式：先通过U的；层求解U的n+1/4层，在通过U的n+1/4层求U的；1的n+1/2层，最后通过U的n+1/2层求解U的；+1层，

4、下为计算u：J的预算 格式：I-2hx2)1n+ u jk 4 =UnjkI-2hy2)Ujk1n+ 2 =u1 n+ 4 jkn+1 jk-unjk1 n+ x2 + y) U jk 2 (5)5)式n+ 其次利用 U；k， Ujk2构造更精确的校正格式即下（3：LOD算法：由第n层到第n+1层计算分为两步:(1)第一步：从 n-n+1，2求对向后差分，U yy向前差分，构造出差分格式为:1n -uj 12kn 12Uj,k2(361)1h2n 22nc nU j 1,k + Uj 1,k 2U j ,khnUj 1,k)(2)第二步：从n+丄-n+1，2式为:n 1 nUj,k -Uj,k

5、(361)2其中 j,k 1,L ,M假定第n层的丄h21nx(Uj,k2n 、Uj,k)1求Uj,k2对u心向前差分,uyy向后差分，构造出差分格1h2n 1n 1 n 1U j,k 1 2Uj,k Uj,k 1 +h21(u-k1 u；/)1nUj,k211n2Uj,k2h21nUj,k21)1 一1,n0,1,2, L ,上表 n - 表示在 t=t 12n _2 2(n12)取值。n 1Ujn,k已求得，则由(361)1求出uj,k2,这只需按行(j 1,L ,M 1)解一些具有三对角系数矩阵的方程组；再由(361)2求出U；，这只需按列(k 1L ,M 1)解一些具有三对角系数矩阵的

6、方程组, 所以计算时容易实现的。4：六点对称格式：将向前差分格式和向后差分格式做算术平均，即可以得到得六点对 称格式：、程序代码1： ADI 法：%交替方向差分格式 ADIclcx_a=0; x_b=1;%X勺区间端点 y_a=0; y_b=1;%y的区间端点N=40;%控制空间区域划分h=1/N;%空间步长x=x_a:h:x_b;y=y_a:h:y_b;T=1600;tao=1/T;%时间步长r=tao/(hA2);% 网比a=1/16;U=o nes(N+1,N+1);%迭代矩阵%按题意将边界点的值取为 0for j=1:N+1U(1,j)=0;U(N+1,j)=0;end%初值条件for

7、 i=2:Nfor j=1:N+1 U(i,j)=sin(pi*x(i)*cos(pi*y(j);endend%差分格式方程组的系数矩阵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);组装系数矩阵 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);f

8、2=zeros(N+1,1);for n=1:T %计算到时间层 t=1%x 方向的迭代for k=2:Nfor 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);endU(2:N,k)=Af;end%y 方向的迭代for j=2:Nfor k=2:Nf2(k)=r*a/2*(U(j,k+1)+U(j,k-1)+(1-r*a)*U(j,k);endU(j,:)=(A2f2);endend%构造 t=1 时精确解网格函数 jingquejie=zeros(N+1,N+1);for i=1

9、:N+1for j=1:N+1jingquejie(i,j)=sin(pi*x(i)*cos(pi*y(j)*exp(-piA2 /8); endend deta=abs(U-j in gquejie);%色对误差 deta_max=max(max(deta);fprintf( 最大误差 %fn,deta_max)figure(1);x,y=meshgrid(x);%生成网格采样点 mesh(x_l,y_l,deta);title(误差网格分布);figure(2);mesh(x_l,y_l,ji ngquejie);%精确值的网格函数值 title(精确解);figure(3);mesh(x

10、_l,y_l,U);%数值解的网格函数 title( 数值解 );U;2:预-校法差分格式 ：%用预-校法解抛物型方程clearclcformat longJ = 40;%x,y方向上的划分个数N = 1600;%t方向上的划分个数，这里只求到 t=1h=1/J;%x和y方向上的步长t=1/N;%t 方向上的步长r=1;%网格比a=1/16; %方程中的系数U=zeros(J+1,J+1,N+1);使用预-校法计算值U1=zeros(J+1,J+1,N+1);真值%计算真值for n=1:N+1for i=1:J+1for j=1:J+1U1(i,j, n)=si n( pi*(i-1)*h)

11、*cos(pi*(j-1)*h)*exp(-piA2* (n-1)*t/ 8);endendend%边值条件 U 在 t=0 层有 U=sin(pi*x(i)cos(pi*y(k)for j=1:J+1for k=1:J+1U(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:Ju(1,i)=1+a*r;endb = zeros(1,J+1);b1 = zeros(1,J+1);y =

12、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+1l(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+1u仁zeros(J+1)*(J+1),1);构造 u 的 1/4 层for k=1:J+1%按行求for i=1:J+1 b(1,i)=U(i,k,n-1);endy(1,1) = b(1,1);for i = 2 :

13、J+1y(1,i) = b(1,i) - l(1,i)*y(1,i - 1);endx(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);endendu2=zeros(J+1)*(J+1),1);%构造 u 的 1/2 层for k=1:J+1 %按列求 for i=1:J+1 b1(1,i)=u1(k+(J+1)*(i-1),1);endy1(1,1) = b1(1

14、,1);for i = 2 : J+1 y1(1,i) = b1(1,i) - l(1,i)*y1(1,i - 1);endx1(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);endend%求解 Ufor i=1:J+1 %边值条件： u(0,k,n)=u(J,k,n)=0,k=0,.,K U(1,i,n)=0;U(J+1,i,n)=0;endfor

15、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);endU(j,1,n)=U(j,2,n);%边值条件 u(j,0,n)=u(j,1,n),j=0,.,JU(j,J+1,n)=U(j,J,n); %边值条件 u(j,K-1,n)=u(j,K,n),j=0,.,J endend%在节点(xi, yj) = (i/4,j/4),j,k=1 2 3的

16、计算结果 for i=1:3for j=1:3UTRUE(i,j)=Ut(i*10+1,j*10+1);%精确解 PrU(i,j)=UU(i*10+1,j*10+1);%lod 差分解 endendErrors=PrU-UTRUE;误差 format long UTRUE PrU format shortErrors3丄0D算法%主%程序 %求解方程 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)%LOD法

17、主函数 function=LOD() clear clcA=4A(-2);%方程右边系数 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%真%实解函数 function fT=True(x,y,t) fT=sin(pi*x)*cos(pi*y)*exp(-piA2*t/ 8); end%LODi分函数 % function=LOD_chafen(A,ax

18、,bx,ay,by,t0,h,tao) ticNX=(bx-ax)/h;%x方向剖分份数 NY=(by-ay)/h;%x方向剖分份数 N=NX+1;Node=NA2;%结点个数r=A*tao/(hA2);% 网比coefM=sparse(eye(Node);% 系数矩阵R=sparse(zeros(Node,1);%不考虑边界条件时的系数矩阵for j=2:N-1for i=2:N-1k=i+(j-1)*N;coefM(k,k-1)=-r/ 2;%对角线左下coefM(k,k)=1+r;% 对角线coefM(k,k+1)=-r/ 2;%对角线右上endend%初始条件Mat= sparse(z

19、eros(Node,1);for i=1:Nfor j=1:N Mat(i-1)*N+j)=sin(pi*(i-1)*h)*cos(pi*(j-1)*h);endend%循环求解for m=1:10%第 n 层到 n+1/2 层for i=1:NcoefM(i,i+N)=-1;endfor i=Node-N+1:NodecoefM(i,i-N)=-1;end%每一列开头和结尾为 0for i=2:N-1coefM(1+(i-1)*N,2+(i-1)*N)=0;coefM(i*N,i*N-1)=0;endfor j=2:N-1for i=2:N-1R(i+(j-1)*N)=r/ 2*Mat(j+

20、(i-2)*N)+r/ 2*Mat(j+i*N)+(1-r)*Mat(j+(i-1)*N);endendMat,z=bicgstab(coefM,R,1e-6,100);%n+1/2 层到 n+1 层%右端项for i=2:N-1for j=2:N-12*Mat(i+j*N);R(j+(i-1)*N)=r/ 2*Mat(j-2)*N+i)+(1-r)*Mat(i+(j-1)*N)+r/ endend %将上次赋值取消 for i=1:NcoefM(i,i+N)=0;endfor i=Node-N+1:Node coefM(i,i-N)=0;end %考虑边界条件for i=2:N-1 coef

21、M(1+(i-1)*N,2+(i-1)*N)=-1; coefM(i*N,i*N-1)=-1;endMat,z=bicgstab(coefM,R,1e-6,100);end %一维转化为二维并求真实解 lod=sparse(zeros(N,N); true=sparse(zeros(N,N);for i=1:Nfor j=1:N lod(i,j)=Mat(j+(i-1)*N); true(i,j)=True(i-1)*h,(j-1)*h,tao*10);endend % %在节点( xi， yj) =(i/4,j/4),j,k=1 2 3 的计算结果 for i=1:3for j=1:3TRU

22、E(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) title(LOD 差分解 ) axis(ax bx ay by -1

23、 1)%画出精确解与差分解之间的误差图figure(2)surf(xx,yy,(true-lod)title( 精确解与差分解的误差 )tocend %精%确解 function true=TRUE(J,K,T) for n=1:T+1for j=1:J+1for k=1:K+1true(j,k,n )=si n( pi*(j-1)*h)*cos(pi*(k-1)*h)*exp(-piA2/8)*( n-1)*te);endendend4:六点对称格式 %主%程序 clc clear%用六点对称格式求解二维抛物方程的初边值问题ut=(42)*(uxx +uyy)%x轴的边值条件 u(0,y,t

24、)=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;b仁1; %x的取值范围0x1;%y的取值范围0y1 m=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+1for k=1:n+1uu(j,k)=uexact(j-1)*h,(k-1)*h,n1*tao); endend %在节点( xj，yk)=(j/4,k/4),j,

25、k=1 2 3的计算结果for j=1:3for k=1:3Exact(j,k)=uu(j*m/4,k*n/ 4);%精确解Differe ncesoluti on (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) title(精确解)%画出差分解图像figure(2) surf(xx,yy,u) title(差分解)%画出精确解与差分解之间的误差图figure(3) surf(

26、xx,yy,(uu-u)title(精确解与差分解的误差) %数%值%差分 function u=six(a,a1,b1,m,n,th,n1)%用六点对称格式求解二维抛物方程的初边值问题ut=(42)*(uxx +uyy)%x轴的边值条件 u(O,y,t)=u(1,y,t)=O%y轴的边值条件 uy(x,0,t)=uy(x,1,t)=0%初值条件 u(x,y,0)=sin(pi*x)*cos(pi*y) m=40;% 对 x 轴分割的份数n=40;% 对 y 轴分割的份数tao=1/th;%时间步长knots=(m+1)*(n+1);% 结点个数h=1/m;%空间步长r=tao/(a*a*h*

27、h);% 其中抛物方程中的 a uO=;%用于存储初值即第0层的值 for j=1:m+1for k=1:n+1 u0(j,k)=uexact(j-1)*h,(k-1)*h,0);endend% 形成系数矩阵和右端项 XS=sparse(eye(knots);Rhs=sparse(zeros(knots,1); solution=sparse(zeros(knots,1);%边界 x=0for j=1:n+1l=j;XS(l,l)=1;Rhs(l)=0;end%边界 x=1for j=1:n+1 l=m*(n+1)+j; XS(l,l)=1; Rhs(l)=0;end%边界 y=0for j=