摘要 1. 前言
2. Crank-Nicolson差分法
2.1） 差分法定义 2.2） 差分格式的建立
2.3） Crank-Nicolson差分格式（六点格式） 2.4） Crank-Nicolson差分格式的向量表示 2.5） Crank-Nicolson差分格式的稳定性 2.6） Crank-Nicolson差分格式的收敛性
3. 数值算例
3.1） 利用Crank-Nicolson方法求解数值算例
4. 总结 5. 参考文献
6. 致谢
1 3 4
4 4 7 9 11 14
17
17
20 21 22
苏州大学本科生毕业设计（论文）
摘要
本篇文章主要介绍了一维热传导方程的混合问题的第一边界值问题，并采用六点（Crank−Nicolson）格式对该问题进行讨论。主要涉及到对该格式的推导与证明，并对Crank−Nicolson格式的稳定性与收敛性进行探讨；然后针对具体的微分方程使用该方法进行运算。最终使用matlab得出答案与画出该微分方程的图像。
关键字：一维热传导方程，第一边界值，六点格式，稳定性，收敛性
1
苏州大学本科生毕业设计（论文）
Abstract
This article mainly introduces the first boundary value problem for the mixed problem of one-dimensional heat conduction equations, and using the Crank-Nicolson scheme to solve this problem. It mainly involves the derivation and proof of this
scheme, and the analysis for the stability and convergence of this method. Finally, we use Matlab to give some numerical experiments to confirm our theoretical analysis.
Keywords: One-dimensional heat conduction equation, first boundary value, six-point format, stability, convergence
2
苏州大学本科生毕业设计（论文）
1 前言
偏微分方程指的是包括未知函数的偏导数或者偏微分的方程。这类方程在诸多领域内被广泛应用，其中在数学和物理以及工程技术中最为广泛，所以我们也在习惯上将这类微分方程称之为数学物理方程。
偏微分方程本质上是研究一个方程或者方程组的解的存在性和解的唯一性（或自由度）的各种性质的方法。其中解的存在性是指存在某个解满足某些补充条件；而解的唯一性则是研究有多少个解能够满足这些特定的条件。通常，我们不仅要对解的性质与求解方法进行探究，而且还要利用偏微分方程对某些自然现象进行解释和预测。
所以说，偏微分方程的形成与发展一直都与物理学和一些自然科学不可分割。并在诸多方面相互促进彼此的发展与进步。例如在数学上，几何，代数，分析，拓扑等学科理论都深受偏微分方程的影响。而物理学上最为常见的分子对的不规则运动都可以依赖偏微分方程进行解答。
因为分子与分子之间的间距较大，所以分子的运动依靠的是无规则的热运动以及分子之间的碰撞，从而实现其内部的能量迁移，并进一步的实现宏观的热量传递。
在物理学上我们知道，热能是由分子的不规则运动产生的，在热能流动的过程中由两个基本的过程：传导和对流。气体分子之间的间距比较大，气体依靠分子的无规则热运动以及分子间的碰撞，在气体内部发生能量迁移，从而形成宏观上的热量传递。而通过偏微分方程，可以利用物体内部温度的分布不随时间的变化而变化这一特性，实现一维定态热传导。
故而在本文中，我们主要研究的是一维热传导的第一边界值问题。第一部分我们通过差分格式的建立，在网格上将微商替换成差商，从而实现将原问题离散化这一操作。第二部分将通过差分格式进一步推导出Crank−Nicolson格式（六点格式），并针对其稳定性，收敛性做进一步探究。第三部分则是通过对数值算例的计算，比较其真解与数值解的误差来实现对Crank−Nicolson格式的正确性与易用性方面的验证。
3
苏州大学本科生毕业设计（论文）
2 Crank-Nicolson差分法 2.1 差分法定义
我们都知道，在研究热传导现象或者气体扩散现象时，将得到一个二阶抛物型偏微分方程。对此类方程我们可以提出初值问题和混合问题。 1）初值问题：
在区域𝑄:*−∞<𝑥<+∞,𝑡≥0+内求函数𝑢(𝑥,𝑡)满足方程：
𝐿𝑢≡
𝜕𝑢𝜕𝑡
−𝜕𝑥2=0,−∞<𝑥<+∞,𝑡>0
𝜕2𝑢
(2.1.1) (2.1.2)
及初始条件：
𝑢|𝑡=0=𝜑(𝑥),−∞<𝑥<+∞
其中𝐿表示微分算子∂𝑡−𝜕𝑥2，𝜑(𝑥)是给定的初始函数。
2）混合问题：
对于方程(2.1.1)我们可以提出第一、第二、第三边值问题[1]；为简单起见，在本论文中我们仅仅讨论如下的第一边值问题。
第一边值问题：
在区域𝑅:*0≤𝑥≤1，0≤𝑡≤𝑇+内求函数𝑢(𝑥,𝑡)满足方程
𝐿𝑢≡
𝜕𝑢𝜕𝑡
∂
𝜕2
−𝜕𝑥2=0,0<𝑥<1,0<𝑡≤𝑇
𝜕2𝑢
(2.1.3) (2.1.4) (2.1.5)
及初始条件：
𝑢|𝑡=0=𝜑(𝑥),0≤𝑥≤1
边界条件：
𝑢(0,𝑡)=𝜇1(𝑡),{0≤𝑡≤𝑇, 𝑢(1,𝑡)=𝜇2(𝑡),
并且𝜑(0)=𝜇1(0),𝜑(1)=𝜇2(1),其中𝜇1(𝛽)，𝜇1(𝛽)是给定的函数。
2.2 差分格式的建立
用差分格式求解(2.1.3)的第一步是将连续问题离散化。假设ℎ=𝛥𝑥是自变量𝑥的改变量，𝜏=𝛥𝑡是自变量𝑡的改变量，如此由𝑥𝑘=𝑘ℎ,𝑡𝑗=(𝑘=0,±1,±2,⋯;𝑗=
0,±1,±2,……)两组平行线组成的长方形网格将会覆盖整个𝑥−𝑡平面，其中ℎ被称为沿𝑥方向的步长，𝜏被称为沿𝑡方向的步长，而网格线的交点则被称为网格的结点。
故而边值问题的网格是：
4
苏州大学本科生毕业设计（论文）
𝑇
𝑡𝑗=𝑗𝜏,𝑗=0,1,2,…,𝑚0;𝑚0=⌈⌉, {𝜏 𝑥𝑘=𝑘ℎ,𝑘=0,1,2,…,𝑁;𝑁ℎ=1.
而在𝑡=0,𝑥=0,𝑥=1上的结点被称为边界结点，其它所有属于*0<𝑥<1,0<𝑡≤𝑇+
内部的结点则被称为内部结点。
所以不难看出差分方法的本质原理就是利用网格的结点的差商来代替微商，从而求出微分方程解的近似值，故这一方法又被称为网格法。而为了方便起见，我们可以将网格点(𝑥𝑘,𝑡𝑗)简记为(𝑘,𝑗)。
接下来我们在点(𝑘,𝑗)处列出部分常用的公式：
(𝜕𝑥2)𝜕2𝑢
(𝑘,𝑗)
=
𝑢(𝑘+1,𝑗)−2𝑢(𝑘,𝑗)+𝑢(𝑘−1,𝑗)
ℎ2−12𝑢𝑥4(𝑥̃,𝑗), |𝑥̃−𝑥𝜅|≤ℎ
ℎ2
(4)
(2.2.1)
其中：
(𝜕𝑡)𝜕𝑢
(𝑘,𝑗)
=
𝑢(𝑘,𝑗+1)−𝑢(𝑘,𝑗)
𝜏
′′
−2𝑢𝑡2(𝑘,𝑡),
𝜏
(2.2.2)
其中：
|𝑡−𝑡𝑗|≤𝜏
(𝜕𝑡)𝜕𝑢
(𝑘,𝑗)
=
𝑢(𝑘,𝑗)−𝑢(𝑘,𝑗−1)
𝜏
𝜏′′̃−𝑢𝑡2(𝑘,𝑡), 2
(2.2.3)
其中：
̃−𝑡|≤𝜏 |𝑡𝑗−1
(𝜕𝑡)𝜕𝑢
(𝑘,𝑗)
=
𝑢(𝑘,𝑗+)−𝑢(𝑘,𝑗−)
𝜏
1
2123
−24𝑢𝑡3(𝑘,𝑡′),
𝜏2
()
(2.2.4)
其中：
|𝑡′−𝑡𝑗−1|≤𝜏/2
2公式(2.2.1)的右侧的第一项是𝑢(𝑥,𝑡)关于自变量𝑥的二阶中心差商，而公式(2.2.2)、(2.2.3)、(2.2.4)的右侧的第一项则分别是𝑢(𝑥,𝑡)关于自变量𝑡的向前差商、向后差商和中心差商。
接下来，我们来简单推导𝑢(𝑥,𝑡)一阶差商和二阶差商：将函数𝑢(𝑘,𝑗+1)以结点𝑢(𝑘,𝑗)为中心做关于𝑡的Lagrange−Taylor级数展开得
𝜕𝑢𝜏2𝜕2𝑢
𝑢(𝑘,𝑗+1)=𝑢(𝑘,𝑗)+𝜏∙()+∙(2),𝑡∈(𝑡𝑗,𝑡𝑗+1)
𝜕𝑡(𝑘,𝑗)2𝜕𝑡(𝑘,𝑡
̃)
故有：
𝜕𝑢𝑢(𝑘,𝑗+1)−𝑢(𝑘,𝑗)𝜏′′()=−∙𝑢𝑡2(𝑘,𝑡) 𝜕𝑡(𝑘,𝑗)𝜏2至此关于变量𝑡的一阶向前差商推导完毕，关于变量𝑡的一阶向后差商以及关于变
量𝑥的一阶向前向后差商均与此类似，故在此不做推导。
将函数𝑢(𝑥𝑘,𝑡𝑗+1)以结点(𝑥𝑘,𝑡𝑗)为中心关于𝑡运用Lagrange−Taylor级数展开
25
苏州大学本科生毕业设计（论文）
𝜏𝜕𝑢1𝜏2𝜕2𝑢
𝑢(𝑥𝑘,𝑡𝑗+1)=𝑢(𝑥𝑘,𝑡𝑗)+∙()+∙()∙(2)
2𝜕𝑡2!2𝜕𝑡()2𝑘,𝑗(𝑘,𝑗)1𝜏3𝜕3𝑢1𝜏4𝜕4𝑢
+∙()∙(3)+·()·(4),𝑡∈(𝑡𝑗,𝑡𝑗+1) 3!2𝜕𝑡(𝑘,𝑗)4!2𝜕𝑡(𝑘,𝑡
̃)
将函数𝑢(𝑥𝑘,𝑡𝑗−1)以结点(𝑥𝑘,𝑡𝑗)为中心关于𝑡运用Lagrange−Taylor级数展开
2𝜏𝜕𝑢1𝜏2𝜕2𝑢
𝑢(𝑥𝑘,𝑡𝑗−1)=𝑢(𝑥𝑘,𝑡𝑗)+(−)∙()+∙(−)∙(2)
2𝜕𝑡(𝑘,𝑗)2!2𝜕𝑡(𝑘,𝑗)21𝜏3𝜕3𝑢1𝜏4𝜕4𝑢
+∙(−)∙(3)+·()·(4),𝑡∈(𝑡𝑗,𝑡𝑗+1) 3!2𝜕𝑡(𝑘,𝑗)4!2𝜕𝑡(𝑘,𝑡
̃)
故：
𝜕𝑢𝜏3𝜕3𝑢∙𝜏=𝑢(𝑥𝑘,𝑡𝑗+1)−𝑢(𝑥𝑘,𝑡𝑗−1)−∙ 𝜕𝑡24𝜕𝑡322则：
𝜕𝑢
=𝜕𝑡𝑢(𝑥𝑘,𝑡𝑗+1)−𝑢(𝑥𝑘,𝑡𝑗−1)
22𝜏𝜏2𝜕3𝑢−∙ 24𝜕𝑡3将函数𝑢(𝑥𝑘,𝑡𝑗+1)以结点𝑢(𝑥𝑘,𝑡𝑗)为中心关于𝑡运用Lagrange−Taylor级数展开
𝜕𝑢12𝜕2𝑢
𝑢(𝑥𝑘,𝑡𝑗+1)=𝑢(𝑥𝑘,𝑡𝑗)+()·𝜏+∙𝜏∙(2)
𝜕𝑡(𝑘,𝑗)2!𝜕𝑡(𝑘,𝑗)13𝜕3𝑢14𝜕4𝑢
+∙𝜏∙(3)+·𝜏·(4),𝑡∈(𝑡𝑗,𝑡𝑗+1) 3!𝜕𝑡(𝑘,𝑗)4!𝜕𝑡(𝑘,𝑡
̃)
将函数𝑢(𝑥𝑘,𝑡𝑗−1)以结点𝑢(𝑥𝑘,𝑡𝑗)为中心关于𝑡运用Lagrange−Taylor级数展开
𝜕𝑢−𝜏2𝜕2𝑢
𝑢(𝑥𝑘,𝑡𝑗−1)=𝑢(𝑥𝑘,𝑡𝑗)+()·(−𝜏)+()∙(2)
𝜕𝑡(𝑘,𝑗)2!𝜕𝑡(𝑘,𝑗)−𝜏3𝜕3𝑢1𝜕4𝑢4
+()∙(3)+·(−𝜏)·(4),𝑡∈(𝑡𝑗,𝑡𝑗+1)
3!𝜕𝑡(𝑘,𝑗)4!𝜕𝑡(𝑘,𝑡
̃)
故：
𝜕2𝑢𝜏4(4)
2
(2)·𝜏=𝑢(𝑥𝑘,𝑡𝑗+1)+𝑢(𝑥𝑘,𝑡𝑗−1)−2𝑢(𝑥𝑘,𝑡𝑗)−𝑢(𝑥𝑘，𝑡𝑗) 𝜕𝑡(𝑘,𝑗)12因此：
6
苏州大学本科生毕业设计（论文）
𝜕2𝑢𝑢(𝑥𝑘,𝑡𝑗+1)−2𝑢(𝑥𝑘,𝑡𝑗)+𝑢(𝑥𝑘,𝑡𝑗−1)𝜏2(4)
=−𝑢𝑡4(𝑥𝑘,𝑡𝑗), 𝜕𝑡2𝜏212
至此，关于变量𝑡的二阶中心差商推导完毕，关于变量𝑥的二阶中心差商与此类似，
故在此不做推导。
2.3 Crank-Nicolson差分格式（六点格式）：
建立该差分格式的主要思想是在某些离散点上用差商代替微分方程中的
𝜕2𝑢𝜕𝑢𝜕𝑥2𝜕𝑡
,项利用不同的差商格式跌取相应的离散点，从而得到不同的差分方法；下
面我们将着重介绍Crank−Nicolson六点格式。
在点(𝑘,𝑗+2)处利用𝑢(𝑥,𝑡)关于𝑡的中心差商得到：
(𝜕𝑡)(𝑘,𝑗+1)=
21
𝜕𝑢𝑢(𝑘,𝑗+1)−𝑢(𝑘,𝑗)
𝜏(3)
+24𝑢𝑡3(𝑘,𝑡)
𝜏2
(2.3.1)
其中：
|𝑡−𝑡𝑗+1|≤𝜏⁄2
2利用Taylor公式：
𝜕2𝑢𝜕2𝑢𝜏𝜕3𝑢𝜏2𝜕4𝑢(2)=(2)+()+(), 𝜕𝑥(𝑘,𝑗+1)𝜕𝑥(𝑘,𝑗+1)2𝜕𝑥2𝜕𝑡(𝑘,𝑗+1)8𝜕𝑥2𝜕𝑡2(𝑘,𝜂)
221
𝜕2𝑢𝜕2𝑢𝜏𝜕3𝑢𝜏2𝜕4𝑢(2)=(2)−()+(), 𝜕𝑥(𝑘,𝑗)𝜕𝑥(𝑘,𝑗+1)2𝜕𝑥2𝜕𝑡(𝑘,𝑗+1)8𝜕𝑥2𝜕𝑡2(𝑘,𝜂)
222
其中：
𝜏𝜏
𝜂1=𝑡𝑗+1+𝜃1,𝜂2=𝑡𝑗−𝜃2
2220<𝜃1<1,0<𝜃2<1
故有：
𝜕2𝑢1𝜕2𝑢1𝜕2𝑢𝜏2𝜕4𝑢𝜏2𝜕4𝑢
(2)=()+()−()−() 𝜕𝑥(𝑘,𝑗+1)2𝜕𝑥2(𝑘,𝑗+1)2𝜕𝑥2(𝑘,𝑗)16𝜕𝑥2𝜕𝑡2(𝑘,𝜂)16𝜕𝑥2𝜕𝑡2(𝑘,𝜂)
21
2
=
𝑢(𝑘+1,𝑗+1)−2𝑢(𝑘,𝑗+1)+𝑢(𝑘−1,𝑗+1)+𝑢(𝑘+1,𝑗)−2𝑢(𝑘,𝑗)+𝑢(𝑘−1,𝑗)
2ℎ2ℎ2(4)ℎ2(4)𝜏2(4)𝜏2(4)
̃−𝑢𝑥4(𝑥,̃𝑗+1)−𝑢𝑥4(𝑥̃,𝑗)−𝑢𝑥2𝑡2(𝑘,𝜂1)−𝑢𝑥2𝑡2(𝑘,𝜂2),
24241616 (2.3.2)
7
苏州大学本科生毕业设计（论文）
由(2.3.1), (2.3.2)得到：
(𝐿𝑢)
1(𝑘,𝑗+)
2𝜕𝑢𝜕2𝑢1(3)(3)=(−2)=𝐿ℎ,𝜏𝑢(𝑘,𝑗+)+𝑅ℎ,𝜏,
𝜕𝑡𝜕𝑥(𝑘,𝑗+1)22其中：
1𝑢(𝑘,𝑗+1)−𝑢(𝑘,𝑗)(3)
𝐿ℎ,𝜏𝑢(𝑘,𝑗+)=−
2𝜏1𝑢(𝑘+1,𝑗+1)−2𝑢(𝑘,𝑗+1)+𝑢(𝑘−1,𝑗+1)
· 2ℎ21𝑢(𝑘+1,𝑗)−2𝑢(𝑘,𝑗)+𝑢(𝑘−1,𝑗)−· 2ℎ2
(3)𝑅ℎ,𝜏
𝜏2(3)ℎ2(4)ℎ2(4)
̃=𝑢3(κ,𝑡)+𝑢𝑥4(𝑥,̃𝑗+1)+𝑢𝑥4(𝑥̃,𝑗)
24𝑡2424𝜏2(4)𝜏2(4)
+𝑢𝑥2𝑡2(𝑘,𝜂1)+𝑢𝑥2𝑡2(𝑘,𝜂2) 1616=Ο(ℎ2+𝜏2)
因为：
(𝐿𝑢)
1(𝑘,𝑗+)2
=0,
(2.3.3)
所以：
1(3)(3)
𝐿ℎ,𝜏𝑢(𝑘,𝑗+)+𝑅ℎ,𝜏=0.
2(3)
略去𝑅ℎ,𝜏不计，则得差分方程：
3
𝐿ℎ,𝜏𝑢𝑘,𝑗+1=
2()
𝑢𝑘,𝑗+1−𝑢𝑘,𝑗
−
𝜏1𝑢𝑘+1,𝑗+1−2𝑢𝑘,𝑗+1+𝑢𝑘−1,𝑗+1+𝑢𝑘+1,𝑗−2𝑢𝑘,𝑗+𝑢𝑘−1,𝑗· 2ℎ2=0
以上方程还可以改写为：
令ℎ2=𝑟，则：
𝑟𝑟
−𝑢𝑘+1,𝑗+1+(1+𝑟)𝑢𝑘,𝑗+1−𝑢𝑘−1,𝑗+1 22𝑟𝑟
=(1−𝑟)𝑢𝑘,𝑗+𝑢𝑘+1,𝑗+𝑢𝑘−1,𝑗.
22
(2.3.4)
𝜏
8
苏州大学本科生毕业设计（论文）
从上述方程可知，在点(𝑘,𝑗+2)处列方程时我们要用到第𝑗+1层上的三个点
(𝑘,𝑗+1),(𝑘+1,𝑗+1),(𝑘−1,𝑗+1)及第𝑗层上的三个点(𝑘,𝑗),(𝑘+1,𝑗),(𝑘−1,𝑗)。因此
1
这一格式被称为六点格式，可以利用下面的图形来表示：
(𝑘−1,𝑗+1) (𝑘,𝑗+1) (𝑘+1,𝑗+1) j+1 1(𝑘,𝑗+) 2j (𝑘−1,𝑗) (𝑘,𝑗) (𝑘+1,𝑗)
这一差分格是隐式差分格式的一类，且由式(2.3.3)可知，当ℎ→0,𝜏→0时，差分方程(2.3.4)将会逼近微分方程(2.1.3)，故而其截断误差为𝛰(ℎ2+𝜏2)。
2.4 Crank-Nicolson差分格式的向量表示
下面介绍Crank-Nicolson差分格式的向量表示。假设长方形网格是：
{
𝑥𝑘=𝑘ℎ,𝑘=0,1,2,…,𝑁;𝑁ℎ=1. 𝑡𝑗=𝑗𝜏,𝑗=0,1,2,…,𝑚0;𝑚0=⌈𝜏⌉,
𝑇
为了把所述的边值问题化为近似的差分方程问题，我们必须对初始条件和边界条件建立相应的差分相似。由于我们考虑的区域𝑅:*0≤𝑥≤1,0≤𝑡≤𝑇+是一个矩形，该网格的结点全部在边𝑥=0,𝑥=1和𝑡=0上，所以对于第一边值问题可以根据给定的条件(2.1.4), (2.1.5)直接得到该差分方程所相应的初始条件：
𝑢𝑘,0=𝜑(𝑘ℎ),𝑘=0,1,2,…,𝑁 (2.4.1) 和边界条件：
𝑢0,𝑗=𝜇1(𝑗𝜏),𝑢𝑁,𝑗=𝜇2(𝑗𝜏),𝑗=0,1,…,𝑚0 (2.4.2) 故，Crank−Nicolson差分格式相应的差分方程问题为：
2(1+𝑟)𝑢𝑘,𝑗+1−𝑟𝑢𝑘+1,𝑗+1−𝑟𝑢𝑘−1,𝑗+1=2(1−𝑟)𝑢𝑘,𝑗+𝑟𝑢𝑘+1,𝑗+𝑟𝑢𝑘−1,𝑗 (2.4.3) 𝑘=1,2,…,𝑁−1;𝑗=0,1,2,…,𝑚0−1
()𝑢=𝜑𝑘ℎ,𝑘=0,1,…,𝑁, 𝑘,0(2.4.4)
𝑢=𝜇(𝑗𝜏),𝑢=𝜇(𝑗𝜏),𝑗=0,1,2,…,𝑚. (2.4.5)
1𝑁,𝑗20{0,𝑗
根据初始条件(2.4.4)和边界条件(2.4.5)，通过方程组(2.4.3)可逐层求得𝑢𝑘,𝑗 令𝑗=0,则有：
9
苏州大学本科生毕业设计（论文）
2(1+𝑟)𝑢𝑘,1−𝑟𝑢𝑘+1,1−𝑟𝑢𝑘−1,1=2(1−𝑟)𝑢𝑘,0+𝑟𝑢𝑘+1,0+𝑟𝑢𝑘−1,0,
𝑘=1,2,…,𝑁−1,
𝑢𝑘,0=𝜑(𝑘ℎ),𝑘=0,1,…,𝑁, {𝑢0,0=𝜇1(0),𝑢𝑁,0=𝜇2(0),
可以看出方程的左端是一个以𝑢1,1,𝑢2,1,…,𝑢𝑁−1,1为未知数的线性代数方程组，其系数矩阵为：
2(1+𝑟)−𝑟00…00
−𝑟2(1+𝑟)−𝑟0…00
0−𝑟2(1+𝑟)−𝑟…00 ………………… 0000…−𝑟2(1+𝑟)]
方程的右端的是一个以𝑢1,0,𝑢2,0,…,𝑢𝑁−1,0为系数的线性代数方程组，其系数矩阵为：
2(1−𝑟)𝑟00…00
𝑟2(1−𝑟)𝑟0…00 𝐴′= 0𝑟2(1−𝑟)𝑟…00
……………… …
[0000…𝑟2(1−𝑟)]
为方便起见，可将方程组(2.4.3)写成向量形式。令：
𝑟𝑢0,𝑗𝑢0,𝑗
𝑢1,𝑗𝜑(ℎ) 0 0 𝑢2,𝑗(2ℎ) 𝜑 𝒖𝒋=[⋮],𝒇𝒋= ⋮ =𝑟 ⋮ ,𝝋=
⋮ 0 0 𝑢𝑁−1,𝑗[𝜑((𝑁−1)ℎ)][𝑟𝑢𝑁,𝑗][𝑢𝑁,𝑗]
𝐵′=
[
矩阵𝐴为：
1−2𝑟 𝑟 𝐴= 0
…[0
𝑟1−2𝑟𝑟…0
0𝑟1−2𝑟…0
00𝑟…000−𝑟…0
…00…00
…00 ……… …𝑟1−2𝑟]…0…0…0………−𝑟
00
0 … 1+2𝑟]
矩阵𝐵为：
1+2𝑟−𝑟0 −𝑟1+2𝑟−𝑟 𝐵= 0−𝑟1+2𝑟
…… …
[000
则方程组(2.4.3),(2.4.3),(2.4.4)可以简记为：
(𝐼+𝐵)𝒖𝒋+𝟏=(𝐼+𝐴)𝒖𝒋+𝒇𝒋+𝟏+𝒇𝒋,𝑗=0,1,2,…,𝑚0−1{ 𝒖𝟎=𝝋,
(2.4.5)
至此我们完成了热传导方程的第一边值问题相应差分格式的建立，通过这一格式在接下来的问题中我们就能够将微分方程问题转化为线性代数方程问题。由于这一过程的特殊操作，所以我们又可以将其称之为微分方程的离散化过程。在接下来我们所要求解的方程中，由于所得的方程组的系数矩阵都是三对角矩阵，故我们可以使用追赶法对其求解。
10
苏州大学本科生毕业设计（论文）
2.5 Crank-Nicolson差分格式的稳定性
我们在建立差分格式的时候可以看到，计算数值𝑢𝑘,𝑗我们是在所考虑的网格区域内是分层进行的。所以我们如果在第𝑗层如果产生了误差，那么势必会影响到第𝑗+1层以及以后各个层的数值结果。但是又因为我们在解差分方程组时引入误差也是不可避免的，例如诸多的问题都含有舍入误差和初始误差，故而我们需要研究这些误差在计算过程中对结果所产生的具体影响，这也就是所谓的稳定性问题。
如果我们计算开始时引入的误差在以后逐层计算过程中的对结果所造成的影响逐渐消失或保持有界，那么我们可以称该差分格式是稳定的；否则我们就称该差分格式是不稳定的。通常情况下不稳定的差分格式即便是收敛也是没有实用价值的，因为此时误差积累会随步数的增多而逐渐扩大，进而在结果中将会彻底淹没真解。
而为了我们讨论方便起见，我们假定边界值的计算是精确的，这样就可以在初始层引入误差𝜀0（同样我们也可以假定是在计算过程中的某一层上产生了误差，但这样讨论的结果也是一样的），而在以后各层的计算未引入其他误差。
设
𝜺𝟎=[𝜀1,0,𝜀2,0,…,𝜀𝑁−1,0],𝜺𝒋=[𝜀1,𝑗,𝜀2,𝑗,…,𝜀𝑁−1,𝑗].
𝑇
𝑇
则前面我们所讨论的差分格式可以写成以下向量形式：
𝒖=𝐻𝒖𝒋−𝟏+𝒇𝒋,{𝒋 𝒖𝟎=𝝋,
其中
𝒖𝒋=[𝑢1,𝑗,𝑢2,𝑗,…,𝑢𝑁−1,𝑗]
𝐻是(𝑁−1)×(𝑁−1)矩阵。
̃𝒋应满足 差分方程的近似解𝒖
̃𝒋=𝐻𝒖̃𝒋−𝟏+𝒇𝒋,𝒖{ ̃𝟎=𝝋+𝜺𝟎. 𝒖
̃𝒋−𝒖𝒋所满足的方程是 则误差𝜺𝒋=𝒖
{
𝜺𝒋=𝐻𝜺𝒋−𝟏,𝒋=𝟏,𝟐,…,𝒎𝟎,𝜺𝟎为初始误差.
𝜺𝒋=𝐻𝑗𝜺𝟎
𝑻
由此方程可知 因此
‖𝜺𝒋‖≤‖𝐻‖𝑗‖𝜺𝟎‖,𝑗=1,2,…,𝑚0
由此可见，无论𝑗如何，若‖𝐻‖𝑗=𝐶(𝐶为与ℎ,𝜏无关的常数)，则有 ‖𝜺𝒋‖≤𝐶‖𝜺𝟎‖
这样就可以说明误差𝜺𝟎对以后各层的影响是有限的，因为‖𝐻‖小于1，所以在以后各层上的误差必定是在逐渐减小的；所以在这样的两种情况下，我们可以认为差分格式是稳定的，因此一般可以如下定义：
11
苏州大学本科生毕业设计（论文）
若𝜺𝒋在一定范数下满足不等式
‖𝜺𝒋‖≤𝐶‖𝜺𝟎‖,（𝑗≥1）
我们就称该差分格式是稳定的，其中𝐶为与ℎ,𝜏无关的常数。
定理一：对于初始误差（或是某一层产生的误差）𝜺𝟎，若存在小于1的常数𝑪满足不等式：‖𝜺𝒋‖≤𝑪‖𝜺𝟎‖,（𝒋≥𝟏），则称该差分格式为稳定差分格式。 接下来我们开始证明：
易知由隐式差分格式的方程组可以推断出其误差方程，满足：
𝐵𝒖=𝒖𝒋−𝟏+𝒇𝒋,𝑗=1,2,…,𝑚0{𝒋 𝒖𝟎=𝝋 ̃𝒋−𝒖𝒋满足的方程是： 所以误差𝜺𝒋=𝒖
𝐵𝜺𝒋=𝜺𝒋−𝟏,
{ 𝜺𝟎为初始误差
(2.5.1)
由上述方程我们可知六点格式的误差满足下列方程：
{
(𝐼+𝐵)𝜺𝒋+𝟏=(𝐼+𝐴)𝜺𝒋,j=0,1,…,𝑚0−1𝜺𝟎为初始误差
1−2𝑟𝑟0 𝑟1−2𝑟𝑟 𝐴= 0𝑟1−2𝑟
………[0001+2𝑟−𝑟0 −𝑟1+2𝑟−𝑟 𝐵= 0−𝑟1+2𝑟
…… …
[000
(2.5.2)
其中𝐼是单位矩阵，
0
0𝑟…000−𝑟…0
…00…00
…00 ……… …𝑟1−2𝑟]…00…00
…00 ……… …−𝑟1+2𝑟]
易知矩阵𝐴的特征多项式为：
1−2𝑟−𝜆
𝑟
𝐷𝑁−1(𝜆)=|
⋯0
𝑟
1−2𝑟−𝜆
⋯⋯
0⋯0𝑟⋯0
|
⋯⋯⋯0𝑟1−2𝑟−𝜆
按第一行展开，则有： 这是𝐷𝑚(𝜆)的递推公式
𝐷𝑁−1(𝜆)=(1−2𝑟−𝜆)𝐷𝑁−2(𝜆)−𝑟2𝐷𝑁−3(𝜆) 𝐷𝑚−1(𝜆)=(1−2𝑟−𝜆)𝐷𝑚−2(𝜆)−𝑟2𝐷𝑚−3(𝜆)
𝑚=3,4,5,…,𝑁
其中
𝐷0(𝜆)=1,𝐷1(𝜆)=1−2𝑟−𝜆
由以上的公式可以依次求得𝐷2(𝜆),𝐷3(𝜆),…,𝐷𝑁−1(𝜆)。另外，我们也可以将上述递推公式看成是𝐷𝑚(𝜆)在𝜆暂时固定的情况下的关于𝑚的二阶线性差分方程，通过求解差分方程我们可以得到𝐷𝑁−1(𝜆)的表达式。
12
苏州大学本科生毕业设计（论文）
若令𝐷𝑚(𝜆)=𝜇𝑚,那么特征方程为：
𝜇2−(1−2𝑟−𝜆)𝜇+𝑟2=0
此特征方程的两个根依次为：
11
𝜇1=(1−2𝑟−𝜆)+√(1−2𝑟−𝜆)2−4𝑟2,
2211
𝜇2=(1−2𝑟−𝜆)−√(1−2𝑟−𝜆)2−4𝑟2.
22
所以
𝑚𝑚
𝐷𝑚(𝜆)=𝐶1𝜇1+𝐶2𝜇2
设𝜌=2(1−2𝑟−𝜆),则有：
𝜇1=𝜌+√𝜌2−𝑟2,𝜇2=𝜌−√𝜌2−𝑟2 1
由
𝐷0(𝜆)=𝐶1+𝐶2=1,
𝐷1(𝜆)=𝐶1(𝜌+√𝜌2−𝑟2)+𝐶2(𝜌−√𝜌2−𝑟2)=1−2𝑟−𝜆=2𝜌
解得：
𝐶1=
𝜌+√𝜌2−𝑟22√𝜌2−𝑟2,𝐶2=
𝜌−√𝜌2−𝑟22√𝜌2−𝑟2 所以
𝐷𝑚(𝜆)=
12√𝜌2−𝑟2{(𝜌+√𝜌2−𝑚+12𝑟)
−(𝜌−√𝜌2−𝑚+1
2𝑟)}
现在求方程𝐷𝑁−1(𝜆)=0的根，即求矩阵𝐴的特征值
令
1
(1−2𝑟−𝜆)=𝜌=𝑟cos𝜃 2则
√𝜌2−𝑟2=√𝑟2(𝑐𝑜𝑠2𝜃−1)=𝑖𝑟sin𝜃
所以
1sin𝑁𝜃
{𝑟𝑁𝑒𝑖𝑁𝜃−𝑟𝑁𝑒−𝑖𝑁𝜃}=𝑟𝑁−1
𝑖2𝑟sin𝜃sin𝜃故𝐷𝑁−1(𝜆)=0的根，即sin𝑁𝜃=0的根。因此
𝑁𝜃=(𝑁−𝑘)𝜋,𝑘=1,2,…,𝑁−1
(𝑁−𝑘)𝜋θ=,𝑘=1,2,…,𝑁−1
𝑁𝐷𝑁−1(𝜆)=
由
1
(1−2𝑟−𝜆)=𝑟cos𝜃 2可得：
𝜆𝑘=1−2𝑟(1+cos𝜃)=1−4𝑟𝑐𝑜𝑠2
(𝑁−𝑘)𝑘𝜋
𝜋=1−4𝑟𝑠𝑖𝑛2()， 2𝑁2𝑁13
苏州大学本科生毕业设计（论文）
𝑘=1,2,…,𝑁−1
此为矩阵𝐴的特征值。
又因为矩阵𝐵和矩阵𝐴之间的关系为：
𝐵=𝐼−(𝐴−𝐼)=2𝐼−𝐴
因此矩阵𝐵的特征值为：
𝑘𝜋𝑘𝜋
𝜇𝑘=2−(1−4𝑟𝑠𝑖𝑛2())=1+4𝑟𝑠𝑖𝑛2(),
2𝑁2𝑁𝑘=1,2,…,𝑁−1
且𝐼+𝐵的特征值是𝐵的特征值加1，即
𝑘𝜋
2+4𝑟𝑠𝑖𝑛(),𝑘=1,2,…,𝑁−1
2𝑁2
并且特征向量与𝐵相同；
𝐼+𝐴的特征值是𝐴的特征值加1，即
𝑘𝜋
2−4𝑟𝑠𝑖𝑛2(),𝑘=1,2,…,𝑁−1
2𝑁并且特征向量与𝐴相同。
因为𝐼+𝐵的特征值均不为零，所以𝐼+𝐵的逆矩阵存在，误差方程可以改写为：
{
𝜺𝟎为初始误差
特征值是𝐼+𝐴的特征值与𝐼+𝐵的特征值的商，即
𝑘𝜋
2−4𝑟𝑠𝑖𝑛2(2𝑁)𝑘𝜋2+4𝑟𝑠𝑖𝑛2()2𝑁,𝑘=1,2,…,𝑁−1
𝜺𝒋+𝟏=(𝐼+𝐵)−𝟏(𝐼+𝐴)𝜺𝒋,j=0,1,…,𝑚0−1
故其绝对值无论𝑟>0如何都小于1，因此Crank−Nicolosn格式（六点格式）也是无条件稳定的。
2.6 Crank-Nicolson差分格式的收敛性
接下来我们将要讨论当ℎ→0,𝜏→0时，差分方程的精确解𝑢𝑘,𝑗是否会趋近于微分方程解𝑢(𝑘,𝑗)的问题。
(𝟑)𝟒)(𝟒)
定理二：若第一边值问题的解在𝑹内有连续偏导数𝒖(,𝒖,𝒖，不论𝒓>𝟎如何，𝟒𝟐𝟐𝒙𝒙𝒕𝒙𝟑
六点格式是一致收敛的。
设微分方程的解和差分方程精确解的差为：𝜌𝑘,𝑗=𝑢(𝑘,𝑗)−𝑢𝑘,𝑗,𝑟=𝜏⁄ℎ2=
(1)(1)
常数,𝑅ℎ,𝜏=𝛰(𝜏+ℎ2),其中𝑒𝑘,𝑗=−𝑅ℎ,𝜏。
14
苏州大学本科生毕业设计（论文）
设
𝜌1,𝑗𝑒1,𝑗𝜌1,𝑗𝑒1,𝑗
𝝆𝒊=[⋮],𝒆𝒋=[⋮]
𝜌𝑁−1,𝑗𝑒𝑁−1,𝑗
故六点格式的误差所满足的方程为：
(𝐼+𝐵)𝝆𝒋+𝟏=(𝐼+𝐴)𝝆𝒋+𝜏𝒆𝒋,
{
𝝆𝟎=0
其中𝒆𝒋的分量𝑒𝑘,𝑗=𝛰(𝜏2+ℎ2)。 前面证明(𝐼+𝐵)−1存在，故有：
𝝆𝒋+𝟏=(𝐼+𝐵)−1(𝐼+𝐴)𝝆𝒋+𝜏(𝐼+𝐵)−1𝒆𝒋,𝑗=0,1,…,𝑚0−1
{
𝝆𝟎=0
因为(𝐼+𝐵)−1的特征值是：
𝑘𝜋
(2+4𝑟𝑠𝑖𝑛2())2𝑁−1
,𝑘=1,2,…,𝑁−1
所以最大特征值为：
𝜋
(2+4𝑟𝑠𝑖𝑛2())2𝑁−1
𝜋ℎ
=(2+4𝑟𝑠𝑖𝑛2())
2−1
而
lim
𝜋𝑠𝑖𝑛2(2ℎ)
ℎ2𝜋2=() 2−1
ℎ→0
所以当ℎ→0,𝜏→0时
𝜋ℎ
(2+4𝑟𝑠𝑖𝑛2())
2≤1
故：
又因(𝐼+𝐵)−1(𝐼+𝐴)的特征值为：
‖(𝐼+𝐵)−1‖2≤1
𝑘𝜋ℎ
2−4𝑟𝑠𝑖𝑛2(2)
𝑘𝜋ℎ22+4𝑟𝑠𝑖𝑛(2)
当ℎ→0,𝜏→0时以下恒等式恒成立
𝑘𝜋ℎ
2−4𝑟𝑠𝑖𝑛2(2)||≤1
𝑘𝜋ℎ2+4𝑟𝑠𝑖𝑛2(2)
所以
‖(𝐼+𝐵)−1(𝐼+𝐴)‖2≤1
故有：
‖𝝆𝒋‖2≤𝑗𝜏
1·𝑁2·𝛰(𝜏+ℎ
2
2)
≤
1𝑇𝑁2·𝛰(𝜏2+ℎ2)
15
苏州大学本科生毕业设计（论文）
当ℎ2=𝑟=常数时，由于𝑁ℎ=1，所以
|𝜌𝑘,𝑗|≤
7 𝛰(ℎ2𝜏
+
3
ℎ2),(对一切𝑘,𝑗成立)
当ℎ=𝑠=常数时，则有
|𝜌𝑘,𝑗|≤
3
𝛰(ℎ2),(对一切𝑘,𝑗成立)
𝜏
(4)(4)(3)
如此证得，若第一边值问题的解在𝑅内有连续偏导数𝑢𝑥4,𝑢𝑥2𝑡2,𝑢𝑥3，不论𝑟>0如
何，六点格式是一致收敛的。
16
苏州大学本科生毕业设计（论文）
3数值算例
3.1 利用Crank-Nicolson方法求解数值算例
已知，微分方程为：
{𝑢(𝑥,0)=sin(𝜋𝑥),0≤𝑥≤1 𝑢(0,𝑡)=0,𝑢(1,𝑡)=0,0≤𝑡≤1
𝜕𝑢𝜕𝑡
=
𝜕2𝑢
,𝑥𝜖(0,1),𝑡𝜖(0,1] 𝜕𝑥2 (3.1.1)
其真解为：𝑢(𝑥,𝑡)=𝑒−𝜋𝑡𝑠𝑖𝑛(𝜋𝑥)
2
接下来进行探究当ℎ=0.01,𝑟=1与𝑟=1/2时数值解情况，并对𝒉取不同值
时对真解与数值解的误差进行探究。
当ℎ=0.01,𝑟=1时 调用Real函数
Real(0.01,1);
得出真解图像：
调用PDEParabolicCN函数：
PDEParabolicCN(1,1,inline(′sin(pi∗x)′),inline(′0′),inline(′0′),0.01,1)
17
苏州大学本科生毕业设计（论文）
得出数值解图像：
当ℎ=0.01,𝑟=1/2时
调用PDEParabolicCN函数：
PDEParabolicCN(1,1,inline(′sin(pi∗x)′),inline(′0′),inline(′0′),0.01,1/2) 得出数值解图像：
18
苏州大学本科生毕业设计（论文）
接下来我们来探究在ℎ取不同值时真解与误差解之间的误差图像： 其中ℎ的取值为：
11111,,,, 2223242526
调用Deviation中代码我们得到如下图像：
由此可见，在𝒉不断变小的情况下，真解与数值解的误差逐渐缩小。这与我们之前推断的𝐂𝐫𝐚𝐧𝐤−𝐍𝐢𝐜𝐨𝐥𝐬𝐨𝐧差分法的收敛性完全一致。因此在日常生活计算中，我们可以使用数值解来替代真解进行一些生产学习操作。
19
苏州大学本科生毕业设计（论文）
4总结
差分法是我们在学习生活上经常涉及到的数学方法，差分思想也是我们经常触及的思维方式。因为定解问题在许多情况下时不具有解析解的，亦或者说其解析解是不易计算的，这对学习与生产将会带来诸多的问题。所以在这样的情况下我们要采用可行的数值解法。
有限差分方法就是生活学习中常见易用的一种数值解法，它的核心基本思想是先对问题的定义域进行网格上的剖分，随后在网格点上，采用适当的数值微分公式实现利用差商对定解问题中的微商进行替换，从而实现离散化原问题，进而求出数值解。
此外，对差分格式的解的研究上，其存在性与唯一性、求解的方法与解法的数值稳定性、真解与数值解的误差分析等问题都需要我们进行耐心细致的研究与分析。
故而在对Crank−Nicolosn差分法研究时我们对其稳定性与收敛性进行进一步的探究，并对其数值解与真解的误差进行考量；以确保该方法可以获得正确易用的结果。
并且因为有限差分方法具有通用性强，构建简单，使用灵活等特点，并容易在计算机上进行实现。所以我们在计算过程中借用matlab数学工具对自己的计算进行了快速有效的验证。
20
苏州大学本科生毕业设计（论文）
5参考文献
[1] 《偏微分方程数值解法·南京大学数学系计算数学专业》--李荣华
21
苏州大学本科生毕业设计（论文）
6致谢
在本篇论文撰写的过程中，深受卢培培老师的帮助。无论是从最开始的选题，主体的确定，以及后期的文章编写，错误的查找、发现、修改，卢培培老师一直不辞辛苦。再次我要表达我深深的谢意。卢老师的治学之道是在我即将离开培育我四年土地之前又一意义深远的榜样。
并且非常感谢数学院各位老师的帮助。在我在苏大求学的这几年，对我的问题一直是知无不答，并且每一次回答都细致入微。老师们对待知识认真严谨的态度我定会牢记于心，以后无论走到哪里，经历了什么；这一份严谨治学，谦逊待人的人生态度永远不敢忘记。身为苏大的学生，数科院的孩子，学生在以后的学习中一定不断完善自己，不断提高对知识的严谨态度。
同时因为学生能力有限，所以在论文编写的过程依然存在诸多问题，在今后学习的过程中我一定更加努力，不断的完善。不辜负各位老师与苏大的培育之恩。
22
因篇幅问题不能全部显示，请点此查看更多更全内容