六点对称法ADI法预校法和LOD法解二维抛物线方程资料.docx
- 文档编号:2280214
- 上传时间:2022-10-28
- 格式:DOCX
- 页数:21
- 大小:127.76KB
六点对称法ADI法预校法和LOD法解二维抛物线方程资料.docx
《六点对称法ADI法预校法和LOD法解二维抛物线方程资料.docx》由会员分享,可在线阅读,更多相关《六点对称法ADI法预校法和LOD法解二维抛物线方程资料.docx(21页珍藏版)》请在冰豆网上搜索。
六点对称法ADI法预校法和LOD法解二维抛物线方程资料
六点对称法-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
forj=1:
N+1
U(1,j)=0;
U(N+1,j)=0;
end
%初值条件
fori=2:
N
forj=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);
forn=1:
T%计算到时间层t=1
%x方向的迭代
fork=2:
N
forj=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方向的迭代
forj=2:
N
fork=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);
fori=1:
N+1
forj=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
formatlong
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);%真值
%计算真值
forn=1:
N+1
fori=1:
J+1
forj=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))
forj=1:
J+1
fork=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);
fori=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);
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);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求解U,求解到t=1,按层
forn=2:
N+1
u1=zeros((J+1)*(J+1),1);%构造u的1/4层
fork=1:
J+1%按行求
fori=1:
J+1
b(1,i)=U(i,k,n-1);
end
y(1,1)=b(1,1);
fori=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);
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);
end
end
u2=zeros((J+1)*(J+1),1);%g构造u的1/2层
fork=1:
J+1%按列求
fori=1:
J+1
b1(1,i)=u1(k+(J+1)*(i-1),1);
end
y1(1,1)=b1(1,1);
fori=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);
fori=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
fori=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
forj=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));
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=123的计算结果
fori=1:
3
forj=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;%误差
formatlong
UTRUE'
PrU'
formatshort
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
%%%%%%%%%%%真实解函数
functionfT=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(N
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 对称 ADI 法预校法 LOD 二维 抛物线 方程 资料