六点对称法ADI法预校法和LOD法解二维抛物线方程资料Word文档格式.docx
《六点对称法ADI法预校法和LOD法解二维抛物线方程资料Word文档格式.docx》由会员分享，可在线阅读，更多相关《六点对称法ADI法预校法和LOD法解二维抛物线方程资料Word文档格式.docx（21页珍藏版）》请在冰豆网上搜索。
先通过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:
六点对称格式：
将向前差分格式和向后差分格式做算术平均，即可以得到得六点对称格式：
二、程序代码
%交替方向差分格式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:
y_b];
T=1600;
tao=1/T;
%时间步长
r=tao/（h^2）;
%网比
a=1/16;
U=ones（N+1,N+1）;
%迭代矩阵
%按题意将边界点的值取为0
forj=1:
N+1
U（1,j）=0;
U（N+1,j）=0;
end
%初值条件
fori=2:
N
forj=1:
U（i,j）=sin（pi*x（i））*cos（pi*y（j））;
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）;
forn=1:
T%计算到时间层t=1
%x方向的迭代
fork=2:
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）;
U（2:
N,k）=A\f;
%y方向的迭代
forj=2:
f2（k）=r*a/2*（U（j,k+1）+U（j,k-1））+（1-r*a）*U（j,k）;
U（j,:
）=（A2\f2）'
%构造t=1时精确解网格函数
jingquejie=zeros（N+1,N+1）;
fori=1:
jingquejie（i,j）=sin（pi*x（i））*cos（pi*y（j））*exp（-pi^2/8）;
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'
%精确值的网格函数值
精确解'
figure（3）;
mesh（x_l,y_l,U'
%数值解的网格函数
数值解'
U;
%用预-校法解抛物型方程
clear
formatlong
J=40;
%x,y方向上的划分个数
N=1600;
%t方向上的划分个数，这里只求到t=1
h=1/J;
%x和y方向上的步长
t=1/N;
%t方向上的步长
r=1;
%网格比
%方程中的系数
[U]=zeros（J+1,J+1,N+1）;
%使用预-校法计算值
[U1]=zeros（J+1,J+1,N+1）;
%真值
%计算真值
fori=1:
J+1
U1（i,j,n）=sin（pi*（i-1）*h）*cos（pi*（j-1）*h）*exp（-pi^2*（n-1）*t/8）;
%边值条件U在t=0层有U=sin（pi*x（i））cos（pi*y（k））
fork=1:
U（j,k,1）=sin（pi*（（j-1）*h））*cos（pi*（（k-1）*h））;
%U1（:
:
1）-U（:
1）%验证初值条件
%追赶法
l=ones（1,J+1）;
l=l*（-a*r/2）;
v=l;
u=ones（1,J+1）;
J
u（1,i）=1+a*r;
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）;
fori=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）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求解U,求解到t=1，按层
forn=2:
u1=zeros（（J+1）*（J+1）,1）;
%构造u的1/4层
J+1%按行求
b（1,i）=U（i,k,n-1）;
y（1,1）=b（1,1）;
fori=2:
y（1,i）=b（1,i）-l（1,i）*y（1,i-1）;
x（1,J+1）=y（1,J+1）/u（1,J+1）;
u1（（J+1）*k,1）=x（1,J+1）;
fori=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）;
u2=zeros（（J+1）*（J+1）,1）;
%g构造u的1/2层
J+1%按列求
b1（1,i）=u1（k+（J+1）*（i-1）,1）;
y1（1,1）=b1（1,1）;
y1（1,i）=b1（1,i）-l（1,i）*y1（1,i-1）;
x1（1,J+1）=y1（1,J+1）/u（1,J+1）;
u2（（J+1）*k,1）=x1（1,J+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）;
%求解U
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;
forj=2:
J%按列求解
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））;
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
%在节点（xi，yj）=（i/4,j/4）,j,k=123的计算结果
3
UTRUE（i,j）=Ut（i*10+1,j*10+1）;
%精确解
PrU（i,j）=UU（i*10+1,j*10+1）;
%lod差分解
Errors=PrU-UTRUE;
%误差
UTRUE'
PrU'
formatshort
Errors
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（）
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）
%%%%%%%%%%%真实解函数
functionfT=True（x,y,t）
fT=sin（pi*x）*cos（pi*y）*exp（-pi^2*t/8）;
%%%%%LOD差分函数%%%%%
function[]=LOD_chafen（A,ax,bx,ay,by,t0,h,tao）
tic
NX=（bx-ax）/h;
%x方向剖分份数
NY=（by-ay）/h;
N=NX+1;
Node=N^2;
%结点个数
r=A*tao/（h^2）;
coefM=sparse（eye（Node））;
%系数矩阵
R=sparse（zeros（N
升级会员 