微分方程数值解第二次上机报告Word文档下载推荐.docx
- 文档编号:22336383
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:27
- 大小:170.57KB
微分方程数值解第二次上机报告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)'
末时刻向前差分方法相对误差随时间步长变化对数-对数图'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 第二次 上机 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)