本科毕业设计--求解热传导方程的高精度隐式差分格式

1.2.2 隐式差分格式的稳定性

证:先用差分格式（1.2.1）为：

?1nnnun?unjj??(uj?1?2uj?uj?1) 其中????h2

nikjh利用稳定性的Fourier方法令，un，并将它代入上式就得到 j?vevn?1eikjh?vneikjh??(eikj?2?e?ikj)vn

消去公因子有

vn?1?[1??(eikj?2?e?ikj)]vn?vn?由此得到增长因子

1??(eikj1vn?1 ?ikj?2?e)G(?,k)?11?1??(eikj?2?e?ikj)(1?2?)??(eikj?e?ikj)

11??(1?2?)?2?coskh1?2?(1?coskh)11??1|1?2?(1?coskh)|1?4?sin2kh2kh??h?0,k?0,??0且0?sin2?12

khkh?4?sin2?0?1?4?sin2?1221???1kh1?4?sin22显然这个格式是相容的。它是无条件稳定，因为按照Lax定理可知，该格式收敛的。

|G(?,k)|?1.3 、 Crank-Nicolson格式

我们在在前面讨论的显格式和隐格式，即：

n?1un?ujj?1n?1n?1un?2u?uj?1jj?1?n?1unj?uj??h2?0 （1.3.1）

???nnunj?1?2uj?uj?1h2?0 （1.3.2）

用?乘（1.3.2），用(1??)乘（1.3.1），把其结果相加就得到一个差分格式

6

n?1unj?uj

nn?1n?1?1?un?un?unj?1?2uj?uj?1j?1?2ujj?1?????(1??)??0 （1.3.3） 22?hh????其中0???1，我们乘差分格式公式(1.3.3）为加权隐式格式。

从上述可以看到，当???1un?unjj1时的情况，此时我们把它单独写出 2???2hn?1n?1n?1nnn[(u?2u?u)?(u?2u?uj?1jj?1j?1jj?1)]?0 （1.3.4） 2此格式一般称作Crank-Nicolson格式。此外我们注意到，当??1时，公式

(1.3.3）为向后差分格式（隐式格式）；当??0时，公式(1.3.3）为向前差分格式（显式格式）。

1.3.1 Crank-Nicolson差分格式的截断误差

证：（用taylor展开）

n?1un?u(x,t),u?u(xj,tn?1)jjnjuun?1j?un?2?2un?3?3un??1?u(x,t??t)?u(xj,tn)??[]j?[2]j?[3]j?t2!?t3!?t?unh2?2unh3?3unh4?4un?u(x??x,t)?u(xj,tn)?h[]j?[2]j?[3]j?[4]j??2?x2!?x3!?x4!?x?unh2?2unh3?3unh4?4un?u(x??x,t)?u(xj,tn)?h[]j?[2]j?[3]j?[4]j??3?x2!?x3!?x4!?x(0??1?1)(0??2?1)(0??3?1)nj?1unj?1

?把上述代入差分格式中，得截断误差为：

?un??2un?2?3un??1??2un?1h2?4un?1k(x,t)?[]j?[2]j?[]j?{([2]j?[]j??4?t2?t6?t32?x24?x4h2?4un?1?2unh2?4unh2?4un?[4]j??5)?([2]j?[4]j??2?[4]j??3)}24?x?x24?x24?xnj?un??2un?1?2un??2un?2?3un??1?{[]j?([2]j?[2]j)}?{{[2]j?[3]j}?t2?x?x2?t6?th2?4un?1?4un?1?4un?4un??{([4]j??4?[4]j??5)?([4]j??2?[4]j??3)}}

24?x?x?x?x?o(?2?h2).(0??1,?2,?3,?4,?5?1)

22从上述可知，截断误差为kn(x,t)?o(??h)，它对空间方向为二阶截断误j差，而对时间方向为二阶截断误差，则此隐格式的精度为2。

7

1.3.2 Crank-Nicolson差分格式的稳定性

证：由差分格式公式（1.3.4）可以写成如下形式

?1111nnn?1n?1?1??un?(1???)u???u???u?(1???)u???unj?1jj?1j?1jj?12222其中???h2消去公因子有

11??eikj?(1???)???e?ikj2vn?2vn?1

11???eikj?(1???)???e?ikj22由此得到增长因子

111??eikj?(1???)???e?ikj(1???)???(eikj?e?ikj)22G(?,k)?2?111???eikj?(1???)???e?ikj(1???)???(eikj?e?ikj)222kh1?2??sin21???(1?coskh)2??1???(1?coskh)1?2??sin2kh2h?0,k?0,??0且1?2??sin2kh2?1?|G(?,k)|?kh1?2??sin221?2??sin2khkh?1?2??sin222

显然这个格式是相容的。它是无条件稳定，因为按照Lax定理可知，

Crank-Nicolson差分格式收敛的。

8

2． 高精度格式的构造

本文热传导方程对空间变量应用四阶紧致格式离散，对时间变量应用梯形方

24法，构造扩散方程的精度为O??h的绝对稳定的隐式差分格式

??2.1梯形方法[4]

求解常微分方程初值问题

?dy??f(x,y)?dx?y(x0)?y0

?

对方程

dy?f(x,y)两边从xn到xn?1积分，得 dxxn?1xny(xn?1)?y?xn???f(x,y(x))dx (2.1)

用左矩形公式计算上式右侧积分，即

?xn?1xnf(x,y(x))dx?hf(xn,y(xn))

并用yn作为的近似值y(xn)，得

yn?1?yn?hf?xn,yn? (2.2） 故欧拉法也称为矩形法。

欧拉法形式简单,但精度低，为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算(2.1)式右侧积分，即

?xn?1xnf(x,y(x))dx?h?f?xn,y?xn???f?xn?1,y?xn?1???2

并用yn作为y(xn)的近似值，得到梯形公式:

yn?1?yn?h?f?xn,yn??f?xn?1,yn?1??2 (2.3)

(2.3)式称为梯形公式,此方法称为梯形法。

9

2.2本文格式的构造

下面我们考虑非齐次边界条件的热传导方程

??u?2u???0,a?x?b,t?0?2?t?x??a?x?b?u(x,0)?f(x),?u(a,t)?g1(t)t?0?t?0??u(b,t)?g2(t)（2.2.1）的右边对x变量 四阶紧致格式离散

(2.2.1)

?x2?2u?2??ui?o(h4)12?xh2(1??x)12 (2.2.2)

把(2.2.2)代入(2.2.1)后得到下面的常微分方程组

dui(t)12?1222 ?(1??x)?xui，令A??2(1?1?x),B??xdt1212h

dui(t)?A?1Bui?A?1C(t) (2.2.3) dtTU(t)?[u1(t),u2(t),...,um?2(t),um?1(t)]T,U(0)?[g(x1),g(x2),...,g(xm?2),g(xm?1)].C(t)?[（-

1'?1'?u0(t)+(2)u0(t)），0，,0，（-uN(t)+(2)uN(t)）]' 12h12h

?5?6??1其中：A??12?????11256112112???21???1?2?B????，

?????5??6?11???? ??2???dU(t)?A?1BU(t)+A?1C(t),?应用梯形方法得到本文格式 ?dt??U(0)?U0un?1?(I??2A?1B)?1(I??1??A?1B)un?（I??A?1B)?1(A?1C(tn)?A?1C(tn?1)2222 （2.2.4）

10