常微分方程两点边值问题的差分求解 | Dirichlet 边界条件

两点Dirichlet边值问题

对于一个二阶常微分方程的 Dirichlet 边值问题，其一般形式可以表示为：

这里表示关于的二阶导数，是已知函数，和分别是在区间端点和处给定的边界值。

为了讨论方便，这里只考虑如下的两点 Dirichlet 边值问题：

其中，和为已知函数。

差分法的基本原理

差分法是一种数值方法，它通过将连续的微分方程离散化为一系列代数方程组，从而使得原本难以直接求解的问题变得易于处理。

首先对求解区域等距剖分为份，得到个节点，，步长。然后将原方程弱化，使之仅在内部节点成立，即有

继续利用中心差分格式对二阶导数值进行精细处理，可得在节点处的离散方程：

其中，，。忽略上式中的高阶小项，并用数值解作为精确解的近似值，并结合边界条件，可得以下差分格式：

这是一个由个未知量、个方程构成的线性方程组，写成矩阵形式即为

其中

注意到为对称的三对角矩阵，因此，可以使用追赶法、高斯消元法或其他适当的算法求解该方程组，从而获得每个内部节点上的值。

例题

用差分法求解 Dirichlet 边值问题：

将计算区域剖分成 16 份，即取步长，计算出节点处的数值解，并给出与精确解的误差。

Fortran 实现

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

program main
  implicit none
  integer:: m, i
  real(8),parameter::Pi=3.14159265359d0
  real(8):: a, b, h, alpha, beta
  real(8),allocatable:: matrix_a(:), matrix_b(:), rhs(:), ans(:)
  real(8),allocatable:: x(:), y(:)
  m =16!总的剖分数
  a =0.0d0!边界左端点
  b =Pi/2.0d0!边界右端点
  h =(b - a)/dble(m)!步长
  alpha =0.0d0!左端点值
  beta =1.0d0!右端点值

! 节点坐标
  allocate(x(0:m))
  x =[(a + i * h, i =0, m)]

! 矩阵的右端项，含m-1个元素的数组rhs
  allocate(rhs(m-1))
  rhs =[(h*h*f(x(i)), i =1, m-1)]
! 考虑边界条件
  rhs(1)= rhs(1)- alpha
  rhs(m-1)= rhs(m-1)- beta

! 对称三对角矩阵，主、次对角线元素数组a、b，计算结果存储数组ans
  allocate(matrix_a(m-1), matrix_b(m-1), ans(m-1))
  matrix_a =1.0d0
  matrix_b =[(h*h*q(x(i))-2.0d0, i =1, m-1)]

! 方程组求解
  call chase_algorithm(matrix_a, matrix_b, matrix_a, m-1, rhs, ans)
  deallocate(matrix_a, matrix_b, rhs)

! y为数值解，组装计算结果
  allocate(y(0:m))
  y(0)= alpha
  y(1:m-1)= ans
  y(m)= beta
  deallocate(ans)

! 保存结果到文件
  open(10,file="res.txt")
  write(10,*)"#     x          y         yexact       err"
  do i =0, m
    write(10,"(4e12.5)") x(i), y(i), exact(x(i)),abs(exact(x(i))- y(i))
  enddo
  close(10)

contains
!原方程右端项函数f(x)
  function f(x) result(res)
  real(8):: res, x
    res =-(x*x +1)*sin(x)
  end function f
!原方程中的函数q(x)
  function q(x) result(res)
  real(8):: res, x
    res =-x*x
  end function q
!精确解
  function exact(x) result(res)
  real(8):: res, x
    res =sin(x)
  end function exact
!解决三对角方程组的追赶法程序
!a,b,c分别为系数矩阵的下次、主、上次对角线元素组
!n为方程组阶数，d为矩阵右端项向量
!ans为计算结果
  subroutine chase_algorithm(a, b, c, n, d, ans)
  integer,intent(in):: n
  real(8),intent(in):: a(n), b(n), c(n), d(n)
  real(8),intent(out):: ans(n)
  real(8),allocatable:: g(:), w(:)
  integer:: i
  real(8):: p

  allocate(g(n), w(n))

  g(1)= d(1)/ b(1)
  w(1)= c(1)/ b(1)

  do i =2, n
      p = b(i)- a(i)* w(i-1)
      g(i)=(d(i)- a(i)* g(i-1))/ p
      w(i)= c(i)/ p
  enddo
  ans(n)= g(n)
  do i = n-1,1,-1
      ans(i)= g(i)- w(i)* ans(i+1)
  enddo

  deallocate(g, w)
  end subroutine chase_algorithm
end program main

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

Gnuplot 绘图

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

set terminal png
set output "res.png"

set title "Finite Difference Method"
set xlabel "x"
set ylabel "y"
set key top left spacing 2

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

结语

本文介绍了如何利用差分法解决常微分方程的两点 Dirichlet 边值问题，并通过一个具体例子展示了从理论到实践的全过程。差分法作为一种基本而有效的数值方法，在处理边界值问题时展现了强大的适用性和灵活性。