基于径向基函数的线方法求解热传导方程

The Method of Line for Heat Equation Based on Radial Basis Function

摘要: 文章借助线方法求解热传导方程的初边值问题，该方法是基于径向基函数的真正无网格方法，数值实验的结果表明了该方法的可行性和有效性。

Abstract: In this paper, initial and boundary value problem of heat equation is solved by the method of line which is a truly meshless method based on radial basis function. Numerical result shows that the approach is feasible and effective.

1. 引言

热传导方程是一类重要的拋物型偏微分方程。其一维形式如下

$\frac{\partial u\left(x,t\right)}{\partial t}=\alpha \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f(x)$

其中t是时间变量，x是空间变量， $u\left(x,t\right)$ 表示t时刻在位置x处的温度， $\alpha$ 为常数，它取决于材料的热传导率、密度和热容， $f\left(x\right)$ 为给定的函数，表示热源。作为模型方程，除了描述一个区域内温度如何随时间变化以外，热传导方程常也被应用于金融、自然环境、工程设备及生物机体某些物理现象的研究中 [1] 。由于其精确解往往不容易求得 [2] ，因此研究其数值求解方法无疑具有非常重要的理论意义和工程应用价值。

由现有文献知，热传导方程的数值求解方法主要有有限差分法、有限元法、谱方法等 [3] [4] [5] [6] ，其中应用最广泛的是差分法和边界元方法。差分法在求解问题时要求对整个区域作网格剖分，所以对复杂区域生成质量较好的网格非常耗时，特别是边界不规则时计算精度也不高。边界元方法要作复杂的奇异积分，费时费力。传统方法在一定程度上给出了满足精度要求的数值解，但它们都是网格方法，求解过程依赖于区域的网格剖分，本文考虑使用真正无网格方法数值求解热传导方程，希望能达到计算精度高，计算量小，使用方便，适应性强等好的效果。

2. 热传导方程的径向基无网格线方法

径向基线方法 [7] 是一种不依赖网格剖分的真正无网格方法。其基本思想是首先在离散点上运用径向基函数来近似空间导数，将偏微分方程问题转化为的常微分方程组问题，再利用常微分方程组的高精度数值方法进行求解，最终得到偏微分方程的高精度数值解。

下面考虑热传导方程的初边值问题

$\frac{\partial u\left(x,t\right)}{\partial t}=\alpha \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2},0\le x\le l,0\le t\le \infty$ (1)

$u\left(x,0\right)=\phi \left(x\right)$ (2)

$u\left(0,t\right)={\mu }_1\left(t\right),u\left(l,t\right)={\mu }_2\left(t\right)$ (3)

其中已知函数 $\phi \left(x\right)$ 表示初始温度分布情况，已知函数 ${\mu }_1\left(t\right)$ ， ${\mu }_{\text{2}}\left(t\right)$ 表示区间 $\left[0,l\right]$ 两个端点处温度随时间的变化情况。

首先，将求解区间进行离散，得到N个离散节点 $x_i,i=1,2,\cdots ,N$ ，其中 $x_i,i=2,3,\cdots ,N-1$ 是内部点，而 $x_1$ 和 $x_N$ 是边界点，借助径向基函数得到方程解得近似函数，用 $u^n\left(x\right)$ 表示为

$u^n\left(x\right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}{\lambda }_j^n\left(t\right){\psi }_j\left(x\right)}={\Psi }^{\text{T}}\left(x\right)\lambda$ (4)

其中 $\lambda$ 是与时间相关的未知量。

令

$u^n\left(x_i\right)=u_i$ (5)

引入向量符号 $u$ ，

$A\lambda =u$ (6)

其中

$u={\left[u_1\left(t\right),u_2\left(t\right),\cdots ,u_N\left(t\right)\right]}^{\text{T}}$ (7)

$\lambda ={\left[{\lambda }_1\left(t\right),{\lambda }_2\left(t\right),\cdots ,{\lambda }_N\left(t\right)\right]}^{\text{T}}$ (8)

则

$A=\left(\begin{array}{c}{\Psi }^{\text{T}}\left(x_{\text{1}}\right)\\ {\Psi }^{\text{T}}\left(x_{\text{2}}\right)\\ ⋮\\ {\Psi }^{\text{T}}\left(x_N\right)\end{array}\right)=\left(\begin{array}{cccc}{\psi }_1\left(x_1\right)& {\psi }_2\left(x_1\right)& \cdots & {\psi }_N\left(x_1\right)\\ {\psi }_1\left(x_{\text{2}}\right)& {\psi }_2\left(x_{\text{2}}\right)& \cdots & {\psi }_N\left(x_{\text{2}}\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\psi }_1\left(x_N\right)& {\psi }_2\left(x_N\right)& \cdots & {\psi }_N\left(x_N\right)\end{array}\right)$ (9)

由(5)和(6)可以得到

$u^n\left(x\right)={\Psi }^{\text{T}}\left(x\right)A^{-1}u=W\left(x\right)u$ (10)

其中

$W\left(x\right)={\Psi }^{\text{T}}\left(x\right)A^{-1}$ (11)

将(10)代入(1)，按照顺序将 $x_i$ 代入得到

$\frac{\text{d}{u}_i}{\text{d}t}-\alpha W_{xx}\left(x_i\right)u=0,i=1,2,\cdots ,N$ (12)

其中

$u_i\left(t\right)=u_i$ (13)

$W_x\left(x_i\right)=\left[W_{1x}\left(x_i\right)W_{2x}\left(x_i\right)\cdots W_{Nx}\left(x_i\right)\right]$ (14)

$W_{jxx}\left(x_i\right)=\frac{{\partial }^2}{\partial x^2}{W}_j\left(x_i\right),j=1,2,\cdots ,N$ (15)

为了把上面的系统写成列向量的形式，令

$U={\left[u_1{u}_2\cdots u_N\right]}^{\text{T}}$ (16)

$W_{xx}={\left[W_{jxx}\left(x_i\right)\right]}_{N\text{×}N}$ (17)

则(12)可改写为

$\frac{\text{d}U}{\text{d}t}-\alpha \left(W_{xx}U\right)=0$ (18)

令

$H\left(U\right)=\alpha \left(W_{xx}U\right)$ (19)

则

$\frac{\text{d}U}{\text{d}t}=H\left(U\right)$ (20)

初始条件为

$U^{\text{0}}={\left[\phi \left(x_1\right),\phi \left(x_2\right),\cdots ,\phi \left(x_N\right)\right]}^{\text{T}}$ (21)

又由(2)式给出边界条件

$u\left(0,t\right)={\mu }_1\left(t\right),u\left(l,t\right)={\mu }_2\left(t\right)$ (22)

通过上述方法可将热传导方程问题转化为常微分方程组问题，可以借助如下Runge-Kutta方法进行求解。

$\{\begin{array}{l}{U}^{n+1}=U^n+\frac{\Delta t\left(K_1+2\left(K_2+K_3\right)+K_4\right)}{6}\\ K_1=H\left(U^n\right),K_2=H\left(U^n+\frac{\Delta t}{2}{K}_1\right)\\ K_3=H\left(U^n+\frac{\Delta t}{2}{K}_2\right),K_4=H\left(U^n+\frac{\Delta t}{2}{K}_3\right)\end{array}$ (23)

其中 $\Delta t$ 是时间步长.在Runge-Kutta方法中，只要时间步长 $\Delta t$ 足够小，就能保证求解的稳定性 [2] 。

3. 数值实验

我们考虑热传导方程的初边值问题(1)~(3)，该问题的真解是

$u_z\left(x,t\right)={\text{e}}^{-{\pi }^{2}t}\sin \text{π}x$ .

为衡量数值解的精度，定义误差如下：

$ABE=\underset{i}{\max }\left|\frac{u_z\left(x_i,t_i\right)-u\left(x_i,t_i\right)}{u_z\left(x_i,t_i\right)}\right|$

其中 $u_z\left(x_i,t_i\right)$ 和 $u\left(x_i,t_i\right)$ 分别是真解和数值解。

为了考察基于径向基函数的无网格线方法在热传导方程数值求解中的可行性和有效性，我们下面分别讨论T时刻的数值解情况(图像、误差)以及不同参数对数值解的影响情况。

首先，给出时刻T = 1，33个离散节点时数值解(图1)及其误差(图2)图像。

Figure 1. Comparison of exact and numerical solutions

图1. 真解数值解对比图

Figure 2. The distribution of maximum relative error

图2. 误差分布图

由图1可得，计算所得数值解与真解贴合程度较好,说明所求数值解精确度较高，能一定程度上满足人们对热传导方程解的需求。分析图2可知，数值解的最大相对误差在区间的两个端点处最大，在区间内部，误差较小且变化幅度不大，说明数值解在区间内部可以达到更高的精确。

其次，我们固定所有参数，考虑不同时刻的数值解的图像(图3)和最大相对误差(表1)情况。

Figure 3. Numerical solutions at different times

图3. 不同时刻的数值解图像

通过观察图3可以知道，数值解完美反映了不同时刻解曲线的特征：在区间 $x\in \left[0,1\right]$ 上的变化趋势为先增大后减小，关于 $x=0.5$ 对称，随着考虑时刻T的增大，曲线趋于平缓，弯曲程度逐渐减小。

表1的数据说明，当时刻T逐渐增大时，数值解的误差随之增大。但当时刻T小于1时，误差小于0.01，此时具有较高的精度。

最后，讨论比例系数和离散节点数目对数值解精度的影响情况。比例系数RAT定义为

$\text{RAT}=\frac{\Delta t}{h^2}$

其中 $\Delta t$ 是时间步长，h是等分节点的间距。

Table 1. Maximum relative errors at different times

表1. 不同时刻的误差情况

通过多次数值实验，我们发现RAT对数值解的稳定性和精度有较大影响，并当其超过某一定值时，数值解的稳定性无法保证，误差也随之骤然上升(如图4)。

Figure 4. The changing errors with proportional coefficient RAT

图4. 误差随比例系数RAT的变化

分析实验结果可得，RAT突变点在0.3附近取得，约为0.2818。故当T = 0.5时，RAT在 $\left(0,0.2818\right)$ 范围内变动时误差较小，数值解结果较精确。

表2列出了固定比例系数RAT的值为0.02，时间T的值为1，进行数值实验后得到的最大相对误差随节点数目的变化情况。由表2可知，当N以10为步长增加时，数值解误差随之减小，并当N = 31时误差小于0.01，达到较精确的程度。

4. 结论

文章借助径向基函数无网格线方法数值求解一维热传导方程的初边值问题，在空间上使用径向基函

Table 2. Maximum relative errors with different nodes

表2. 不同离散节点对应的误差

数进行逼近，将偏微分方程问题转化为常微分方程组问题，进而利用Runge-Kutta方法得到高精度数值解。数值实验表明径向基无网格方法步骤简洁、易于操作，实验结果说明了该方法的可行性和有效性。因此，有必要进一步推广和完善该方法。

基金项目

大学生创新创业训练计划项目(201810452015)，山东省自然科学基金项目(ZR2018LA008)。

NOTES

^*通讯作者。