六点对称法-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法：
%交替方向差分格式 ADI
clc
x_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/(h^2);%网比
a=1/16;
U=ones(N+1,N+1);%迭代矩阵
%按题意将边界点的值取为0
for 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));
end
end
%差分格式方程组的系数矩阵
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)=A\f;
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,:)=(A2\f2)';
end
end
%构造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(-pi^2/8);
end
end
deta=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:预-校法差分格式：
%用预-校法解抛物型方程
clear
clc
format long

J = 40;%x,y方向上的划分个数
N = 1600;%t方向上的划分个数，这里只求到t=1
h=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(-pi^2*(n-1)*t/8);
end
end
end

%边值条件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));
end
end
% 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;
end
b = 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
end
end
%在节点（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差分解
end
end
Errors=PrU-UTRUE;%误差
format long
UTRUE'
PrU'
format short
Errors
3: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()
clear
clc

A=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(-pi^2*t/8);
end

%%%%%LOD差分函数%%%%%
function[]=LOD_chafen(A,ax,bx,ay,by,t0,h,tao)
tic
NX=(bx-ax)/h;%x方向剖分份数
NY=(by-ay)/h;%x方向剖分份数
N=NX+1;
Node=N^2;%结点个数
r=A*tao/(h^2);%网比

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;%对角线右上
end
end
%初始条件
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);
end
end
%循环求解
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);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%在节点（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);%差分解
end
end
error=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('精确解与差分解的误差')
toc
end
%%%%%%%%%精确解
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((-pi^2/8)*(n-1)*te);
end
end
end
4:六点对称格式
%%%%%%%%%%%%%主程序
clc
clear
%用六点对称格式求解二维抛物方程的初边值问题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<1
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+1
for k=1:n+1
uu(j,k)=uexact((j-1)*h,(k-1)*h,n1*tao);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%在节点（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);%差分解
end
end
Exact
Differencesolution

% 作图
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);%其中抛物方程中的a
u0=[];%用于存储初值即第0层的值
for j=1:m+1
for k=1:n+1
u0(j,k)=uexact((j-1)*h,(k-1)*h,0);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 形成系数矩阵和右端项
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=0
for j=1:m+1
l=(j-1)*(n+1)+1;
XS(l,l)=1;
XS(l,l+1)=-1;
end
%边界y=1
for 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;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%开始迭代计算%%
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
收集于网络，如有侵权请联系管理员删除