微分方程数值解第二次上机报告Word文档下载推荐.docx
《微分方程数值解第二次上机报告Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《微分方程数值解第二次上机报告Word文档下载推荐.docx（27页珍藏版）》请在冰豆网上搜索。
逐渐降低，在边界
上满足零边界值条件，符合热传导物理规律。
从图2可以看出向前差分格式在空间步长
时间步长
时，有很好的收敛性。
末时刻差分方法解与精确解较为贴合。
经过多次放大后（如图3）才可以看出细微差别。
这是因为此时网格比
满足稳定性条件。
而五种方法原理不同，取不同空间步长、时间步长时的收敛性也不同。
下面研究差分方法的收敛阶，以深入研究收敛性。
3.收敛阶理论推导
将方程
在节点
离散化：
（1）
时间步长为
，空间步长为
.
*下以向前差分差分格式为例：
对充分光滑的解
，由Taylor展式：
（2）
（3）
（4）
对
（2）式移项得：
（5）
（3）、（4）式相加得：
（6）
（5）、（6）式带入
（1）式得：
（7）
其中：
（8）
（7）式舍去
即得到逼近
（1）式的向前格式差分方程：
（9）
其中，
从（8）式来看，网格化近似后余项
对时间步长
的局部误差阶为2，对空间步长
的局部误差阶为3.于是对时间、空间两种步长的整体误差阶为1和2.
4.收敛阶编程探讨
以向前差分格式为例，先固定时间步长，变动空间步长：
固定时间步长
，分别取四个空间步长
，这样取值是为了避免在右边界处数值计算时有时为
，有时取不到
，以影响末时刻结果精确性。
计算末时刻相对误差，误差与步长分别取对数后绘图如图4。
图4：
末时刻向前差分方法相对误差随空间步长变化对数-对视图
图中其实有三条直线，程序中用矩阵U存放了三条直线的斜率.分别约为：
1.403、2.135、2.056.符合理论讨论中的关于空间步长达到2阶收敛精度。
再固定空间步长，变动时间步长：
固定空间步长
，取两个空间步长
，再计算末时刻相对误差。
误差与步长分别取对数后绘图如图5。
同时计算了这条直线的斜率约为
符合理论讨论中的关于时间步长达到1阶收敛精度。
图5：
末时刻向前差分方法相对误差随时间步长变化对数-对视图
5.附录
所有Matlab程序如下：
forward.m文件：
clc
clear
l=1;
%参数l的值
a=1;
%系数a的值
tmax=0.1;
%t最大值
k=0.0002;
%时间t增量
h=0.02;
%x增量
s=a^
（2）*k/（h^2）;
%网格比
x=0:
h:
l;
t=0:
k:
tmax;
o=length（x）;
p=length（t）;
u=zeros（o,p）;
[X,T]=meshgrid（x,t）;
%u（x,0）初始层
fori=1:
o
ifx（i）<
=1/2
u（i,1）=2*x（i）;
else
u（i,1）=2-2*x（i）;
end
end
%u（0,t）和u（l,t）边界条件
u（1,:
）=0;
u（o,:
%向前差分
forj=1:
（p-1）
fori=2:
（o-1）
u（i,j+1）=s*u（i+1,j）+（1-2*s）*u（i,j）+s*u（i-1,j）;
%末时刻精确解
u_exact=zeros（60,o）;
fork=1:
60%取60项
u_exact（k,i）=8*sin（k*pi/2）*exp（-0.1*（k*pi）^2）*sin（k*pi*x（i））/（pi^2）;
u_e（i）=sum（u_exact（:
i））;
%解三维网格绘图
figure（）
mesh（X'
T'
u）
xlabel（'
x'
）;
ylabel（'
t'
zlabel（'
u（x,t）'
title（'
向前差分求解三维图'
）
%末时刻解比较
plot（x,u（:
p）,'
-'
）%差分方法求解
holdon
plot（x,u_e,'
-.'
）%精确解
末时刻向前差分方法解与精确解比较图'
legend（'
向前差分方法解'
'
精确解（n=60）'
backward.m文件：
k=0.00004;
h=0.01;
N=o-2;
%每一层内点个数
[m,n]=meshgrid（x,t）;
E=zeros（N,1）;
%边界不为零在此改动
F=zeros（N,1）;
%隐式方程组右端先声明
%隐式差分方程组系数矩阵A
A=zeros（N,N）;
A（1,1）=1+2*s;
A（1,2）=-s;
A（N,N）=1+2*s;
A（N,N-1）=-s;
fori=2:
N-1
A（i,i）=1+2*s;
A（i,i+1）=-s;
A（i,i-1）=-s;
%对时间层逐层求解
p-1
F=u（2:
N+1,j）+E;
%方程组右端更新
u（2:
N+1,j+1）=A\F;
mesh（m'
n'
隐式差分求解三维图'
）%隐式差分方法求解
末时刻隐式差分方法解与精确解比较图'
隐式差分方法解'
Crank_Nicolson.m文件：
s=a^
（2）*k/（2*h^2）;
%C-N网格比
%非齐次边界问题在此改动
%方程组隐式一端系数矩阵A
%方程组显式一端系数矩阵T
T=zeros（N,N）;
T（1,1）=1-2*s;
T（1,2）=s;
T（N,N）=1-2*s;
T（N,N-1）=s;
T（i,i）=1-2*s;
T（i,i+1）=s;
T（i,i-1）=s;
N+1,j+1）=inv（A）*T*F;
C-N差分求解三维图'
末时刻C-N差分方法解与精确解比较图'
C-N差分方法解'
AFB.m文件：
%一维交替显隐式格式
u=zeros（o,2*p-1）;
%共2p-1层前p层为初始层+隐式矫正层后p-1层为显式预测层
%交替显隐式格式网格比
%时间层循环
u（i,j+p）=s*u（i+1,j）+（1-2*s）*u（i,j）+s*u（i-1,j）;
%显式预估
N+1,j+p）+E;
%隐式矫正
u（:
1:
p））
交替显隐式格式求解三维图'
末时刻交替显隐式格式解与精确解比较图'
交替显隐式格式解'
skip.m文件：
%一维跳点格式
%跳点格式网格比
%跳点格式循环
ifrem（i+j+1,2）==0
%（i+j+1）先在偶数格点显式
ifrem（i+j+1,2）~=0
u（i,j+1）=u（i+1,j）/（1+2*s）+s*u（i+1,j+1）+s*u（i-1,j+1）;
%（i+j+1）再在奇数格点隐式（实际上显式）
跳点格式求解三维图'
）%跳点格式求解
末时刻跳点格式解与精确解比较图'
跳点格式解'
errforwardx.m文件:
formatlong
%向前差分格式收敛阶
h0=0.0125;
%初始空间步长
err=zeros（1,4）;
%存放相对误差
ke=1;
m=log（[1248]*h0）;
forhh=[1248]%対空间步长循环
h=hh*h0;
u_e=zeros（1,o）;
forii=1:
forkk=1:
u_exact（kk,ii）=8*sin（kk*pi/2）*exp（-0.1*（kk*pi）^2）*sin（kk*pi*x（ii））/（pi^2）;
u_e（ii）=sum（u_exact（:
ii））;
err（ke）=norm（u（:
p）-u_e'
2）/norm（u_e,2）;
%末时刻相对误差
ke=ke+1;
end;
%各段斜率计算
len=length（err）;
U=zeros（1,len-1）;
forflag=1:
len-1
U（flag）=（log（err（flag+1））-log（err（flag）））/（m（flag+1）-m（flag））;
%末时刻相对误差
plot（m,log（err）,'
）
log（h）'
log（error）'
末时刻向前差分方法相对误差随空间步长变化对数-对数图'
向前差分方法'
errforwardt.m文件:
%x增量
k0=0.00001;
%初始空间步长
err=zeros（1,2）;
%存放末时刻误差
k1=1;
m=log（[11.1]*k0）;
forkk=[11.1]
k=kk*k0;
%t增量
forke=1:
u_exact（ke,i）=8*sin（ke*pi/2）*exp（-0.1*（ke*pi）^2）*sin（ke*pi*x（i））/（pi^2）;
err（k1）=norm（u（:
k1=k1+1;
log（k）'
末时刻向前差分方法相对误差随时间步长变化对数-对数图'
升级会员 