一维热传导方程的高效解法：Matlab谱方法原理与实现

考虑一维热传导问题：

其中a为正常数，f(x,t)、φ(x)、α(t)和β(t)为已知函数。u(x,0)=φ(x)给出了函数u在t=0时刻的取值，即初始条件。u(0,t)=α(t)和u(1,t)=β(t)给出了函数u在边界上的取值，即边界条件。

首先，以空间步长h=1/N、时间步长τ=T/M分别将x轴上区间[0,1]、t轴上区间[0,T]分成N、M等分，得x_i=ih、0≤i≤N和t_k=k_τ、0≤k≤M。再根据x_i和t_k对计算区域做网格剖分，如图2-6所示。横线与竖线的交点即为结点，称在t=0、x=0或x=1上的结点为边界结点，其余结点为内结点。

图2-6 计算区域的网格剖分

然后在结点处考虑热传导问题：

根据泰勒公式，有：

将上面两式代入式（2-28），并将小量O(τ+h²)作为局部截断误差忽略掉，同时引入步

长比r=aτ/h²，并使用简写u^k_i=u(x_i,t_k)，得差分格式：

图2-7显示了差分格式（2-31）所涉及结点的相对位置。要计算函数u在(x_i,t_k+1)处的取值，就需要用到其在(x_i-1,t_k)、(x_i,t_k)和(x_i+1,t_k)处的取值。(www.xing528.com)

图2-7 差分格式（2-31）所涉及结点的相对位置

最后，将式（2-31）写成矩阵形式（2-32），注意最右端向量的首尾元素已包含了边界条件。这样，通过一次矩阵相乘和相加就可以由(u^k₁,u^k₂,…,u^k_N-2,u^k_N-1)^T得到(u^k+1₁,u^k+1₂,…,u^k+1_N-2,u^k+1_N-1)^T。

需要强调的是，当步长比r=aτ/h²≤1/2时，差分格式（2-31）是稳定的，反之则是不稳定的，计算结果也不可信。这一稳定性条件称为CFL条件，该条件是根据提出者Courant、Friedrichs、Lewy的名字命名的，它保证了误差和噪声在迭代的每一步中都不会放大，证明从略。下面看一个实例：

根据有限差分法的矩阵形式（2-32）计算热传导问题（2-33），代码如下：

程序2-3

程序输出如图2-8所示，左图为数值结果，右图为数值结果与解析解u=e^x+t之间的误差，程序使用的步长比r=1/2，所以计算过程是稳定的。

图2-8 数值结果（左图）和误差（右图）