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

1、精选优质文档-倾情为你奉上 偏微分数值解法实验报告 实验名称：六点对称格式，ADI法，预校法，LOD法解二维抛物线方程的初值问题实验成员： 吴兴 杨敏 姚荣华 于潇龙 余凡 郑永亮实验日期：2013年5月17日指导老师：张建松一、 实验内容用六点对称格式，ADI法，预校法和LOD法求解二维抛物线方程的初值问题：已知（精确解为：）设差分解为，则边值条件为：初值条件为：取空间步长，时间步长网比。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） 第一步： ，构造出差分格式为：（2） 第二步：，构造出差分格式为：其中。假定第n层的已求得，则由求出，这只需按行解一些具有三对角系数矩阵的方程组；再由求出，这只需按列解一些具有三对角系数矩阵的方程组，所以计算时容易实现的。4:六点对称格式： 将向前差分格式和

3、向后差分格式做算术平均，即可以得到得六点对称格式：二、程序代码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/(h2);%网比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%差分格式方程组的系数矩阵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);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%构造t=1时精确解网格函数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);%生成网格采样点mesh(x_l,y_l,deta);title('误差网格分布');figure(2);mesh(x_l,y_l,jingquejie');%精确值的网格函数值title('精确解');figure(3);mes

7、h(x_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+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);%构造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构造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%在节点（xi，yj）=（i/4,j/4),j,k=1 2 3的计算结果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法主函数function=LOD()clearclcA=4(-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(-pi2*t/8);end%LOD差分函数%function=LOD_ch

15、afen(A,ax,bx,ay,by,t0,h,tao)ticNX=(bx-ax)/h;%x方向剖分份数NY=(by-ay)/h;%x方向剖分份数N=NX+1;Node=N2;%结点个数r=A*tao/(h2);%网比coefM=sparse(eye(Node);%系数矩阵R=sparse(zeros(Node,1);%不考虑边界条件时的系数矩阵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%循环求解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 %每一列开头和结尾为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%一维转化为二维并求真实解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%在节点（xi，yj）=（i/4,j/4),j,k=1 2 3的计算结果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%在节点（xj，yk）=（j/4,k/4),j,k=1 2 3的计算结果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('精确解与差分解的误差')%数值差分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轴分割的份数n=40;% 对y轴分割的份数tao=1/th;%时间步长knots=(m+1)*(n+1);% 结点个数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% 形成系数矩阵和右端项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;