拉格朗日内插通用版

网上有类似的程序,但是我觉得都不完美,于是自己参考现有的代码,写了一个自己满意的。


//拉格朗日内插,n为阶数,fx,fy为整体,函数内部会对其筛选,fx,fy长度>=n
function lagrange2(fx,fy: array of Extended; xx: Extended; n: Integer): Extended;
var i,j: Integer;
    yy: Extended;
    a: array of Extended;
    fx2,fy2: array of Extended;
    sum1,st_n,half_n: Integer;
    len: Integer;
begin
  SetLength(a,n+1);
  SetLength(fx2,n+1);
  SetLength(fy2,n+1);
  sum1 := 0;
  half_n := n div 2;
  len := Length(fx)-1;   //在这里是95

//===================BEGIN=========取nn个数参与插值,使要求的数在其中间========
    st_n := Trunc(xx-fx[0]) div (15*60); //时间处于第几个15分钟
    if(st_n>=len) then    //超出所有
    begin
      for j := len-n to len do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
    end
    else if(st_n<=0) then   //低于所有
    begin
      for j := 0 to n do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
    end
    else if(len-st_n< half_n) then      //右边不够
    begin
      for j := st_n to len do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
      for j := (len-n+1) to st_n do
      begin
        fx2[sum1] := fx[j-1];
        fy2[sum1] := fy[j-1];
        sum1 := sum1+1;
      end;
    end
    else if(st_n<(n div 2)) then      //左边不够
    begin
      for j := 0 to st_n do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
      for j := (st_n+1) to (n) do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
    end
    else      //两边都够
    begin
      for j := (st_n-(n div 2)) to st_n do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
      for j := (st_n+1) to (st_n+n-(n div 2)) do
      begin
        fx2[sum1] := fx[j];
        fy2[sum1] := fy[j];
        sum1 := sum1+1;
      end;
    end;
//========================END======取nn个数参与插值,使要求的数在其中间========
  yy := 0;
  for I := 0 to n do
  begin
    a[i] := fy2[i];
    for j := 0 to n do
    begin
      if(j<>i) then
        a[i] := a[i]*(xx-fx2[j])/(fx2[i]-fx2[j]);
    end;
    yy := yy+a[i];
  end;
  Result := yy;
end;