常微分方程初值问题的区域分解逼近方法

doi:10.12677/aam.2025.145251

期刊菜单

常微分方程初值问题的区域分解逼近方法
Domain Decomposition Approximation of Initial Value Problems for Ordinary Differential Equations

DOI: 10.12677/aam.2025.145251, PDF, HTML, XML, 科研立项经费支持
作者: 宋灵杉, 梁同洋, 刘工业：河南科技大学数学与统计学院，河南洛阳
关键词: 常微分方程；初值问题；拉格朗日插值法；多区域配置法；Ordinary Differential Equations； Initial Value Problem； Lagrange Interpolation Method； Collocation of Multi-Domain

摘要: 本文主要借助拉格朗日插值逼近来求解常微分方程初值问题的近似解，利用插值基函数展开数值解，将问题转化为线性或非线性方程组，对非线性方程组采用不动点迭代法求解。数值实验表明，多区域方法比单区域方法误差更小。

Abstract: This article mainly uses Lagrangian interpolation approximation to solve the approximation solution of the initial value problem of ordinary differential equations. The numerical solution is expanded by using interpolation basis functions, and the problem is transformed into a system of linear or nonlinear equations. The fixed-point iteration method is used to solve the system of nonlinear equations. Numerical experiments show that the multi-domain method has a smaller error than the single-region method.

文章引用：宋灵杉, 梁同洋, 刘工业. 常微分方程初值问题的区域分解逼近方法[J]. 应用数学进展, 2025, 14(5): 228-237. https://doi.org/10.12677/aam.2025.145251

1. 引言

常微分方程初值问题是数学和工程学的一个基本问题，它描述了系统状态随时间的演化规律。这类问题在自然科学和工程技术中有着广泛的应用，如物理、化学、生物学、经济学等领域。因此，常微分方程的背景意义不仅仅体现在理论的研究当中，还被广泛地作用于对生活和工作的实际问题的建模和求解过程中，同样在科学研究与工程技术也都具有非常重要的价值。随着科学和技术的进步，对这类问题的求解方法和效率提出了更高的要求。常微分方程初值问题仍有许多问题需解决、算法也有待改进，因此在数学界关注度也比较大。通常学者们会使用数值方法如有限元方法、辛普森法则、微分方程求解器来解决常微分方程的初值问题。有学者采用Milstein方法[1]，也有学者采用Birkhoff配点法[2]，利用泰勒展开式求解二阶常微分方程[3]。最近，有人利用Lagrange插值[4]函数的优点，采用等距节点来求解偏微分方程[5]单区域时空方向上的数值解，他们的研究结果表明，这些方法都能够有效地进行初值问题的求解，但是仍然没有做出最好的结果，好在初值问题目前备受关注，学者们也都有提出不同的效果、更好的解决方法，拉格朗日插值法求近似解是误差较小的方法之一。我们常用的插值方法[6]主要有：拉格朗日插值法、牛顿插值法、埃尔米特插值法、分段插值法、样条插值法。插值法多种多样，不同插值法使用的条件和适用场景都不同，本文主要采用拉格朗日插值法来解决常微分方程的初值问题并且通过区域分解的方式进而提高精度、减小误差。

先通过构造单区域上的Lagrange插值逼近算法来逼近上述方程的精确解，再将区间(0, 2)分解成两部分，构造多区域上的Lagrange插值逼近算法，通过数值实验发现，该算法所求的近似解能够很好地逼近精确解，而且该算法实施较为容易，且易于理解。所提方法同样适用于求解刚性方程和二阶问题，见例2和例3。

2. 基于等距节点的拉格朗日多项式的微分矩阵

对于任意次数以 $a_{i}$ 为节点， $n$ 次插值基函数为：

$\begin{matrix} l_{i} (a_{k}) = \frac{(a_{k} - a_{0}) (a_{k} - a_{1}) \dots (a_{k} - a_{i - 1}) (a_{k} - a_{i + 1}) \dots (a_{k} - a_{n})}{(a_{i} - a_{0}) (a_{i} - a_{1}) \dots (a_{i} - a_{i - 1}) (a_{i} - a_{i + 1}) \dots (a_{i} - a_{n})} \\ = \frac{(a_{k} - a_{0}) (a_{k} - a_{1}) \dots (a_{k} - a_{n})}{(a_{i} - a_{0}) (a_{i} - a_{1}) \dots (a_{i} - a_{i - 1}) (a_{i} - a_{i + 1}) \dots (a_{i} - a_{n}) (a_{k} - a_{i})}, (k \neq i), k = 0, 1, \dots, N . \end{matrix}$

满足：

$l_{i} (a_{k}) = {\begin{matrix} 0, k \neq i; \\ 1, k = i . \end{matrix}$

引入 $ω_{n + 1} (a) = (a - a_{0}) (a - a_{1}) \dots (a - a_{n})$ ，

则 $ω'_{n + 1} (a_{i}) = (a_{i} - a_{0}) \dots (a_{i} - a_{k - 1}) (a_{i} - a_{k + 1}) \dots (a_{i} - a_{n})$ 。

所以 $n$ 次插值多项式为：

$L_{n} (a) = \sum_{i = 0}^{n} y_{i} \frac{ω_{n + 1} (a)}{(a - a_{i}) ω'_{n + 1} (a_{i})}$

$\partial_{a} l_{i} (a_{k}) = \frac{(a_{k} - a_{0}) (a_{k} - a_{1}) \dots (a_{k} - a_{k - 1}) (a_{k} - a_{k + 1}) \dots (a_{k} - a_{N})}{(a_{i} - a_{0}) (a_{i} - a_{1}) \dots (a_{i} - a_{i - 1}) (a_{i} - a_{i + 1}) \dots (a_{i} - a_{N}) (a_{k} - a_{i})}$

设 $x$ 为等距节点， $a_{k + 1} - a_{k} = h, h = \frac{1}{N}, k = 0, 1, \dots, N .$

$\partial_{x} L_{N} (a_{k}) = \sum_{j = 0}^{N} d_{k i} y_{i}$

这里， $D = (d_{k i})$ 是 $(N + 1) \times (N + 1)$ 级的矩阵，且 $d_{k i} = \partial_{a} l_{i} (a_{k})$ ，由文献[5]可知：

$d_{k i} = {\begin{cases} \frac{{(- 1)}^{i - k} k! (N - k)!}{h (k - i) i! (N - i)!}, k \neq i, \\ \frac{1}{h} ((1 + \frac{1}{2} + \dots + \frac{1}{i}) - (1 + \frac{1}{2} + \dots + \frac{1}{N - i})), 0 < k = i < N, \\ - \frac{1}{h} (1 + \frac{1}{2} + \dots + \frac{1}{N}), k = i = 0, \\ \frac{1}{h} (1 + \frac{1}{2} + \dots + \frac{1}{N}), k = i = N . \end{cases}$

对于 $d_{k i}$ 的具体推导在文献[5]中有详细推导，由文献[7]可知，m阶微分矩阵与一阶微分矩阵的关系是： $D^{m} = D^{(m)}$ 。

3. 常微分方程多区域Lagrange插值逼近算法格式

这一章主要考虑常微分方程初值问题多区域配置法。为简单利用等距节点的Lagrange插值逼近，将所求问题转化为线性或非线性方程组求解。

3.1. 单区域

首先考虑单区域方法。利用插值基函数展开数值解，得到以展开系数为未知量的方程组。然后，根据方程组是线性或非线性选择相应的求解方法。

3.1.1. 基于单区域常微分方程拉格朗日插值逼近的算法推导

例1：

${\begin{cases} y^{2} - y' - 1 = 0, x \in (0, 2) \\ y (0) = 0 \end{cases}$ (1)

下面是关于例1的基于单区域一阶非线性常微分方程拉格朗日插值逼近的算法推导：

令 $x_{j} = j h, h = \frac{1}{M}, j = 0, 1, \dots, M$ ，记 $l_{j} (x)$ 为 $x \in (0, 2)$ 上的等距插值基函数，用插值多项式 $y_{M} (x)$ 逼近精确解 $y (x)$ ，利用插值基函数将数值解，即插值多项式 $y_{M} (x)$ 展开为：

$y_{M} (x) = \sum_{j = 0}^{M} y (x_{j}) l_{j} (x) = \sum_{j = 0}^{M} {\hat{y}}_{j} l_{j} (x), y (x_{j}) = {\hat{y}}_{j}$ . (2)

代入(2)式得：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{j = 0}^{M} d_{k j} {\hat{y}}_{j} - 1 = 0, k = 1, 2, 3, \dots, M \\ y (0) = \sum_{j = 0}^{M} l_{j} (0) {\hat{y}}_{j} = 0 \end{array}$ (3)

展开得：

${\hat{y}}_{k}^{2} - 2 ({\hat{y}}_{0} d_{k 0} + {\hat{y}}_{1} d_{k 1} + \dots + {\hat{y}}_{M - 1} d_{k, M - 1} + {\hat{y}}_{M} d_{k, M}) - 1 = 0$ (4)

初值条件展开为：

${y_{0} = 0}$ (5)

将(5)式代入(4)式，化简并移项可得：

${\hat{y}}_{k}^{2} - 2 \sum_{j = 1}^{M} d_{k j} {\hat{y}}_{j} - 1 = 0$ (6)

可以得到：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{j = 1}^{M} d_{k j} {\hat{y}}_{j} - 1 = 0, k = 1, 2, 3, \dots, M \\ y_{0} = 0 \end{array}$ (7)

取 $k = 1, 2, 3, \dots, M$ ，则(7)式可以写成：

$2 (\begin{matrix} d_{1, 1} & \dots & d_{1, M} \\ ⋮ & ⋱ & ⋮ \\ d_{M, 1} & \dots & d_{M, M} \end{matrix}) (\begin{matrix} {\hat{y}}_{1} \\ ⋮ \\ {\hat{y}}_{M} \end{matrix}) - (\begin{matrix} {\hat{y}}_{1}^{2} \\ ⋮ \\ {\hat{y}}_{M}^{2} \end{matrix}) + b = 0$

不妨记：

$A = (\begin{matrix} d_{1, 1} & \dots & d_{1, M} \\ ⋮ & ⋱ & ⋮ \\ d_{M, 1} & \dots & d_{M, M} \end{matrix}), \bar{y} = (\begin{matrix} {\hat{y}}_{1} \\ ⋮ \\ {\hat{y}}_{M} \end{matrix}), {\bar{y}}^{2} = (\begin{matrix} {\hat{y}}_{1}^{2} \\ ⋮ \\ {\hat{y}}_{M}^{2} \end{matrix}), b = (\begin{matrix} 1 \\ ⋮ \\ 1 \end{matrix})$

可以得到如下非线性方程组：

$\begin{array}{l} \bar{y} = \frac{1}{2} A^{- 1} [{\bar{y}}^{2} - b] \\ {\bar{y}}_{n + 1} = \frac{1}{2} A^{- 1} [{\bar{y}}_{n}^{2} - b] \end{array}$ (8)

下面考虑一个刚性问题。可以证明下面例2为刚性问题。利用配置方法可以求得数值解，但误差没有一般方程求解的效果好，但仍不失为一个有效方法。节点取12时，误差达到10的负3次幂。

例2：

${\begin{array}{l} y^{'} = - 1000 (y - \cos x) + \sin x, x \in (0, 2) \\ y (0) = 0 \end{array}$ (9)

同理，可以得到：

令

$A_{1} = (\begin{matrix} d_{1, 1} & \dots & d_{1, M} \\ ⋮ & ⋱ & ⋮ \\ d_{M, 1} & \dots & d_{M, M} \end{matrix}), {\bar{y}}_{1} = (\begin{matrix} {\hat{y}}_{1} \\ ⋮ \\ {\hat{y}}_{M} \end{matrix}), c_{1} = 1000 (\begin{matrix} \cos x_{1} \\ ⋮ \\ \cos x_{M} \end{matrix}) + (\begin{matrix} \sin x_{1} \\ ⋮ \\ \sin x_{M} \end{matrix})$

则可以得到方程：

${\bar{y}}_{1} = {(1000 E + A_{1})}^{- 1} c_{1}$ (10)

其中， $E$ 表示为 $M$ 阶单位阵。

例3：

${\begin{array}{l} y^{″} - y = 3 e^{- 2 x}, x \in (- 1, 1) \\ y (- 1) = e^{2} \\ y (1) = e^{- 2} \end{array}$ (11)

同理可得：令

$A_{2} = (\begin{matrix} d_{1, 1}^{(2)} & \dots & d_{1, M - 1}^{(2)} \\ ⋮ & ⋱ & ⋮ \\ d_{M - 1, 1}^{(2)} & \dots & d_{M - 1, M - 1}^{(2)} \end{matrix}), {\bar{y}}_{2} = (\begin{matrix} {\hat{y}}_{1} \\ ⋮ \\ {\hat{y}}_{M - 1} \end{matrix}), c_{2} = 3 (\begin{matrix} e^{- 2 x_{1}} \\ ⋮ \\ e^{- 2 x_{M - 1}} \end{matrix}), H = y_{0} (\begin{matrix} d_{1, 0}^{(2)} \\ ⋮ \\ d_{M - 1, 0}^{(2)} \end{matrix}) + y_{M} (\begin{matrix} d_{1, M}^{(2)} \\ ⋮ \\ d_{M - 1, M}^{(2)} \end{matrix})$

则解得方程为：

${\bar{y}}_{2} = {(A_{2} - E)}^{- 1} (c_{2} - H)$ (12)

其中， $E$ 表示为 $M - 1$ 阶单位阵。

3.1.2. 基于单区域常微分方程拉格朗日差值逼近的数值解及误差分析

用 $L^{\infty}$ 来衡量精确解与数值解之间的误差：

$E_{M} = \max | y (x_{k}) - y_{M} (x_{k}) |, 1 \leq k \leq M .$

对例1，由(8)式，使用不动点迭代法可得出近似解 $y_{M} (x)$ ，而方程(1)的精确解为：

$y = \frac{1 - e^{x}}{1 + e^{x}}$

可以得到误差图(图1)，对最大误差 $E_{M}$ 取对数可以得到随着插值次数 $M$ 的增大，最大误差 $E_{M}$ 呈下降趋势，且误差效果良好，可以说明上述所提算法格式良好，具有高精度。

Figure 1. Single area error plot

图1. 单区域误差图

同理，对例2，由(10)式求得其数值解。方程(9)的精确解为：

$y_{1} = - \frac{999999}{1000001} e^{- 1000 x} + \frac{999999}{1000001} \cos x + \frac{2000}{1000001} \sin x$

图2是例2刚性方程单区域误差图，可以看到随着节点的增大，误差在减小。单区域实验显示，数值解在初始阶段能快速捕捉解的剧烈衰减特性。

Figure 2. Stiff Equation single area error plot

图2. 刚性方程单区域误差图

同理，对例3，由(12)式求得其数值解。方程(11)的精确解为：

$y_{2} = e^{- 2 x}$ .

Figure 3. Single area error plot 1

图3. 单区域误差图1

图3是例3二阶常微分方程单区域误差图，对最大误差 $E_{M}$ 取对数可以得到随着插值次数 $M$ 的增大，最大误差 $E_{M}$ 呈下降趋势，且误差效果良好，可以说明上述所提算法格式良好，具有高精度。

3.2. 多区域

3.2.1. 基于多区域常微分方程拉格朗日差值逼近的算法推导

为了能够得到更高的误差精度，我们将整个求解区域 $(0, 2)$ 分解成两个区域 $Ω_{1} = (0, 1), Ω_{2} = (1, 2)$ 。数值解展开为：

${y_{M} (x) |}_{Ω_{i}} = y_{M}^{i} (x) = \sum_{n = 0}^{M_{i}} y_{M}^{(i)} l_{n} (x), x \in Ω_{i}, i = 1, 2.$

由已知 $y \in H^{2} (Ω)$ ，在公共边界处有数值解的函数值和导数值相等，即：

函数值： $\sum_{n = 0}^{M_{1}} {\hat{y}}_{n}^{(1)} l_{n} (x_{M_{1}}) = \sum_{n = 0}^{M_{2}} {\hat{y}}_{n}^{(2)} l_{n} ({\tilde{x}}_{0})$ 。

由于未知量为 $M_{1} + M_{2} + 2$ 个可列方程，区域一为 $M_{1}$ 个、区域二为 $M_{2}$ 个函数值，故未知量大于方程数还需一个条件可加入：导数值相等。

导数值： $\sum_{n = 0}^{M_{1}} {\hat{y}}_{n}^{(1)} \partial_{x} l_{n} (x_{M_{1}}) = \sum_{n = 0}^{M_{2}} {\hat{y}}_{n}^{(2)} \partial_{x} l_{n} ({\tilde{x}}_{0})$ 。

由微分矩阵定义，将上述两式联立可以得到：

${\hat{y}}_{M_{1}}^{(1)} = {\hat{y}}_{0}^{(2)} = \frac{\sum_{n = 1}^{M_{2}} {\hat{y}}_{n}^{(2)} {\tilde{d}}_{0, n} - \sum_{n = 1}^{M_{1} - 1} {\hat{y}}_{n}^{(1)} d_{M_{1}, n}}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} .$ (13)

其中， $d_{k, m}$ 表示区域 $Ω_{1}$ 上的微分矩阵， ${\tilde{d}}_{k, m}$ 表示区域 $Ω_{2}$ 上的微分矩阵， $x_{k} \in Ω_{1}, k = 1, 2, \dots, M_{1},$

${\tilde{x}}_{k} \in Ω_{2}, k = 1, 2, \dots, M_{2} .$

在区域 $Ω_{1}$ 用数值解逼近区域 $Ω_{2}$ 的精确解得到：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{n = 0}^{M_{1}} {\hat{y}}_{n}^{(1)} \partial_{x} l_{n} (x_{k}) - 1 = 0, k = 1, 2, \dots, M_{1}, \\ \sum_{n = 0}^{M_{1}} {\hat{y}}_{n}^{(1)} l_{n} (x_{0}) = 0 \end{array}$ (14)

对(14)利用微分矩阵定义并将端点值提出可以得到：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{n = 1}^{M_{1}} {\hat{y}}_{n}^{(1)} d_{k, n} - 2 {\hat{y}}_{0}^{(1)} d_{k, 0} - 1 = 0 \\ k = 1, 2, \dots, M_{1} - 1 \end{array}$ (15)

将(13)式代入到(15)式中有：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{n = 1}^{M_{1} - 1} {\hat{y}}_{n}^{(1)} d_{k, n} - 1 = 2 y_{0}^{(1)} d_{k, 0} + 2 \frac{\sum_{n = 1}^{M_{2}} {\hat{y}}_{n}^{(2)} {\tilde{d}}_{0, n} - \sum_{n = 1}^{M_{1} - 1} {\hat{y}}_{n}^{(1)} d_{M_{1}, n}}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} d_{k, M_{1}} \\ - 2 \frac{{\hat{y}}_{0}^{(1)} d_{M_{1}, 0}}{d_{M_{1,} M_{1}} - {\tilde{d}}_{0, 0}} d_{k, M_{1}} \\ k = 1, 2, \dots, M_{1} - 1 \end{array}$ (16)

将(16)式展开，并记：

$G = (\begin{matrix} d_{1, 1} & \dots & d_{1, M_{1} - 1} \\ ⋮ & ⋱ & ⋮ \\ d_{M_{1} - 1, 1} & \dots & d_{M_{1} - 1, M_{1} - 1} \end{matrix}), x_{1} = (\begin{matrix} {\hat{y}}_{1}^{(1)} \\ ⋮ \\ {\hat{y}}_{M_{1} - 1}^{(1)} \end{matrix}), x_{1}^{2} = (\begin{matrix} {\hat{y}}_{1}^{2} \\ ⋮ \\ {\hat{y}}_{M_{1} - 1}^{2} \end{matrix}), x_{2} = (\begin{matrix} {\hat{y}}_{1}^{(2)} \\ ⋮ \\ {\hat{y}}_{M_{2}}^{(2)} \end{matrix})$

$\begin{array}{l} E_{1} = (\begin{matrix} d_{M_{1}, 1} & \dots & d_{M_{1}, M_{1} - 1} \end{matrix}), E_{2} = (\begin{matrix} {\tilde{d}}_{0, 1} & \dots & {\tilde{d}}_{0, M_{2}} \end{matrix}), \\ F_{1} = (\begin{matrix} d_{1, M_{1}} \\ ⋮ \\ d_{M_{1} - 1, M_{1}} \end{matrix}), b_{1} = (\begin{matrix} 1 \\ ⋮ \\ 1 \end{matrix}) . \end{array}$

则(16)式可以转化成下列矩阵方程：

$x_{1}^{2} - 2 G x_{1} - b_{1} = 2 \frac{E_{2} x_{2} - E_{1} x_{1}}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} F_{1}$ (17)

同理，对区域 $Ω_{2}$ 仿照区域 $Ω_{1}$ 可以得到：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{n = 0}^{M_{2}} {\hat{y}}_{n}^{(2)} \partial_{x} l_{n} ({\tilde{x}}_{k}) - 1 = 0 \\ k = 1, 2, \dots, M_{2} \end{array}$ (18)

对(18)式利用微分矩阵定义并将端点值提出可以得到：

${\begin{array}{l} {\hat{y}}_{k}^{2} - 2 \sum_{n = 1}^{M_{2}} {\hat{y}}_{n}^{(2)} {\tilde{d}}_{k, n} - 2 {\hat{y}}_{0}^{(2)} {\tilde{d}}_{k, 0} - 1 = 0 \\ k = 1, 2, \dots, M_{2} \end{array}$ (19)

将(13)代入(19)式得：

${\begin{cases} {\hat{y}}_{k}^{2} - 2 \sum_{n = 1}^{M_{2}} {\hat{y}}_{n}^{(2)} {\tilde{d}}_{k, n} - 1 = 2 \frac{\sum_{n = 1}^{M_{2}} {\hat{y}}_{n}^{(2)} {\tilde{d}}_{0, n} - \sum_{n = 1}^{M_{1} - 1} {\hat{y}}_{n}^{(1)} d_{M_{1}, n}}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} {\tilde{d}}_{k, 0} - 2 \frac{{\hat{y}}_{0}^{(1)} d_{M_{1}, 0}}{d_{M_{1,} M_{1}} - {\tilde{d}}_{0, 0}} d_{k, 0} \\ k = 1, 2, \dots, M_{2} \end{cases}$ (20)

将(20)式展开，并记：

$P = (\begin{matrix} {\tilde{d}}_{1, 1} & \dots & {\tilde{d}}_{1, M_{2}} \\ ⋮ & ⋱ & ⋮ \\ {\tilde{d}}_{M_{2}, 1} & \dots & {\tilde{d}}_{M_{2}, M_{2}} \end{matrix}), x_{2}^{2} = (\begin{matrix} {\hat{y}}_{1}^{2} \\ ⋮ \\ {\hat{y}}_{M_{2}}^{2} \end{matrix}), F_{2} = (\begin{matrix} {\tilde{d}}_{1, 0} \\ ⋮ \\ {\tilde{d}}_{M_{2}, 0} \end{matrix}), b_{2} = (\begin{matrix} 1 \\ ⋮ \\ 1 \end{matrix})$

转化成如下矩阵方程：

$x_{2}^{2} - 2 P x_{2} - b_{2} = 2 \frac{E_{2} x_{2} - E_{1} x_{1}}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} F_{2}$ (21)

联立(17)和(21)可以得到：

$(\begin{matrix} x_{1}^{2} \\ x_{2}^{2} \end{matrix}) - 2 (\begin{matrix} G & 0 \\ 0 & P \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) - (\begin{matrix} b_{1} \\ b_{2} \end{matrix}) = 2 \frac{1}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} (\begin{matrix} F_{1} \\ F_{2} \end{matrix}) (\begin{matrix} - E_{1} & E_{2} \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \end{matrix})$ (22)

(22)式可以等价于：

$(\begin{matrix} x_{1}^{2} \\ x_{2}^{2} \end{matrix}) - (\begin{matrix} b_{1} \\ b_{2} \end{matrix}) = 2 \frac{1}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} (\begin{matrix} F_{1} \\ F_{2} \end{matrix}) (\begin{matrix} - E_{1} & E_{2} \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) + 2 (\begin{matrix} G & 0 \\ 0 & P \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \end{matrix})$

提取公因式得：

$(\begin{matrix} x_{1}^{2} \\ x_{2}^{2} \end{matrix}) - (\begin{matrix} b_{1} \\ b_{2} \end{matrix}) = 2 (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) [\frac{1}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} (\begin{matrix} F_{1} \\ F_{2} \end{matrix}) (\begin{matrix} - E_{1} & E_{2} \end{matrix}) + (\begin{matrix} G & 0 \\ 0 & P \end{matrix})]$

不妨令 $I = \frac{1}{d_{M_{1}, M_{1}} - {\tilde{d}}_{0, 0}} (\begin{matrix} F_{1} \\ F_{2} \end{matrix}) (\begin{matrix} - E_{1} & E_{2} \end{matrix}) + (\begin{matrix} G & 0 \\ 0 & P \end{matrix}), \bar{X} = (\begin{matrix} x_{1} \\ x_{2} \end{matrix}), {\bar{X}}_{2} = (\begin{matrix} x_{1}^{2} \\ x_{2}^{2} \end{matrix}), \bar{B} = (\begin{matrix} b_{1} \\ b_{2} \end{matrix})$ ，

可以得到以下方程：

$\begin{array}{l} \bar{X} = 0.5 * I^{- 1} * [{\bar{X}}^{2} - \bar{B}] \\ {\bar{X}}_{n} = 0.5 * I^{- 1} * [{\bar{X}}_{n + 1}^{2} - \bar{B}] \end{array}$ (23)

3.2.2. 基于多区域常微分方程拉格朗日插值逼近的数值解及误差分析

对 ${\bar{X}}_{n} = 0.5 * I^{- 1} * [{\bar{X}}_{n + 1}^{2} - \bar{B}]$ 进行迭代法求出近似解 ${\bar{X}}_{n}$ ，然后用MATLAB编程计算多项式次数 $N$ 和误差的关系，误差表示为 $E_{k}$ ，近似解为 $y$ ，真解为 $x$ ，仍然用 $L^{\infty}$ 来衡量精确解与数值解之间的误差：

$E_{M} = \max | x - y |, 1 \leq k \leq M - 1$

而方程的精确解为：

$y = \frac{1 - e^{x}}{1 + e^{x}}$

能够取得误差图，对最大误差 $E_{M}$ 取对数可以得到随着插值次数 $M$ 的增大，最大的误差 $E_{M}$ 呈下降趋势，并且误差的结果比较好，能够说明以上所用的算法格式比较好，有高精度。其中， $E_{M 1}$ 表示区域一上的误差， $E_{M 2}$ 表示区域二上的误差， $E_{M}$ 使整个区域上的误差 $M = M_{1} + M_{2}$ ：

$\begin{array}{l} E_{M 1} = \max | u_{1} (x_{k}) - u_{M 1} (x_{k}) |, 1 \leq k \leq M_{1} - 1 \\ E_{M 2} = \max | u_{2} (x_{k}) - u_{M 2} (x_{k}) |, 1 \leq k \leq M_{2} - 1 \\ E_{M} = \max (E_{M 1}, E_{M 2}) \end{array}$

Figure 4. Multi region error plot 2

图4. 多区域误差图2

图4是多区域误差图，可以看出，随着次数M的增大，误差 $E_{k}$ 减小。由图可知，随着插值次数M的增加，最大误差减小，单区域精度在−8，多区域误差精度在−9，所以多区域明显好于单区域，故多区域拉格朗日插值效果好于单区域拉格朗日插值效果。

为节省篇幅，这里略去例2、例3的多区域数值结果，有兴趣的读者可以利用本文例1所提供的多区域方法去自行求解。

4. 结论

本文采用拉格朗日插值法对区间分解逼近求解常微分方程的数值解，将常微分方程化成矩阵求解，由数值实验结果可知，算法格式较好。为了继续提高算法精度，对区域进行分解逼近，然后分别在单区域和多区域上进行数值解求取，将数值解与原方程精确解作对比，利用MATLAB画出结果图，通过实验结果可知单区域上的误差要远大于多区域上的误差，故可以体现出多区域的优越性。

基金项目

河南省大学生创新创业训练计划项目(No. 2024232)；河南省大学生创新创业训练计划项目(No. 202310464051)。

参考文献

[1]	李焕荣. 随机常微分方程的几种数值求解方法及其应用[J]. 重庆工商大学学报(自然科学版), 2021, 38(6): 82-88.
[2]	庄清渠, 王金平. 四阶常微分方程Birkhoff配点法[J]. 华侨大学学报(自然科学版), 2018, 39(2): 306-311.
[3]	户永清, 陈正文. 一种二阶常微分方程的数值解法[J]. 四川文理学院学报, 2023, 33(5): 39-44.
[4]	唐仁献. Lagrange基本插值多项式的性质及应用[J]. 零陵学院学报, 1995(S1): 49-51.
[5]	乔炎, 王川, 王秦. Burgers方程混合问题的Lagrange插值逼近[J]. 应用数学进展, 2022, 11(5): 2507-2514.
[6]	李庆扬, 王能超, 易大义. 数值分析[M]. 第5版. 北京: 清华大学出版社, 2001: 1-326.
[7]	Shen, J., Tang, T. and Wang, L.L. (2011) Spectral Methods: Algorithms, Analysis and Applications. Springer. https://doi.org/10.1007/978-3-540-71041-7

为你推荐

友情链接