求解热传导方程的crank_nicolson方法

《求解热传导方程的crank_nicolson方法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、2012年10月枣庄学院学报Oct．2012第29卷第5期JOURNALOFZAOZHUANGUNIVERSITYVol．29NO．5求解热传导方程的Crank－Nicolson方法陶燕燕(青岛科技大学数理学院，山东青岛266061)22［摘要］给出了数值求解热传导方程的一种Crank－Nicolson格式，其截断误差为O(τ+h)，并且分析了该差分格式的稳定性．在最后的数值例子中，验证了该格式求解出的数值解可以很好的逼近精确解，以及当空间步长和时间步长11同时缩小倍时，最大误差约缩小为原来的．24［关键词］热传导方程;差分方法;Crank－Nicolson

2、格式;无条件稳定［中图分类号］O175．26［文献标识码］A［文章编号］1004－7077(2012)05－0004－050引言求解热传导方程的差分方法有古典显式格式、古典隐式格式和Crank－Nicolson格式，单纯地从并行计算层面来看，古典显式格式计算量小，便于编程，但是这种格式是条1件稳定的(r≤，其中r是步长比);而古典隐式格式和Crank－Nicolson格式是绝对稳定2的，但是需要解整体的线性代数方程组，不能直接实现并行运算．有限差分区域分解算法综合了显格式和隐格式的优点，是一种高效实用的方法．目前已有许多关于区域分解方法的结果，文献［2］的内

3、边界点实用大空间步长H=mh的古典显格式，在每个子区域内实用古典隐格式求解，发展了有限差分区域分解算法，文献［3］用Saul＇yev非对称格式和全隐格式导出了一种有限差分格式，它的稳定性为约束条件是r≤1．本文主要研究热传导方程的Crank－Nicolson差分格式，并验证了格式的无条件稳定性，在最后的数值试验中充分验证了这一点．1Crank－Nicolson差分格式的构造给出要解决的问题2uu=0＜x＜l，0＜t≤Tt2x(1){u(x，0)=φ(x)0≤x≤lu(0，t)=1(t)u(1，t)=2(t)0＜t≤T为了对初边值问题(1)作差分逼

4、近，首先对于求解区域作网格剖分，取空间步长h=lTn和时间步长τ=，其中J、N为自然数，x=xj=jh(j=1，2…J)Jh=l，t=t=nτ(nJN=1，2…N)Nτ=T．将求解区域{0≤x≤l，0≤t≤T}分割成矩形网格，网格节点为nnnnn(xj，t)，记Uj=u(xj，t)为热传导方程的真解，uj为热传导方程(1)在网格节点(xj，t)τ处的数值解，r=2为步长比．h［收稿日期］2012－05－01［作者简介］陶燕燕(1987－)，女，青岛科技大学数理学院硕士研究生，研究方向:应用数学．·4·陶燕燕求解热传导方程的Crank－Nicolson方法方程

5、(1)在(xj，tn+1)处满足关系式21n+1nn+2u－uujj2{}=+O(τ)(2)xjτ12n+u2=12(un+1+un)+O(τ2+h2){2}2δxjjxj2hn+1n+1n+1nnn1uj+1－2uj+uj－1uj+1－2uj+uj－1=(+)(3)222hh将(2)、(3)式代入(1)舍去误差项得(4)式12n+112n(1－rδx)uj=(1+rδx)uj(4)2222上式(4)被称为Crank－Nicolson差分格式，截断误差为O(τ+h)n+11n+1n+1也可将(4)式改写为(1+r)uj－r(uj+1+uj－1)2n1

6、nn=(1－r)uj+r(uj+1+uj－1)(5)22Crank－Nicolson格式的稳定性分析由热传导方程的Crank－Nicolson格式，可将原来方程写为下面的差分方程组:n+1nAnU=BnU+enn=0，1…N－1{(6)0U=φ－1以下讨论差分方程系数不依赖于时间层面，即Ai=A，Bi=B，故Ci=AB=C上面(6)式可以写为n+1nAU=BU+enn=0，1…N－1{(7)0U=φnT这时稳定性条件为:存在常数K使得C≤K0≤n≤k令λ1，λ2…λM－1为C的特征值，用ρ(C)表示λi的最大值，即C的谱半径，则有:定理1差分格式(7)稳定的

7、必要条件是，存在与k无关的常数c0，使得矩阵C=－1AB的谱半径满足ρ(C)≤1+c0k(8)nnn证明:因为对于所有的n＞0，有ρ(C)=ρ(C)≤C故若差分格式(7)稳定，必有nTρ(C)≤K0＜n≤(9)k容易算出(8)式和(9)式是等价的．事实上，若(8)式成立，则有Tnnc0Tρ(C)≤(1+c0k)≤(1+c0k)k≤e=K(T－k)TkklnK反之，若(9)式成立，特别取≤n≤，则ρ(C)≤KT－k=eT－kkk2lnKklnK2=1+k()+()+…＜1+c0kT－k2!T－klnKk0lnK这里取c0=eT－k00＜k＜k0，证毕．T－k0

8、定理2若A为正规矩阵，则A=ρ(A)．2［4］－1*