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

六点对称法-ADI法-预校法-和LOD法解二维抛物线方程精品文档 偏微分数值解法实验报告 实验名称：六点对称格式，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,上标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:六点对称格式： 将向前差分格式和向后差分格式做算术平均，即可以得到得六点对称格式：二、程序代码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)*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 %计算到时间层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(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);mesh(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)=sin(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(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,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),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(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 endend%在节点（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 longUTRUEPrUformat 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,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_chafen(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= sparse(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-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 %考虑边界条件 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 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) 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:六点对称格式%主程序clcclear%用六点对称格式求解二维抛物方程的初边值问题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的取值范围0x1;%y的取值范围0y1m=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 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) 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%初值条件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; XS(l,l+1)=-1;end%边界y=1for j=1:m+1 l=j*(n+1); XS(l,l)=1; XS(l,l-1)=-1;end%内点编号for j=2:m for k=2:n l=(j-1)*(n+1)+k; XS(l,l)=1+2*r; XS(l,l-1)=-r/2; XS(l,l+1)=-r/2; XS(l,l-(n+1)=-r/2; XS(l,l+(n+1)=-r/2; endend%开始迭代计算%for t=1:n1 for j=2:m for k=2:n l=(j-1)*(n+1)+k; Rhs(l)=r/2*(u0(j+1,k)+u0(j,k+1)+u0(j-1,k)+u0(j,k-1)+(1-2*r)*u0(j,k); end end solution = bicgstab(XS,Rhs,1.0e-6,400); % 将一维结果改成对应的二维 for j=1:m+1 for k=1:n+1 l=(j-1)*(n+1)+k; u(j,k)=solution(l); end end u0=u;end%精确解function f=uexact(x,y,t) f=sin(pi*x)*cos(pi*y)*exp(-pi*pi*t/8);end三、实验结果 四：实验心得 我没有直接参与编程，但是我知道如果紧紧一个人要在这一段时间内把这四个程序编完是要花费很多精力和时间的，而这时候就要靠团队协作了，每个同学都编写一个方法这样的话就会快很多，也能够在规定时间内完成任务。本次试验后我们会发现，这四种方法都是恒稳定的，通过以图形的方式比较它们的精确解和数值解之间的误差，可以直观的看出这四种方法的收敛性都很好。 五、实验分工吴兴 ： 预-校法差分格式编程者 学号：10072117杨敏 ： LOD法编程者 学号：10072118姚荣华: 整理报告 学号：10072119余凡 ：ADI法的编程者 学号：10072121郑永亮：六点对称法编程者 学号：10072122收集于网络，如有侵权请联系管理员删除