偏微分方程数值解法Word文档下载推荐.doc
《偏微分方程数值解法Word文档下载推荐.doc》由会员分享，可在线阅读，更多相关《偏微分方程数值解法Word文档下载推荐.doc（12页珍藏版）》请在冰点文库上搜索。
用六点对称差分格式计算如下热传导方程的初边值问题
已知其精确解为
二、理论
1．考虑的问题
考虑一维模型热传导方程
（1.1），
其中为常数。
是给定的连续函数。
（1.1）的定解问题分两类：
第一，初值问题（Cauchy问题）：
求足够光滑的函数，满足方程（1.1）和初始条件：
（1.2）,
第二，初边值问题（也称混合问题）：
及边值条件
，
假定和在相应的区域光滑，并且于，两点满足相容条件,则上述问题有唯一的充分光滑的解。
现在考虑边值问题（1.1），（1.3）的差分逼近
取为空间步长，为时间步长，其中，是自然数，
；
将矩形域分割成矩形网格。
其中表示网格节点；
表示网格内点（位于开矩形中的网格节点）的集合；
表示位于闭矩形中的网格节点的集合；
表示-网格边界点的集合。
表示定义在网点处的待求近似解，，。
注意到在节点处的微商和差商之间的下列关系（）：
2．区域网格剖分
取空间步长和时间步长，其中都是正整数。
用两族平行直线和将矩形域分割成矩形网格，网格节点为。
以表示网格内点集合，即位于开矩形的网点集合；
表示所有位于闭矩形的网点集合；
是网格界点集合。
其次，用表示定义在网点的函数，
3．建立相应差分格式
数值分析中，Crank-Nicolson方法是有限差分方法中的一种，用于数值求解热方程以及形式类似的偏微分方程。
它在时间方向上是隐式的二阶方法，数值稳定。
该方法诞生于20世纪，由JohnCrank与PhyllisNicolson发展。
①向前差分格式
，==0
②向后差分格式
将向前差分格式和向后差分格式做算术平均，得到的差分格式称之为六点对称格式，也称为Grank-Nicholson格式：
，==0
进一步，
+=++
按层计算：
首先，取，则利用初值和边值==0，来确定出第一层的，，即求解方程组：
+=++
，==0。
求出，在由,取，可利用，解出，。
如此下去，即可逐层算出所有，。
若记
三、截断误差
=。
注意：
=
又
两式相加
而
+
故有
。
四、稳定性与收敛性
抛物方程的两层差分格式可以统一写成向量形式:
（1.7）
其中，和是阶矩阵。
我们假定可逆，即（1.7）是唯一可解的。
对于显格式，等于单位矩阵。
三层格式可以通过引入新变量化成两层格式。
假设差分解的初始值（其实可以是任一层的值）有误差，以后各层计算没有误差，让我们来考察初始误差对以后各层的影响。
令和分别是以和为初始值由差分格式（1.7）得到的两组差分解，则满足
（1.8）
因此，按初值稳定应该意味着。
这就导致如下定义：
假设，我们称差分格式（2.1）按初值稳定，如果存在正常数和，使得以下不等式成立：
（1.8），
这里是上的某一个范数，例如
类似地，假设，我们称差分格式（2.1）按右端稳定，如果存在正常数和，使得以下不等式成立：
（1.8），
可以证明，差分格式若按初值稳定，则一定按右端稳定。
因此，这时我们简单地称差分格式稳定。
前面讨论的向前差分格式①当网比时稳定，当时不稳定。
这就意味着给定空间步长以后，时间步长必须足够小，才能保证稳定。
而向后差分格式②和Grank-Nicholson格式（1.6）则对任何网比都是稳定的，时间步长可以取得大一些，从而提高运算效率。
如果某个差分格式的截断误差当和趋于0时随之趋于0，则称这个差分格式是相容的。
可以证明：
若差分格式是相容的和稳定的，则它是收敛的，并且差分解与微分解之间误差的阶等于截断误差的阶。
因此，对任何网比，向后差分格式（1.6）有收敛阶。
五、结论
对于扩散方程（包括许多其他方程），可以证明Crank-Nicolson方法无条件稳定。
但是，如果时间步长与空间步长平方的比值过大（一般地，大于1/2），近似解中将存在虚假的振荡或衰减。
基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉方法进行计算，这样即可以保证稳定，又避免了解的伪振荡。
由本题可以总结出，抛物型方程的六点对称差分法所得的数值解能够较好地逼近方程的精确解，且区域剖分得越细，即步长越小，数值解与精确解的误差就越小，数值解越逼近精确解。
六、附录
取，，则，满足稳定性条件
另取，，则，亦满足稳定性条件
程序代码
formatlong
a=2;
l=1;
T=1;
N=5;
M=100;
h=l/N;
to=T/M;
r=（a*to）/h^2;
forj=1:
N+1
x（j）=（j-1）*h;
fork=1:
M+1
t（k）=（k-1）*to;
u（j,k）=exp（x（j）+2*t（k））;
end
end
u%求解精确解
us（j,1）=exp（x（j））;
fork=1:
us（1,k）=exp（2*t（k））;
us（N+1,k）=exp（1+2*t（k））;
fork=2:
forj=2:
N
us（j,k）=r*us（j-1,k-1）+（1-2*r）*u（j,k-1）+r*us（j+1,k-1）;
end
us%求解数值解
forj=1:
R（j,k）=abs（u（j,k）-us（j,k））;
R%计算误差
Rmax=max（max（R））%求误差的最大值
精确解与数值解的比较：
x=0:
0.1:
1;
holdon
plot（x,u（:
M+1）,'
b'
）;
plot（x,us（:
y'
title（'
t=1，h=1/10，τ=1/400时精确解和数值解的比较'
）
text（0.05,21,'
蓝：
精确解'
text（0.05,20,'
黄：
数值解'
holdoff
取不同步长时的误差比较：
1/5:
y=0:
1/10:
z=0:
1/20:
plot（x,R（:
N分别取5,10，20
①当N=5时，
误差最大值：
Rmax=0.002346916056785
②当N=10时，
Rmax=5.911936126903328e-004
③当N=20时，
升级会员 