第四十九篇 二阶常微分方程的打靶法

边界值问题

当我们试图用自变量的不同值所提供的信息来解二阶或更高阶的常微分方程时，我们必须使用与前面描述的不同的数值方法。
前面提到的初值问题经常涉及到“时间”作为自变量，解决技术要求我们按步骤“前进”，直到达到所需的解。在这种情况下，解的定义域不是有限的，因为原则上我们可以无限地沿着正方向或负方向前进。
边值问题涉及一个有限解的区间，在这个区间内去求解。这类问题的自变量通常是空间中测量距离的坐标。一个典型的二阶边值问题可能下面的形式。

解的定义域(假设为B >a)是由x在A≤x≤B的范围内的值给出的，我们需要的是找到对应于x在这个范围内的y的值。
本篇将考虑的数值解法分为三类
(a)用有限差分等效代替原来的微分方程的技术。
(b)“打靶法”，试图用一个等价初值问题代替边值问题。
“加权残差”法，即猜测一个满足边界条件的试解，并调整解内的某些参数以使误差最小化。

有限差分法

在这种方法中，微分方程中的导数项被有限差分近似代替。正如微分方程是线性的，将会显示的那样。这个过程会导致一个线性联立方程，可以使用前面描述的方法求解。
首先，需要定义一些有限差分近似来处理经常遇到的导数。考虑下图的解曲线，其中x轴被细分为规则的网格点，距离为h。

可以用下面任意一种方法来表示y对x的一阶导数

这些有限差分公式通过计算函数在xi附近的连接点的直线的斜率来逼近一阶导数。正向和反向公式有偏差，因为它们只包括xi的右侧或左侧的点。中心差分形式将xi左右的点考虑在内。
高阶导数也可以用这种方法近似。例如，在xi处的二阶导数可以通过取一阶导数的逆向差来估计，如下所示

如果我们把一阶导数的前向差分公式代入，得到

类似上面可以得到yi’‘的中心差分公式

类似地，三阶导数的中心差分公式可以由

对x的四阶导数

这些高阶导数的正向和逆向差分也很容易得到，并且通过在差分公式中包含更多的点，可以获得更高的精度。
逆向差分公式是正向差分公式的简单镜像，除了奇数导数(即，y’，y " '等)，其中符号必须相反。
一阶导数在x0处的四点正向差分公式为

逆向差分公式为

“两点边值问题的解决方案包括把x的范围分裂到n个相等的部分,每个宽度h。n = 4,如果给定的边界条件,x = A和x = B,令ξ= x0 + ih, i = 1, 2,……，其中x0 = A, xn = B。
在每个点i = 1,2，…n−1处，微分方程以有限差分形式表示。如果微分方程是线性的，这将导致n−1个联立线性方程，在未知值yi, i = 1,2，…，n−1，其中yi表示在xi处的解。细分得越多，解的细节和精度就越高，但必须解决的联立方程也就越多。

计算实例

给定y’’ = 3x + 4y，边界条件y(0) = 0, y(1) = 1，利用h = 0.2有限差分求解0≤x≤1范围内的方程。

首先将微分方程写成有限差分形式。选择二阶导数公式。通常使用最低阶中心差分形式，除非有特殊的情况去使用较不准确的正向或逆向差分形式，因此

微分方程可以写成

求解域被分割成5条相等的条带，每条宽度h = 0.2。
然后将有限差分方程写在每一个需要解的网格点上。需要时引入已知的边界条件，得到下表所示的四个方程。
这些(非对称)方程可以用线性方程求解中的任何合适的程序求解。第二个表将数值解与精确解进行了比较




显然，有限差分解的精度可以通过在0≤x≤1范围内进行更多的细分来进一步提高，但代价是求解更多的方程。
上个例子显示了一个问题，那就是如果使用了y "的高阶有限差分表示。例如，如果对y’‘使用五点中心差分公式，那么为了表示y1’‘和y4’’，就需要y在解域之外的值。为了解这个系统，还需要更多的信息，因为未知数要比方程多。
在解域之外引入点(至少是暂时引入点)的要求经常会碰到，并且可以通过采取合并适当的边界条件来解决，如下一个示例所示。
考虑一个边值问题

在这种情况下,y (B)是未知的,所以微分方程的有限差分需要写在x = B,即使使用最简单的中心差分公式y”，相应的“解决方案”在x5,将如下图所示。在这种情况下，我们有5个未知数，但只有4个方程。

第五个方程来自导数边界条件，它也可以写成有限差分形式，即使用中心差分

计算实例

给定x2y’’ + 2xy’−2y = x2，边界条件y(1) = 1, y’(1.5) =−1。用h = 0.1有限差分求解1≤x≤1.5范围内的方程。
微分方程以有限差分形式，对两个导数项使用中心差分写成如下

并应用于5个需要解的x值，如下图所示。
得到的方程如下表所示。从下图可以看出，第五个方程在解域外引入了第六个未知数y6 = y(1.6)。中心差分形式的导数边界条件



用于提供第六个方程。
这六个线性方程的解，连同精确解，给出在下表中


本节给出的有限差分解在网格点相对较少的情况下表现得相当好。但并不会总是这种情况，特别是当导数在解域中发生突然的变化时。在没有精确解可供比较的情况下，建议使用两个或三个不同层次的解域来解决问题。解对网格尺寸参数h的灵敏度往往表明解的精度。
从我们的角度来看，应该使h足够小，以达到所需的精度，但不能小于必要的值，因为这会导致过度大的方程组。

打靶法

打靶法试图像解决初值问题一样解决边值问题。考虑下面的二阶方程

在一个打靶法中，将解决一系列的初值问题的形式

通过有条不紊地改变初始梯度ai，我们最终可以在x = B处得到一个足够接近所需边值yB的解。有几种不同的策略来收敛所需的解，这些方法类似于寻找非线性代数方程的根。
这里描述的方法如下图所示是选择两个初始梯度y ’ 0 (A) = a0和y ’ 1 (A) = a1,通过一步法之后，得到y0(B) = y0和y1 (B) = y₁分别在所需的边界条件yB的两边,即,y0 < yB < y1(或y1< yB < y0)。

通过线性插值，给出了初始梯度的改进估计

从而得出y(B) = y∗。
根据下面的测试，其中一个初始梯度被∗取代，如果

然后用y代替y0和a代替a0
或者用y代替y1和a代替a1
在整个迭代过程中，我们的目标是保持一个高估目标边界条件的初始梯度，以及一个低估目标边界条件的初始梯度。随着迭代的进行，y∗趋于目标值，当根据准则满足收敛容差时，计算停止

迭代过程本质上就是之前描述的“假位置法”。详情可见试位法求解非线性方程

程序如下

程序计算的二阶常微分方程为

import numpy as np
a0=np.zeros((2))
k0=np.zeros((2))
k1=np.zeros((2))
k2=np.zeros((2))
k3=np.zeros((2))
y=np.zeros((2))
nsteps=5;xa=0.0;ya=0.0;xb=1.0;yb=1.0
a0[:]=(0.0,1.0)
tol=0.0001;limit=25
y0=np.zeros((nsteps+1,2))
ystar=np.zeros((nsteps+1))
print('--二阶常微分方程的打靶法--')
h=(xb-xa)/nsteps
def f74(x,y):
    f74=np.zeros((2))
    f74[0]=y[1]
    f74[1]=3.0*x**2+4.0*y[0]
    return f74
for j in range(1,3):
    x=xa;y[0]=ya;y[1]=a0[j-1]
    for i in range(0,nsteps+1):
        y0[i,j-1]=y[0]
        k0=h*f74(x,y);k1=h*f74(x+h/2.0,y+k0/2.0)
        k2=h*f74(x+h/2.0,y+k1/2.0);k3=h*f74(x+h,y+k2)
        y=y+(k0+2.0*k1+2.0*k2+k3)/6.0;x=x+h
if (y0[nsteps,0]-yb)*(y0[nsteps,1]-yb)>0:
    print('尝试一个新梯度')
iters=0
while(True):
    iters=iters+1
    astar=a0[0]+(yb-y0[nsteps,0])*(a0[1]-a0[0])/(y0[nsteps,1]-y0[nsteps,0])
    x=xa;y[0]=ya;y[1]=astar
    for i in range(0,nsteps+1):
        ystar[i]=y[0]
        k0=h*f74(x,y);k1=h*f74(x+h/2.0,y+k0/2.0)
        k2=h*f74(x+h/2.0,y+k1/2.0);k3=h*f74(x+h,y+k2)
        y=y+(k0+2.0*k1+2.0*k2+k3)/6.0;x=x+h
    if (abs(ystar[nsteps]-yb)/yb)<tol:
        print(' x              y')
        for i in range(0,nsteps+1):
            print('{:9.5e}'.format(xa+i*h),end='  ')
            print('{:9.5e}'.format(ystar[i]))
        print('迭代到收敛次数',iters)
        break
    if (ystar[nsteps]-yb)*(y0[nsteps,0]-yb)>0:
        y0[:,0]=ystar;a0[0]=astar
    else:
        y0[:,1]=ystar;a0[1]=astar

        

终端输出结果如下