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

1偏微分数值解法实验报告实验名称：六点对称格式，ADI法，预校法，LOD法解二维抛物线方程的初值问题实验成员：吴兴杨敏姚荣华于潇龙余凡郑永亮实验日期：2013年5月17日指导老师：张建松一、实验内容用六点对称格式，ADI法，预校法和LOD法求解二维抛物线方程的初值问题：21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sincos.xxyyyyuuuxyGttuytuytytuxtuxtxtuxyxy已知（精确解为：2(,,)sincosexp()8uxytxyt）设(0,1,,),(0,1,,),(0,1,,)jknxjhjJykhkKtnnN差分解为,njku，则边值条件为：0,,,0,1,1,0,0,1,,,,0,1,,nnkJknnnnjjjKjKuukKuuuujJ初值条件为：0,sincosjkjkuxy取空间步长12140hhh，时间步长11600网比。1:ADI法：由第n层到第n+1层计算分成两步：先先第n层到n+1/2层，对uxx用向后差分逼近，对uyy用向前差分逼近，对uyy用向后差分逼近，于是2得到了如下格式：11112222,,1,,1,,1,,1221222,,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）（1）uu1111111222,,1,,1,,1,,12212212,,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhhuuuuuuuu（+）（2)uu其中j,k=1,2,…,M-1,n=0,1,2,…,上标n+1/2表示在t=tn+1/2+(n+1/2)取值。假定第n层的,njku已求得，则由（1）求出12,njku，再由（2）求出1,njku。2:预-校法差分格式：先通过U的n层求解U的n+1/4层，在通过U的n+1/4层求U的n的n+1/2层，最后通过U的n+1/2层求解U的n+1层，下为计算1n+2jku的预算格式：1n+2n4xjkjk211n+n+224yjkjk21n+n2jkjkn+1n1n+jkjk222xyjk2I-u=u2hI-u=u2huuu-u1=+u(5)h（）(3)（）(4)其次利用，构造更精确的校正格式即下（5)式（）3:LOD算法：由第n层到第n+1层计算分为两步：（1）第一步：12,12njkxxyyu从n-n+，求u对向后差分,u向前差分，构造出差分格3式为：1(3.6.1)11112222,,1,,1,1,,1,22122,,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkjkhhhuuuuuuuu（+）uu（2）第二步：12,12njkxxyyu从n+-n+1，求u对向前差分,u向后差分，构造出差分格式为：2(3.6.1)1111111222,,,1,,1,1,,1221212,,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnyjkjkhhhuuuuuuuu（+）uu其中1211,1,,1,0,1,2,,()22njkMnnn上表表示在t=t取值。假定第n层的,njku已求得，则由1(3.6.1)求出12,njku，这只需按行(1,,1)jM解一些具有三对角系数矩阵的方程组；再由2(3.6.1)求出1,njku，这只需按列(1,,1)kM解一些具有三对角系数矩阵的方程组，所以计算时容易实现的。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;4tao=1/T;%时间步长r=tao/(h^2);%网比a=1/16;U=ones(N+1,N+1);%迭代矩阵%按题意将边界点的值取为0forj=1:N+1U(1,j)=0;U(N+1,j)=0;end%初值条件fori=2:Nforj=1:N+1U(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);forn=1:T%计算到时间层t=1%x方向的迭代fork=2:Nforj=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)=A\f;end%y方向的迭代forj=2:Nfork=2:Nf2(k)=r*a/2*(U(j,k+1)+U(j,k-1))+(1-r*a)*U(j,k);endU(j,:)=(A2\f2)';end5end%构造t=1时精确解网格函数jingquejie=zeros(N+1,N+1);fori=1:N+1forj=1:N+1jingquejie(i,j)=sin(pi*x(i))*cos(pi*y(j))*exp(-pi^2/8);endenddeta=abs(U-jingquejie);%绝对误差deta_max=max(max(deta));fprintf('最大误差%f\n',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:预-校法差分格式：%用预-校法解抛物型方程clearclcformatlongJ=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);%真值%计算真值forn=1:N+1fori=1:J+1forj=1:J+1U1(i,j,n)=sin(pi*(i-1)*h)*cos(pi*(j-1)*h)*exp(-pi^2*(n-1)*t/8);end6endend%边值条件U在t=0层有U=sin（pi*x(i)）cos(pi*y(k))forj=1:J+1fork=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);fori=2:Ju(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);fori=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，按层forn=2:N+1u1=zeros((J+1)*(J+1),1);%构造u的1/4层fork=1:J+1%按行求fori=1:J+1b(1,i)=U(i,k,n-1);endy(1,1)=b(1,1);fori=2: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);7u1((J+1)*k,1)=x(1,J+1);fori=J:-1:1x(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);%g构造u的1/2层fork=1:J+1%按列求fori=1:J+1b1(1,i)=u1(k+(J+1)*(i-1),1);endy1(1,1)=b1(1,1);fori=2:J+1y1(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);fori=J:-1:1x1(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%求解Ufori=1:J+1%%边值条件：u（0,k,n）=u(J,k,n)=0,k=0,...,KU(1,i,n)=0;U(J+1,i,n)=0;endforj=2:J%按列求解fork=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,..,Jendend%在节点（xi，yj）=（i/4,j/4),j,k=123的计算结果fori=1:3forj=1:3UTRUE(i,j)=Ut(i*10+1,j*10+1);%精确解PrU(i,j)=UU(i*10+1,j*10+1);%lod差分解8endendErrors=PrU-UTRUE;%误差formatlongUTRUE'PrU'formatshortErrors3: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(