ABSTRACT

In order to investigate the numerical solution approximated to the exact solution of the initial boundary value problems of the heat transfer equations, the cubic spline solution is given here. Using the parametric cubic spline function based on Legendre-Gauss-Lobatto node for the numerical solution of Dirichlet boundary value problems of the heat transfer equations, the spatial direction is discretized by cubic spline function with parameters, and the parameters are appropriately chosen to improve accuracy of numerical errors and use the Crank-Nicolson format to construct the algorithm scheme in the time direction. Finally, the algorithm and numerical examples are given to demonstrate the high efficiency and accuracy of the proposed algorithm.

Keywords:Heat Transfer Equation, Initial Boundary Value Problems, Legendre-Gauss-Lobatto Nodes, Parametric Cubic Spline Function

热传导方程初边值问题的三次样条解

刘丹，王天军

河南科技大学数学与统计学院，河南 洛阳

收稿日期：2019年7月20日；录用日期：2019年8月7日；发布日期：2019年8月14日

摘 要

本文利用Legendre-Gauss-Lobatto节点为配置点，构造含参数的三次样条函数来数值求解带Dirichlet边界条件的热传导方程。其中，空间方向用含参数的三次样条函数离散，选取适当的参数值以提高数值误差的精度；时间方向用Crank-Nicolson格式离散，构造算法格式。并给出相应的算法格式和数值算例，表明该算法格式的有效性和高精度。

关键词 :热传导方程，初边值问题，Legendre-Gauss-Lobatto节点，含参数的三次样条函数

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

求解热传导方程初边值问题的精确解，通常用分离变量法 [1] 。但是寻求精确解往往比较困难，因此，一般是寻求逼近精确解的数值解。文献 [2] 中关于如何用差分法来求解热传导方程初值以及初边值问题的数值解有详细介绍；文献 [3] [4] 是用Lagrange插值来求解非线性热传导方程数值解；文献 [5] [6] [7] 则是用谱方法来逼近求解一些热传导方程：文献 [5] 用Cheyshev谱配置法来求解热传导方程初边值问题；文献 [6] 用Laguerre拟谱方法来求解半无界非线性热传导方程问题；文献 [7] 用Legendre-Hermite混合谱方法数值求解无穷带状区域上的各向异性线性热传导方程。最近，有作者用含参数的三次样条插值来求解常微分方程 [8] [9] [10] ，但是还未曾有人将其用来求解偏微分方程，本文用含参数三次样条插值求线性热传导方程的数值解，通过适当选取参数提高数值误差精度。

2. 算法格式

2.1. Legendre-Gauss-Lobatto节点

n次Legendre多项式的定义如下：

$L_n\left(t\right)=\frac{{\left(-1\right)}^n}{2^{n}n!}\frac{{\text{d}}^n\left({\left(1-t^2\right)}^n\right)}{\text{d}{t}^n},n=0,1,\cdots ,t\in \left[-1,1\right].$

选取 $t_k\left(0\le k\le n\right)$ 为多项式 $\left(1-t^2\right){\partial }_t{L}_n\left(t\right)$ 的 $n+1$ 个互异的零点，也即为Legendre-Gauss-Lobatto节点。

以Legendre-Gauss-Lobatto节点为配置点来构造含参数三次样条函数，然后用含参数的三次样条函数来逼近求解热传导方程初边值问题。

2.2. 含参数的三次样条函数

使用变量变换 $x=x\left(t\right)$ ， $t\in \left[-1,1\right]$ ，将求解区间 $\left[a,b\right]$ 映射到区间 $\left[-1,1\right]$ 。并将区间 $\left[a,b\right]$ 划分成n个小的子区间 $\left[x_i,x_{i+1}\right],x_i=x\left(t_i\right),h_{i+1}=x_{i+1}-x_i,a=x_0<\cdots <x_i<\cdots <x_n=b,i=0,1,\cdots ,n-1$ 。

一般来说，对函数 $u\left(x\right)$ 在插值节点 $x_k\left(0\le k\le n\right)$ 处进行插值形成的含参三次样条插值函数的好坏程度，即 $S\left(x,\xi \right)\left(S\left(x,\xi \right)\in C^2\left[a,b\right]\right)$ 逼近函数 $u\left(x\right)$ 的好坏程度，取决于参数 $\xi \left(\xi >0\right)$ 的适当选取以及当参数 $\xi \to 0$ 时三次样条函数在区间 $\left[a,b\right]$ 上的改变程度。

该样条函数 $S\left(x,\xi \right)=S\left(x\right)$ 在区间 $\left[x_i,x_{i+1}\right]$ 上满足以下的微分方程 [9] ：

$S^{″}\left(x\right)+\xi S\left(x\right)=\left(S^{″}\left(x_i\right)+\xi S\left(x_i\right)\right)\frac{x_{i+1}-x}{h_{i+1}}+\left(S^{″}\left(x_{i+1}\right)+\xi S\left(x_{i+1}\right)\right)\frac{x-x_i}{h_{i+1}},i=0,1,\cdots ,n-1.$ (1)

公式(1)是一类二阶常系数非齐次常微分方程。它的通解与相对应的齐次方程的通解以及它自己本身的一个特解有关。

显然，与公式(1)相对应的齐次方程的通解如下：

$S\left(x\right)=C_1\cos \sqrt{\xi }x+C_2\sin \sqrt{\xi }x.$

同时，有公式(1)的一个特解如下：

$S^{*}\left(x\right)=\left(S^{″}\left(x_i\right)+\xi S\left(x_i\right)\right)\frac{x_{i+1}-x}{\xi h_{i+1}}+\left(S^{″}\left(x_{i+1}\right)+\xi S\left(x_{i+1}\right)\right)\frac{x-x_i}{\xi h_{i+1}}.$

令 ${\omega }_{i+1}=h_{i+1}{\xi }^{1/2}$ ，再由插值条件： $S\left(x_i\right)=u_i,S\left(x_{i+1}\right)=u_{i+1}$ ，则可以得到公式(1)的通解，如下：

$\begin{array}{c}S\left(x\right)=-\frac{h_{i+1}^2}{{\omega }_{i+1}^2\sin {\omega }_{i+1}}\left[S^{″}\left(x_{i+1}\right)\sin \left(\frac{{\omega }_{i+1}\left(x-x_i\right)}{h_{i+1}}\right)+S^{″}\left(x_i\right)\sin \left(\frac{{\omega }_{i+1}\left(x_{i+1}-x\right)}{h_{i+1}}\right)\right]\\ +\frac{h_{i+1}^2}{{\omega }_{i+1}^2}\left[\frac{x-x_i}{h_{i+1}}\left(S^{″}\left(x_{i+1}\right)+\frac{{\omega }_{i+1}^2}{h_{i+1}^2}S\left(x_{i+1}\right)\right)+\frac{x_{i+1}-x}{h_{i+1}}\left(S^{″}\left(x_i\right)+\frac{{\omega }_{i+1}^2}{h_{i+1}^2}S\left(x_i\right)\right)\right].\end{array}$ (2)

对公式(2)两端进行求导并令 $x=x_i$ ，则有：

$S^{\prime }\left(x_i+0\right)=-\frac{S^{″}\left(x_{i+1}\right)-S^{″}\left(x_i\right)\cos {\omega }_{i+1}}{{\xi }^{1/2}\sin {\omega }_{i+1}}+\frac{S^{″}\left(x_{i+1}\right)-S^{″}\left(x_i\right)}{\xi h_{i+1}}+\frac{S\left(x_{i+1}\right)-S\left(x_i\right)}{h_{i+1}};$

当 $x\in \left[x_{i-1},x_i\right]$ 时，有 $h_i=x_i-x_{i-1}$ ， ${\omega }_i=h_i{\xi }^{1/2}$ 。类似上面步骤，可以得到：

$S^{\prime }\left(x_i-0\right)=-\frac{S^{″}\left(x_i\right)\cos {\omega }_i-S^{″}\left(x_{i-1}\right)}{{\xi }^{1/2}\sin {\omega }_i}+\frac{S^{″}\left(x_i\right)-S^{″}\left(x_{i-1}\right)}{\xi h_i}+\frac{S\left(x_i\right)-S\left(x_{i-1}\right)}{h_i}.$

由 $S^{\prime }\left(x\right)$ 在节点 $x_i$ 处的连续性条件，可以得到：

${\alpha }_i{M}_{i-1}+\left({\beta }_{i+1}+{\beta }_i\right)M_i+{\alpha }_{i+1}{M}_{i+1}=\frac{u_{i+1}}{h_{i+1}}-\left(\frac{1}{h_{i+1}}+\frac{1}{h_i}\right)u_i+\frac{u_{i-1}}{h_i},i=1,\cdots ,n-1.$ (3)

其中：

${\alpha }_i=\frac{h_i\left({\omega }_i\csc {\omega }_i-1\right)}{{\omega }_i^2}=\frac{1}{{\xi }^{1/2}\sin {\omega }_i}-\frac{1}{\xi h_i},S\left(x_i\right)=u_i;$

${\beta }_i=\frac{h_i\left(1-{\omega }_i\cot {\omega }_i\right)}{{\omega }_i^2}=\frac{1}{\xi h_i}-\frac{\cot {\omega }_i}{{\xi }^{1/2}},S^{″}\left(x_i\right)=M_i.$

2.3. 热传导方程算法格式

带Dirichlet边界条件热传导方程：

$\{\begin{array}{l}{u}_t=u_{xx}-f\left(x,t\right),0<x<l,t>0;\\ {u|}_{t=0}=\phi \left(x\right),0\le x\le l;\\ {u|}_{x=0}=\mu \left(t\right),{u|}_{x=l}=v\left(t\right),t>0.\end{array}$ (4)

用三次样条函数 $S\left(x\right)\left(x\in \left[0,l\right]\right)$ 来逼近方程(4)的解，由(4)式可得：

$\frac{\partial u_i\left(t\right)}{\partial t}=\frac{{\partial }^2{u}_i\left(t\right)}{\partial x^2}-f_i\left(t\right)=M_i\left(t\right)-f_i\left(t\right),i=1,2,\cdots ,n-1.$

其中： $f_i\left(t\right)=f\left(x_i,t\right);u_i\left(t\right)=u\left(x_i,t\right);M_i\left(t\right)=u_{xx}\left(x_i,t\right),i=0,1,\cdots ,n$ 。

考虑 $\frac{\partial u_i\left(t\right)}{\partial t}$ 前面的系数，有：

$M_i=q_i\frac{\partial u_i\left(t\right)}{\partial t}+f_i\left(t\right).$ (5)

记 $\frac{\partial u_i\left(t\right)}{\partial t}$ 为 ${u^{\prime }}_i\left(t\right)$ ，并将(5)代入(3)中得到：

$\begin{array}{l}\left({\alpha }_i{q}_{i-1}{u^{\prime }}_{i-1}+\left({\beta }_i+{\beta }_{i+1}\right)q_i{u^{\prime }}_i+{\alpha }_{i+1}{q}_{i+1}{u^{\prime }}_{i+1}\right)+\left({\alpha }_i{f}_{i-1}+\left({\beta }_i+{\beta }_{i+1}\right)f_i+{\alpha }_{i+1}{f}_{i+1}\right)\\ =\frac{u_{i-1}}{h_i}-\left(\frac{1}{h_i}+\frac{1}{h_{i+1}}\right)u_i+\frac{u_{i+1}}{h{}_{i+1}},i=1,2,\cdots ,n-1.\end{array}$ (6)

由(4)中给出的边界条件，处理得： ${u^{\prime }}_0\left(t\right)={\mu }^{\prime }\left(t\right),{u^{\prime }}_n\left(t\right)=v^{\prime }\left(t\right)$ ，则(6)可写成如下形式：

$\begin{array}{l}\left(\left({\beta }_1+{\beta }_2\right)q_1{u^{\prime }}_1+{\alpha }_2{q}_2{u^{\prime }}_2\right)+\left({\alpha }_1{f}_0+\left({\beta }_1+{\beta }_2\right)f_1+{\alpha }_2{f}_2+{\alpha }_1{q}_0{u^{\prime }}_0-\frac{u_0}{h_1}\right)\\ =-\left(\frac{1}{h_1}+\frac{1}{h_2}\right)u_1+\frac{u_2}{h_2},i=1;\end{array}$ (7)

$\begin{array}{l}\left({\alpha }_i{q}_{i-1}{u^{\prime }}_{i-1}+\left({\beta }_i+{\beta }_{i+1}\right)q_i{u^{\prime }}_i+{\alpha }_{i+1}{q}_{i+1}{u^{\prime }}_{i+1}\right)\\ +\left({\alpha }_i{f}_{i-1}+\left({\beta }_i+{\beta }_{i+1}\right)f_i+{\alpha }_{i+1}{f}_{i+1}\right)\\ =\frac{u_{i-1}}{h_i}-\left(\frac{1}{h_i}+\frac{1}{h_{i+1}}\right)u_i+\frac{u_{i+1}}{h_{i+1}},i=2,\cdots ,n-2;\end{array}$ (8)

$\begin{array}{l}\left({\alpha }_{n-1}{q}_{n-2}{u^{\prime }}_{n-2}+\left({\beta }_{n-1}+{\beta }_n\right)q_{n-1}{u^{\prime }}_{n-1}\right)+\left({\alpha }_{n-1}{f}_{n-2}+\left({\beta }_{n-1}+{\beta }_n\right)f_{n-1}+{\alpha }_n{f}_n+{\alpha }_n{q}_n{u^{\prime }}_n-\frac{u_n}{h_n}\right)\\ =\frac{u_{n-2}}{h_{n-1}}-\left(\frac{\text{1}}{h_n}+\frac{\text{1}}{h_{n-1}}\right)u_{n-1},i=n-1.\end{array}$ (9)

由(4)的初值条件： $u\left(x,0\right)=\varphi \left(x\right)\left(0\le x\le l\right)$ ，得： $u\left(x_k,0\right)=\varphi \left(x_k\right),k=0,1,\cdots ,n$ 。

为方便处理，记： $X\left(t\right)={\left(u_1\left(t\right),u_2\left(t\right),\cdots ,u_{n-1}\left(t\right)\right)}^{\text{T}}$ ； $X_0={\left(\varphi \left(x_1\right),\varphi \left(x_2\right),\cdots ,\varphi \left(x_{n-1}\right)\right)}^{\text{T}}$ ；

$A=\left[\begin{array}{cccc}{q}_1\left({\beta }_1+{\beta }_2\right)& q_2{\alpha }_2& & \\ q_1{\alpha }_2& q_2\left({\beta }_2+{\beta }_3\right)& q_3{\alpha }_3& \\ \ddots & \ddots & \ddots & \\ & q_{n-3}\alpha {}_{n-2}& q_{n-2}\left({\beta }_{n-2}+{\beta }_{n-1}\right)& q_{n-1}{\alpha }_{n-1}\\ & & q_{n-2}{\alpha }_{n-1}& q_{n-1}\left({\beta }_{n-1}+{\beta }_n\right)\end{array}\right];$

$B=\left[\begin{array}{cccc}-\left(\frac{1}{h_1}+\frac{1}{h_2}\right)& \frac{1}{h_2}& & \\ \frac{1}{h_2}& -\left(\frac{1}{h_2}+\frac{1}{h_3}\right)& \frac{1}{h_3}& \\ \ddots & \ddots & \ddots & \\ & \frac{1}{h_{n-2}}& -\left(\frac{1}{h_{n-2}}+\frac{1}{h_{n-1}}\right)& \frac{1}{h_{n-1}}\\ & & \frac{1}{h_{n-1}}& -\left(\frac{1}{h_{n-1}}+\frac{1}{h_n}\right)\end{array}\right];$

$F\left(t\right)=\left[\begin{array}{c}{\alpha }_1{f}_0+\left({\beta }_1+{\beta }_2\right)f_1+{\alpha }_2{f}_2+{\alpha }_1{q}_0{u^{\prime }}_0-\frac{u_0}{h_1}\\ {\alpha }_2{f}_1+\left({\beta }_2+{\beta }_3\right)f_2+{\alpha }_3{f}_3\\ ⋮\\ {\alpha }_i{f}_{i-1}+\left({\beta }_i+{\beta }_{i+1}\right)f_i+{\alpha }_{i+1}{f}_{i+1}\\ ⋮\\ {\alpha }_{n-2}{f}_{n-3}+\left({\beta }_{n-2}+{\beta }_{n-1}\right)f_{n-2}+{\alpha }_{n-1}{f}_{n-1}\\ {\alpha }_{n-1}{f}_{n-2}+\left({\beta }_{n-1}+{\beta }_n\right)f_{n-1}+{\alpha }_n{f}_n+{\alpha }_n{q}_n{u^{\prime }}_n-\frac{u_n}{h_n}\end{array}\right].$

有了以上的符号说明，则(7)~(9)可以写成如下形式：

$\{\begin{array}{l}A\frac{\text{d}X\left(t\right)}{\text{d}t}=BX\left(t\right)-F\left(t\right),t>0;\\ X\left(0\right)=X_0.\end{array}$ (10)

3. 数值例子

用式(10)求解式(4)，在时间方向用步长为 $\tau$ 的Crank-Nicolson格式离散，则：

$\{\begin{array}{l}\left(I-\left(\frac{\tau }{2}\right)A^{-1}B\right)X\left(t+\tau \right)=\left(I+\left(\frac{\tau }{2}\right)A^{-1}B\right)X\left(t\right)-\left(\frac{\tau }{2}\right)A^{-1}\left(F\left(t\right)+F\left(t+\tau \right)\right),t=0,\tau ,2\tau ,\cdots ,T-\tau ;\\ X\left(0\right)=X_0.\end{array}$

用如下定义的 $L^{\infty }$ -范数来度量误差：

$E_{N,\tau ,\xi }=\underset{0\le m\le N}{\max }\left|u\left(x_m,t\right)-u_N\left(x_m,t\right)\right|.$

$E{1}_{N,\tau ,\xi }=\underset{0\le m\le N}{\max }\left|u^{\prime }\left(x_m,t\right)-{u^{\prime }}_N\left(x_m,t\right)\right|.$

$E{2}_{N,\tau ,\xi }=\underset{0\le m\le N}{\max }\left|u^{″}\left(x_m,t\right)-{u^{″}}_N\left(x_m,t\right)\right|.$

其中 $E_{N,\tau ,\xi }$ 为三次样条误差， $E{1}_{N,\tau ,\xi }$ 为一阶导数误差， $E{2}_{N,\tau ,\xi }$ 为二阶导数误差。下面给出一些数值结果，以验证该格式的有效性和高精度。

例1：考虑以下微分方程：

$\{\begin{array}{l}{u}_t=u_{xx}-7{\text{e}}^{2x-3t},0<x<1,0<t<1;\\ {u|}_{t=0}={\text{e}}^{2x},0\le x\le 1;\\ {u|}_{x=0}={\text{e}}^{-3t},{u|}_{x=1}={\text{e}}^{-3t+2},0<t<1.\end{array}$ (11)

现有微分方程(11)的一个精确解： $u={\text{e}}^{2x-3t}$ 。

当三次样条参数 $\xi$ 变化时，图1给出 $L^{\infty }$ -范数 $E_{N,\tau ,\xi }$ 的常用对数 ${\log }_{10}{E}_{N,\tau ,\xi }$ 及不同 ${\log }_{10}N$ 之间的关系，

其中 $\tau$ 取值为0.001。可以看到数值误差随N的增加而快速衰减，当参数 $\xi$ 取0.1时，可以获得较好的数

值结果。在图2中，给出 $T=1$ 时数值误差常用对数 ${\log }_{10}{E}_{N,\tau ,\xi }$ 随N的增加及不同的 $\tau$ 的变化情况。可以

看到，随着 $\tau$ 的减小，数值误差随N的增加而快速衰减，这表明计算格式的高精度和稳定性。

图1. ${\log }_{10}{E}_{N,\tau ,\xi }:\tau =0.001$

图2. ${\log }_{10}{E}_{N,\tau ,\xi }:\xi =0.1,T=1$

图3给出了一阶导数误差，图4给出二阶导数误差，其中 $\xi =0.1,\tau =0.001$ 。可以看到，用三次样条函数来逼近求解方程(11)时，其一阶导数和二阶导数也分别收敛于精确解的一阶导数与二阶导数，并且具有很好的收敛速度。

图3. ${\log }_{10}E{1}_{N,\tau ,\xi }:\xi =0.1,\tau =0.001$

图4. ${\log }_{10}E{2}_{N,\tau ,\xi }:\xi =0.1,\tau =0.001$

例2：考虑以下微分方程：

$\{\begin{array}{l}{u}_t=u_{xx}-2\sin \left(x\right){\text{e}}^{-3t},0<x<\frac{\text{π}}{2},0<t<1;\\ {u|}_{t=0}=\sin x,0\le x\le \frac{\text{π}}{2};\\ {u|}_{x=0}=0,{u|}_{x=\text{π}/2}={\text{e}}^{-3t},0<t<1.\end{array}$ (12)

现有微分方程(12)的一个精确解： $u=\sin \left(x\right){\text{e}}^{-3t}$ 。

图5给出 $\tau$ 取值为0.0001时 ${\log }_{10}{E}_{N,\tau ,\xi }$ 与 ${\log }_{10}N$ 之间的关系。可以看到数值误差随N的增加而快速衰减，而当参数 $\xi$ 取0.6时，获得较好的数值结果。在图6中，T取值为1：随着 $\tau$ 的减小，数值误差随N的增加而快速衰减，这表明计算格式的高精度和稳定性。图7、图8分别给出了一阶导数误差、二阶导数误差，其中 $\xi =0.6,\tau =0.0001$ 。同样可以看到，用三次样条函数来逼近求解方程(12)时，其一阶导数和二阶导数也分别收敛于精确解的一阶导数与二阶导数，并且具有很好的收敛速度。

图5. ${\log }_{10}{E}_{N,\tau ,\xi }:\tau =0.0001$

图6. ${\log }_{10}{E}_{N,\tau ,\xi }:\xi =0.6,T=1$

图7. ${\log }_{10}E{1}_{N,\tau ,\xi }:\xi =0.6,\tau =0.0001$

图8. ${\log }_{10}E{2}_{N,\tau ,\xi }:\xi =0.6,\tau =0.0001$

4. 结论

本文利用含参数的三次样条插值函数来逼近热传导方程精确解。以Legendre-Gauss-Lobatto节点为配置点，先构造含参数的三次样条插值来对空间进行离散，再在时间方向上应用步长为 $\tau$ 的Crank-Nicolson格式，然后得到带Dirichlet边界条件热传导方程的数值解。由数值实验得到的结果来看，提出的算法不仅具有一定的有效性和高精度，并且对于精确解的导数也有很好的逼近效果。

致谢

首先，我要感谢基金项目：国家自然科学基金项目(11371123)！最后，对于师姐柴果、师妹马亚楠以及实验室的同学们表示感谢！谢谢你们对我的支持和帮助！

基金项目

国家自然科学基金项目(11371123)。

文章引用

刘 丹,王天军. 热传导方程初边值问题的三次样条解
Cubic Spline Solution to Initial Boundary Value Problems of Heat Transfer Equations[J]. 应用数学进展, 2019, 08(08): 1375-1383. https://doi.org/10.12677/AAM.2019.88161

参考文献