milne.m

来自「基于Matlab的计算方法中常用的算法程序」· M 代码 · 共 29 行

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