线性多步法

在逐步推进的求解过程中，在计算yi+1以前，事实上已经求出了一系列的近似值y0,y1,…,yi。如果能充分利用第i+1步前面的多步信息来预测yi+1，那么可能会获得比较高的精度，这就是构造线性多步法的基本思想。

对于方程y=f(x,y)，它的解可表示为

由于对\(\int_{{x}_{i}}^{{x}_{i+1}} {f(x,y(x))dx}\)用矩形法和梯形法进行数值积分，分别获得了求微分方程数值解的Euler公式和梯形公式。因此，若需要再提高精度，可用更精确的求积公式来计算，也就是对被积函数用更高次的插值多项式来代替，选取不同的插值公式，就会得到不同的数值解法，比如基于数值积分的Milne公式、Simpson公式、Admas公式等。这里选用Milne公式进行求解。Milne公式的形式如下所示：

其中，yi表示y(xi)的近似值fi=f(xi,yi)。另外，前面几步的yi可先用别的单步法（如改进的Euler法、Runge-Kutta法）计算出来，然后再利用上述公式进行求解。在下面编写的M文件Milne.m中，首先用改进的Euler法计算出前几步的值，然后用上面的Milne公式进行求解。

function s=Milne(fun,x0,xn,y0,n)
  % 用线性多步法计算常微分方程初值问题
  %x0为初值条件y(x0)=y0中的x0，y0为初值条件中的y0
  %xn为x取值区间的最后一个节点的横坐标值
  %n为把区间分成的等份数
  if nargin<5
      error
      return
  end
  h=(xn-x0)/n;
  y(1)=y0;
  xx=x0;
  % 先用改进的Euler方法计算前4个y的初值
  for i=1:3
      yp=y(i)+h*feval(fun,xx,y(i));
      xx=x0+i*h;
      yc=y(i)+h*feval(fun,xx,yp);
      y(i+1)=(yp+yc)/2;
  end
  % 然后利用Milne公式计算最后的解
  for i=1:n
      xx=x0+i*h;
      t1=2*feval(fun,xx,y(i+1));
      t2=-feval(fun,(xx+h),y(i+2));
      t3=2*feval(fun,(xx+2*h),y(i+3));
      y(i+4)=y(i)+4/3*h*(t1+t2+t3);
  end
  s=y(n);
  return

【例29】求解下列初值问题。

\(\left\{ \begin{aligned} \&\frac{dy}{dx}=-2xy \\ \&y{|}_{x=0}=1 \end{aligned} \right.\) ， \(0\leq x\leq1.2\)

在MATLAB中实现时，按照下面的步骤进行。

首先，编写上述函数的M文件ff3.m。

function y=ff3(x,y)
  y=-2*x*y;

然后，可以用编写的Milne线性多步法的M函数进行求解。

>> Milne('ff3',0,1.2,1,100)
ans =
    0.2438