小弟写的一个拉各朗日三点插值程序（源码），无法控制递减序列，请各位大虾帮忙看看

Question

procedure TForm1.three(num:Integer; a1: Double; var b1: Double; x,y: TDigits);
var
 i,j,k,m:Integer;
 z,s:Double;
begin
 z:=0.0;
 if num<1 then b1:=z
  else if num=1 then b1:=y[0]
   else if num=2 then b1:=(y[0]*(a1-x[1])-y[1]*(a1-x[0]))/(x[0]-x[1])
   else
   begin
   //选三点
    if a1<=x[1] then
     begin
      k:=0;
      m:=2
     end
    else if a1>=x[num-2] then
     begin
      k:=num-3;
      m:=num-1;
     end
     else
       begin
        k:=0;
        m:=num-1;        
        while (m-k)<>1 do
         begin
           i:=trunc((m+k)/2);           
           if a1<x[i-1] then m:=i
           else k:=i;
         end;
           k:=k-1;
           m:=m-1;
           if abs(a1-x[k-1])<abs(a1-x[m+1]) then k:=k-1
           else m:=m+1;
       end;
   end;
   //插值计算
   z:=0.0;
   for i:=k to m do
    begin
     s:=1.0;
     for j:=k to m do
      begin
       if (j<>i) then
        begin
         s:=s*(a1-x[j])/(x[i]-x[j])
        end;
      end;
      z:=z+s*y[i];
    end;
    b1:=z;
end;

Answer 1

从盒子里面摘录一个给你，希望能对你有所帮助procedure POLINT(XA, YA:array of real; N:integer; X:real; var Y, DY:real);
var
    C,D:array[0..10] of real;
    DIF,DIFT,HO,HP,W,DEN:real;
    NS,I,M:integer;
begin
    NS:=1;
    DIF:=Abs(X - XA[1]);
    For I:=1 To N do
    begin
        DIFT:=Abs(X - XA[I]);
        If DIFT < DIF Then
        begin
            NS:=I;
            DIF:=DIFT;
        end;
        C[I]:=YA[I];
        D[I]:=YA[I];
    end;
    Y:=YA[NS];
    NS:=NS - 1;
    For M:=1 To N - 1 do
    begin
        For I:=1 To N - M do
        begin
            HO:=XA[I] - X;
            HP:=XA[I + M] - X;
            W:=C[I + 1] - D[I];
            DEN:=HO - HP;
            If DEN = 0 Then
            begin
                ShowMessage('PAUSE');
                Exit;
            end;
            DEN:=W / DEN;
            D[I]:=HP * DEN;
            C[I]:=HO * DEN;
        end;
        If 2 * NS < N - M Then
            DY:=C[NS + 1]
        Else
        begin
            DY:=D[NS];
            NS:=NS - 1;
        end;
        Y:=Y + DY;
    end;
end;

Answer 2

谢谢楼上的兄弟，我自己找了一下原因，发现问题了，解决了。
主要原因是原程序在处理递减序列时，选三点错误。
现在已解决。谢谢。
procedure Tfrm_dmjs.three(num:Integer; a1: Double; var b1: Double; x,y: TDigits);
var
 i,j,k,m:Integer;
 z,s:Double;
 x_pie,y_pie:array of Double;
begin
 {-------先处理序列-------}
 SetLength(x_pie,num);
 SetLength(y_pie,num);
 if x[0]>x[num-1] then
  begin
   for i:=0 to num-1 do
    begin
     x_pie[i]:=x[num-1-i];
     y_pie[i]:=y[num-1-i];
    end;
  for i:=0 to num-1 do
   begin
    x[i]:=x_pie[i];
    y[i]:=y_pie[i];
   end;
  end; 
 z:=0.0;
 if num<1 then b1:=z
  else if num=1 then b1:=y[0]
   else if num=2 then b1:=(y[0]*(a1-x[1])-y[1]*(a1-x[0]))/(x[0]-x[1])
   else
   begin
   //选三点
    if a1<=x[1] then
     begin
      k:=0;
      m:=2
     end
    else if a1>=x[num-2] then
     begin
      k:=num-3;
      m:=num-1;
     end
     else
       begin
        k:=0;
        m:=num-1;
        while (m-k)<>1 do
         begin
           i:=trunc((m+k)/2);           
           if a1<x[i-1] then m:=i
           else k:=i;
         end;
           k:=k-1;
           m:=m-1;
           if abs(a1-x[k-1])<abs(a1-x[m+1]) then k:=k-1
           else m:=m+1;
       end;
   end;
   //ShowMessage('k'+IntToStr(k));
  // ShowMessage('m'+IntToStr(m));
   //插值计算
   z:=0.0;
   for i:=k to m do
    begin
     s:=1.0;
     for j:=k to m do
      begin
       if (j<>i) then
        begin
         //ShowMessage('j'+IntToStr(j));
         //ShowMessage('i'+IntToStr(i));
         s:=s*(a1-x[j])/(x[i]-x[j]);
        end;
      end;
      z:=z+s*y[i];
    end;
    b1:=z;
end;