预测-修正法(Milne-Simpson和Adams-Bashforth-Moulton)解常微分方程(python)

第四十八篇 预测-修正法

预测-校正方法是使用先前几个已知点的信息来计算下一个点，如下图所示。
这种方法的一个缺点是，他们需要一些已知信息，可能需要一步法来生成一些点才能开始，详情可见常微分方程的一步法求解。而这些方法的优点在于，为了预测下一个点，能够有效地利用现有的信息。这与四阶龙格-库塔法相反，例如，在每一步中，需要计算四个函数，预测-修正法将不会这么麻烦。
预测-校正方法利用两个公式;根据现有数据以估计下一个点的预测公式，以及改进此估计的校正公式。

可以反复应用校正公式，直到满足某一收敛准则;但是，本篇中描述的任何方法都没有实现此选项。
预测公式通过在y ‘- x曲线下使用xi, xi−1,xi−2等样本点积分来估计yi+1的新值。任何不需要对yi+1进行先验估计的数值积分公式都适合作为预测器使用。
修正公式改进了预测值yi+1，再次在y ‘-x曲线下积分，但这次使用了样本点xi+1, xi, xi−1等。校正公式能够在xi+1处取样，因为从预测阶段可以得到y’i+1的值。任何需要对y’i+1进行先验估计的数值积分公式都适合用作校正器。
之前描述的修正欧拉方法是一种预测-校正器。给出一个标准形式的一阶微分方程，y’ = f(x, y)， y(x0) = y0，该方法从Euler方法开始，这是一个使用矩形规则的预测器

接下来是一个使用梯形法则的校正器

注意，预测器不需要预估y’i+1，而校正器需要。
然而，最著名的预测-校正方法使用的公式在算法的两个部分具有相同的精度。如上所示，这是一种简便的方法可以根据预测项和修正项之间的差异来估计修正项的误差。

米尔恩-辛普森方法

这种方法使用米尔恩的公式作为预测器，用辛普森法则作为校正器。该方法是四阶的，即预测器和校正器的主要误差项都包含h5，并且需要4个初始值y来开始。注意，与预测量相比，与校正量公式相关的误差项较小，这是意料之中的，因为预测量涉及的外推过程不那么精确。
给出一个具有四个初始条件y(xi−3)= yi−3,y(xi−2)= yi−2,y(xi−1)= yi−1和y(xi) = yi的一阶微分方程y’ = f(x, y)，该方法从米尔恩预测器开始

接着是辛普森的校正者

米尔恩预测器使用下图所示积分的三个样本点，在xi−3和xi+1之间对y’- x曲线进行积分。

计算实例

已知y ’ = 2x2 + 2y，有初始条件

使用米尔恩-辛普森预测-校正方法估计y(0.2)
给出了步长h = 0.2的四个初始条件。为了便于手算，可以用计数器制作一个初始条件表格，计算右侧导数函数f(x, y)在每一点的值。

首先应用预测器。

这使得在xi+1 = 0.2处的导数可以预测为fi+1(0.2, 1.4945) = 3.0689，并应用校正器。

这种情况下精确到小数点后四位的的精确解是y(0.2) = 1.4977。
使用米尔恩的方法，或者任何使用辛普森法则作为校正器的方法的一个不利是，在计算的一个阶段产生的误差随后可能会增大。由于这个原因，其他的四阶方法，比如下面描述的方法，往往是更受欢迎的。

Adams-Bashforth-Moulton方法

一种更稳定的四阶方法是基于Adams- bashforth预测器和Adams Moulton校正器，在这种方法中误差不会快速增长。
给定一个一阶微分方程，y’ = f(x, y)具有四个初始条件y(xi−3)= yi−3,y(xi−2)= yi−2,y(xi−1)= yi−1,y(xi) = yi = yi，该方法从Adams-Bashforth的预测器开始

接下来是亚当斯-莫尔顿的校正器

adam - bashforth - moulton方法有比Milne Simpson更大的误差项，尽管主要的误差项仍然表明校正器比预测器更准确。
稳定性的取得，要以一些额外的工作为代价，因为这两个公式需要四个样本点，而不是Milne-Simpson方法中的三个。如下图所示，Adams-Bashforth Predictor使用4个样本点在y’ -x曲线下xi和xi+1的界限之间进行积分。亚当斯-莫尔顿校正器是类似的，但样本点向右移动了一个。

计算实例

已知y’ = 2x2 + 2y，有初始条件

使用Adams-Bashforth-Moulton预测-校正器方法估计y(0.2)。
给出了步长h = 0.2的四个初始条件。为了便于手算，可以用计数器制作一个初始条件表格，计算右侧导数函数f(x, y)在每一点的值。

首先应用预测器公式

这使得在xi+1 = 0.2处的导数可以预测为fi+1(0.2, 1.4941) = 3.0682，并应用校正器公式。

这种情况下精确到小数点后四位的精确解是y(0.2) = 1.4977。

程序如下

#线性常微分方程的theta法
import numpy as np
itype=2;nsteps=5;h=-0.05
x=np.zeros((5))
y=np.zeros((5))
x[0:4]=(1.00,0.95,0.90,0.85);y[0:4]=(3.61623,2.99272,2.55325,2.22755)
def f73(x,y):f73=x*y**2+2.0*x**2return f73 
if itype==1:print('**Milne-Simpson 4阶P-C法**')print('x          y    Error')for i in range(1,5):print('{:9.5e}'.format(x[i-1]),end='  ')print('{:9.5e}'.format(y[i-1]))for j in range(0,nsteps+1):x[4]=x[3]+hy4=y[0]+4.0*h/3.0*(2.0*f73(x[1],y[1])-f73(x[2],y[2])+2.0*f73(x[3],y[3]))y[4]=y[2]+h/3.0*(f73(x[2],y[2])+4.0*f73(x[3],y[3])+f73(x[4],y1))e=-(y[4]-y4)/29print(x[4],y[4],e)y[0:4]=y[1:5];x[0:4]=x[1:5]
elif itype==2:print('**Adams-Bashforth-Moulton4阶P-C法**')print('    x           y          Error')for i in range(1,5):print('{:9.5e}'.format(x[i-1]),end='  ')print('{:9.5e}'.format(y[i-1]))for j in range(0,nsteps+1):x[4]=x[3]+hy4=y[3]+h/24.0*(-9.0*f73(x[0],y[0])+37.0*f73(x[1],y[1])-59.0*f73(x[2],y[2])+55.0*f73(x[3],y[3]))y[4]=y[3]+h/24.0*(f73(x[1],y[1])-5.0*f73(x[2],y[2])+19.0*f73(x[3],y[3])+9.0*f73(x[4],y4))e=-(y[4]-y4)/14print('{:9.5e}'.format(x[4]),end='  ')print('{:9.5e}'.format(y[4]),end='  ')print('{:9.5e}'.format(e))y[0:4]=y[1:5];x[0:4]=x[1:5]

终端输出结果如下