常微分方程的数值求解 | 一阶方程组和高阶方程

学术科技 2024-10-15 22:21 山东

点击上方「蓝字」关注我们

在数学与工程领域中，常微分方程（ODEs）是描述动态系统行为的重要工具。从物理系统的动力学模型到化学反应的动力学过程，常微分方程无处不在。然而，许多实际问题中的常微分方程往往难以获得解析解，这促使人们发展了各种数值方法来近似求解这些方程。

我们在前面已经介绍了一阶常微分方程初值问题的若干解法，若作进一步推广，如求解一阶方程组初值问题，或者高阶方程的初值问题，又该如何分析呢？本文将给出相关的答案。

一阶方程组初值问题

为简单起见，设有 2 个一阶方程组成的方程组初值问题：

\left\{\begin{array}{l} y_{1}^{\prime}=f\left(x, y_{1}, y_{2}\right), & a\lt x \leqslant b \tag{1}\\ y_{2}^{\prime}=g\left(x, y_{1}, y_{2}\right), & a\lt x \leqslant b \\ y_{1}(a)=\alpha, \\ y_{2}(a)=\beta. \end{array}\right.

上述方程可以被看作是一阶向量方程的初值问题：

\left\{\begin{array}{l} \boldsymbol{Y}^{\prime}=\boldsymbol{F}(x, \boldsymbol{Y}), \quad a\lt x \leqslant b, \tag{2}\\ \boldsymbol{Y}(a)=\boldsymbol{Y}_{0} . \end{array}\right.

其中，

\boldsymbol{Y}=\binom{y_{1}}{y_{2}},\ \boldsymbol{F}(x, \boldsymbol{Y})=\binom{f\left(x, y_{1}, y_{2}\right)}{g\left(x, y_{1}, y_{2}\right)}

且

\boldsymbol{Y}_{0}=\binom{\alpha}{\beta}

于是，可以将原标量型一阶方程初值问题的数值解法（如欧拉方法、龙格-库塔方法等），直接应用于上述向量形式的方程，这样就得到了一阶方程组初值问题的数值解法。

例如，如果对式 (2) 采用欧拉方法，则其迭代公式为

\boldsymbol{Y}_{i+1}=\boldsymbol{Y}_{i}+h \boldsymbol{F}\left(x_{i}, \boldsymbol{Y}_{i}\right)

也可写成分量形式，即

\left\{\begin{array}{l} y_{1, i+1}=y_{1, i}+h f\left(x_{i}, y_{1, i}, y_{2, i}\right. \\ y_{2, i+1}=y_{2, i}+h g\left(x_{i}, y_{1, i}, y_{2, i}\right) \\ y_{1,0}=\alpha\\ y_{2,0}=\beta \end{array}\right.

同样，若对式 (2) 采用经典的四阶龙格-库塔方法，则其迭代公式为：

\boldsymbol{Y}_{i+1}=\boldsymbol{Y}_{i}+\frac{h}{6}\left(\boldsymbol{K}_{1}+2 \boldsymbol{K}_{2}+2 \boldsymbol{K}_{3}+\boldsymbol{K}_{4}\right)

其中，

\begin{align*} \boldsymbol{K}_{1}&=\boldsymbol{F}\left(x_{i}, \boldsymbol{Y}_{i}\right), \\ \boldsymbol{K}_{2}&=\boldsymbol{F}\left(x_{i}+\frac{h}{2}, \boldsymbol{Y}_{i}+\frac{h}{2} \boldsymbol{K}_{1}\right), \\ \boldsymbol{K}_{3}&=\boldsymbol{F}\left(x_{i}+\frac{h}{2}, \boldsymbol{Y}_{i}+\frac{h}{2} \boldsymbol{K}_{2}\right), \\ \boldsymbol{K}_{4}&=\boldsymbol{F}\left(x_{i}+h, \boldsymbol{Y}_{i}+h \boldsymbol{K}_{3}\right) \end{align*}

对于由个方程构成的方程组的情形也有类似的讨论，只是向量的分量个数为而已。

高阶方程初值问题

对于高阶常微分方程，可以通过引入新的变量将其转化为一阶方程组的形式。

先考察如下的二阶方程初值问题：

\left\{\begin{array}{l} y^{\prime \prime}=f\left(x, y, y^{\prime}\right), \quad a\lt x \leqslant b \tag{3}\\ y(a)=\alpha, \quad y^{\prime}(a)=\beta \end{array}\right.

通过引入新的变量，则有，就可将上述方程转化为形如式（3）的一阶方程组初值问题：

\left\{\begin{array}{l} y^{\prime}=z, \quad a\lt x \leqslant b, \\ z^{\prime}=f(x, y, z), \quad a\lt x \leqslant b, \\ y(a)=\alpha, \quad z(a)=\beta . \end{array}\right.

类似地，对于更高阶的方程，如阶方程初值问题：

\left\{\begin{array}{l} y^{(n)}=f\left(x, y, y^{\prime}, \cdots, y^{(n-1)}\right), \quad a\lt x \leqslant b \tag{4}\\ y(a)=\alpha_{1}, \quad y^{\prime}(a)=\alpha_{2}, \cdots, \quad y^{(n-1)}(a)=\alpha_{n} \end{array}\right.

按照同样的思路，可以定义，，，，，，可以将原问题转化为以下一阶方程组的初值问题：

\left\{\begin{array}{l} y_{1}^{\prime}=y_{2}, \quad a\lt x \leqslant b \tag{5}\\ y_{2}^{\prime}=y_{3}, \quad a\lt x \leqslant b \\ \cdots \\ y_{n-1}^{\prime}=y_{n}, \quad a\lt x \leqslant b, \\ y_{n}^{\prime}=f\left(x, y_{1}, y_{2}, \cdots, y_{n}\right), \quad a\lt x \leqslant b \\ y_{1}(a)=\alpha_{1}, \quad y_{2}(a)=\alpha_{2}, \cdots, \quad y_{n}(a)=\alpha_{n} \end{array}\right.

此式即为情况下的式（2）。

例如，对于三阶方程初值问题

\left\{\begin{array}{l}x y^{\prime \prime \prime}+2 x^{2} y^{\prime 2}+x^{3} y=x^{4}+1, \quad 1\lt x \leqslant 3, \\ y(1)=1, \quad y^{\prime}(1)=-2, \quad y^{\prime \prime}(1)=3 .\end{array}\right.

令，，，则可将原三阶方程问题转化为如下一阶方程组的初值问题：

\left\{\begin{array}{l} y_{1}^{\prime}=y_{2}, \quad 1\lt x \leqslant 3, \\ y_{2}^{\prime}=y_{3}, \quad 1\lt x \leqslant 3, \\ y_{3}^{\prime}=\frac{x^{4}+1-x^{3} y_{1}-2 x^{2} y_{2}{ }^{2}}{x}, \quad 1\lt x \leqslant 3, \\ y_{1}(1)=1, \quad y_{2}(1)=-2, \quad y_{3}(1)=3 . \end{array}\right.

这样，就可以采用一阶方程组初值问题的数值方法进行求解了。

例题

用经典的四阶龙格-库塔法求解一阶方程组初值问题：

\left\{\begin{array}{ll}y_{1}^{\prime}=y_{2}, & 0\lt x \leqslant 3.1, \\ y_{2}^{\prime}=-y_{1}, & 0\lt x \leqslant 3.1, \\ y_{1}(0)=1, & y_{2}(0)=0 .\end{array}\right.

取步长为 0.1。已知该方程组的精确解为，。

Fortran 实现

利用 Fortran 语言，我们可以高效地实现该问题的数值解法。为了便于比较数值算法的效果，我们同时给出原方程的精确解。具体方法参见以下的代码和注释：

program main
  implicit none
  
  integer :: N, i
  integer, parameter :: nvar = 2
  real(kind=8) :: a, b, h
  real(kind=8), dimension(nvar) :: k1, k2, k3, k4, ye, err
  real(kind=8), allocatable, dimension(:) :: x
  real(kind=8), allocatable, dimension(:,:) :: y
  
  a = 0.0D0
  b = 3.1D0
  h = 0.1D0  !步长
  N = nint((b - a)/h) !总的剖分数
  
  !动态分配长度为(N+1)的数组
  !x存放节点坐标
  !y存放对应节点的数值解
  allocate(x(N+1), y(N+1,nvar))
  
  !节点坐标
  x = [(a + i * h, i = 0, N)]
  
  y(1,:) = [1.0D0, 0.0D0] !初值
  err = 0.0D0
  ye = yexact(x(1), nvar)
  
  !打开存储文件
  open(10, file="res.txt")
  write(10,*) "#    x(1)        y(1)     yexact(1)     y(2)      yexact(2)    err_max "
  !保存初值
  write(10,"(7e12.3)") x(1), y(1,1), ye(1), y(1,2), ye(2), maxval(err)

  do i = 1, N !4阶龙格-库塔法
    x(i+1) = x(i) + h 
    k1 = h*f(x(i), y(i,:), nvar)
    k2 = h*f(x(i) + 0.5D0*h, y(i,:) + 0.5D0*k1, nvar)
    k3 = h*f(x(i) + 0.5D0*h, y(i,:) + 0.5D0*k2, nvar)
    k4 = h*f(x(i+1), y(i,:) + k3, nvar)
    y(i+1,:) = y(i,:) + (k1 + 2D0 * k2 + 2D0 * k3 + k4) / 6.0D0
 ye = yexact(x(i+1), nvar) !计算精确解
    err = abs(y(i+1,:) - ye)  !计算节点处的误差
 !保存结果
   write(10,"(7e12.3)") x(i+1), y(i+1,1), ye(1), y(i+1,2), ye(2), maxval(err)
  end do
  close(10)
  
contains
  !右端项函数
  function f(x, y, nvar) result (res)
    integer, intent(in) :: nvar  
    real(kind=8), intent(in) :: x, y(nvar)
    real(kind=8) :: res(nvar)
    res(1) = y(2)
    res(2) = -y(1) 
  end function f
  !精确解y
  function yexact(x, nvar) result (res)
    integer, intent(in) :: nvar
    real(kind=8), intent(in) :: x
    real(kind=8) :: res(nvar)
    res(1) = cos(x)
    res(2) = -sin(x) 
  end function yexact
end program main

以上程序由原来的求解标量型常微分方程的程序修改而来，将部分变量改为了数组，变动不大。感兴趣的读者可自行对比：

常微分方程的数值求解 | 龙格-库塔法

编译并运行，得到的计算结果文件res.txt的内容如下。

观察发现，当以 0.1 为步长时，计算结果大概有 5~6 位的有效数字，可见四阶龙格-库塔方法精度还是相当高的。

Gnuplot 绘图

为了更好地展示计算结果，我们使用 Gnuplot 来读取 Fortran 程序的输出文件并绘图。在 Fortran 程序同一目录下，创建并运行以下 plt 脚本文件，就可以轻松得到 png 格式的图片文件。运行结果和对应的脚本代码如下：

set terminal png
set output "res.png"

set title "Runge-Kutta Method"
set xlabel "x"
set ylabel "y"
set key right top spacing 2

plot "res.txt" using 1:2 title "Numerical y1", \
"res.txt" using 1:3 title "Analytical y1" \
with lines linewidth 2, \
"res.txt" using 1:4 title "Numerical y2", \
"res.txt" using 1:5 title "Analytical y2" \
with lines linewidth 2

小结

本文主要介绍了一阶方程组初值问题数值求解的基本原理，并通过一个具体的例子展示了如何使用 Fortran 语言来实现这一方法。对于高阶方程，本质上可以利用降阶的思想，将高阶问题转化为一阶方程组的问题，然后再进行求解。无论是简单的欧拉方法还是更为精确的龙格-库塔方法，都可以用来求解此类问题。随着计算机技术的发展，这些数值方法的应用范围也将不断扩大，为科学研究和工程技术带来极大的便利。

往期推荐

常微分方程的数值求解 | 龙格-库塔法

常微分方程的数值求解 | 从欧拉方法启航

Gnuplot：数据可视化的理想工具

常微分方程的数值求解 | 一阶方程组和高阶方程

一阶方程组初值问题

高阶方程初值问题

例题

Fortran 实现

Gnuplot 绘图

小结

推荐阅读