unit1.~pas

插值计算--计算方法_delphi.rar 计算方法课程上机作业

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Grids, DBGrids, ComCtrls, ExtCtrls, Mask, Menus;
////////////////////////////////////////////////////////////////////////////////////
// y=sin(x)
const
  PN                = 100;
const
  N                 = 11;
const
  X                 : array[1..N] of double = (-5, -4, -3, -2, -1, 0,
    1, 2, 3, 4, 5);
const
  Y                 : array[1..N] of double = (0.038462, 0.058824, 0.100000, 0.200000, 0.500000, 1.000000, 0.500000, 0.200000, 0.100000, 0.058824, 0.038462);
  ////////////////////////////////////////////////////////////////////////////////////
           {
    ////////////////////////////////////////////////////////////////////////////////////
  // y=1/(1+x^2)
  //区间为：[-5,5]
  const PN=100;
  const N=11; //节点
  const X:array[1..N] of double=(-5,-4,-3,-2,-1,0,1,2,3,4,5);
  const Y:array[1..N] of double=();
  ////////////////////////////////////////////////////////////////////////////////////
      }
type
  TForm_main = class(TForm)
    StatusBar1: TStatusBar;
    Panel1: TPanel;
    Panel2: TPanel;
    PaintBox1: TPaintBox;
    Panel3: TPanel;
    GroupBox1: TGroupBox;
    RadioButton1: TRadioButton;
    RadioButton2: TRadioButton;
    GroupBox2: TGroupBox;
    StringGrid1: TStringGrid;
    Edit_x: TEdit;
    Edit_num: TEdit;
    UpDown1: TUpDown;
    Label3: TLabel;
    Label2: TLabel;
    Label1: TLabel;
    Edit_y: TEdit;
    RadioButton3: TRadioButton;
    RadioButton4: TRadioButton;
    PopupMenu1: TPopupMenu;
    N1: TMenuItem;
    N2: TMenuItem;
    Button1: TButton;
    Memo1: TMemo;
    GroupBox3: TGroupBox;
    Button2: TButton;
    Button_testpic: TButton;
    Button_reset: TButton;
    procedure Button1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure PaintBox1Paint(Sender: TObject);
    procedure PaintBox1MouseMove(Sender: TObject; Shift: TShiftState; Xp,
      Yp: Integer);
    procedure Button_resetClick(Sender: TObject);
    procedure N1Click(Sender: TObject);
    procedure N2Click(Sender: TObject);
  //  procedure Button_testClick(Sender: TObject);
    procedure Button_testpicClick(Sender: TObject);
    procedure FormResize(Sender: TObject);
  private
    function Lagrange(xi: double; dim: integer): double;
    function Newton(xi: double; dim: integer): double;
    function coordinate_pos(): integer;
    function determine_start_point(x_i: double; dimen: integer): integer;
    function Split(xi: double): double;
    function Cut_linear(xi: double): double;
    procedure PaintFrame();
  public
    paint_num: integer;                 //画图的次数，以此标志不同的颜色
    points: array[1..100, 1..3] of real;
  end;

var
  Form_main         : TForm_main;

implementation

uses Unit_func, Unit_func100;

{$R *.dfm}
//------------------------------------------------------------------------------

function TForm_main.Cut_linear(xi: double): double; //分段线性插值
var
  i, j              : integer;
  L                 : array[1..N] of double;
begin
  if (xi >= X[1]) and (xi <= X[2]) then
    L[1] := (xi - X[2]) / (X[1] - X[2])
  else
    L[1] := 0;

  for i := 2 to N - 1 do
  begin
    if (X[i - 1] <= xi) and (xi <= X[i]) then
      L[i] := (xi - X[i - 1]) / (X[i] - X[i - 1])
    else
      if (X[i] <= xi) and (xi <= X[i + 1]) then
        L[i] := (xi - X[i + 1]) / (X[i] - X[i + 1])
      else
        L[i] := 0;
  end;

  if (X[N - 1] <= xi) and (xi <= X[N]) then
    L[N] := (xi - X[N - 1]) / (X[N] - X[N - 1])
  else
    L[N] := 0;

  result := 0;
  for j := 1 to N do
    result := result + L[j] * Y[j];

  Cut_linear := result;
end;

//------------------------------------------------------------------------------

function TForm_main.Split(xi: double): double; //三次样条插值函数
var
  i, j, k           : integer;
  a_                : array[1..N] of double;
  b_                : array[1..N] of double;
  h                 : array[1..N - 1] of double;
  A                 : array[1..N] of double;
  B                 : array[1..N] of double;
  m                 : array[1..N] of double;
  s                 : double;
begin
  for i := 1 to N - 1 do
  begin
    h[i] := X[i + 1] - X[i];
  end;

  for i := 2 to N - 1 do
  begin
    a_[i] := h[i - 1] / ((h[i - 1] + h[i]));
    b_[i] := 3 * ((1 - a_[i]) * (Y[i] - Y[i - 1]) / h[i - 1] + a_[i] * (Y[i + 1]
      - Y[i]) / h[i]);
  end;

  //-----------------------以下方程假设满足第二类边界条件----------------------
  a_[1] := 1;
  a_[N] := 0;
  b_[1] := 3 * (Y[2] - Y[1]) / h[1];
  b_[N] := 3 * (Y[N] - Y[N - 1]) / h[N - 1];
  //---------------------------------------------------------------------------

  A[1] := -a_[1] / 2;
  B[1] := b_[1] / 2;

  for j := 2 to N - 1 do
  begin
    A[j] := -a_[j] / (2 + (1 - a_[j]) * A[j - 1]);
    B[j] := (b_[j] - (1 - a_[j]) * B[j - 1]) / (2 + (1 - a_[j]) * A[j - 1]);
  end;

  //追赶法解线性方程
  m[N] := (b_[N] - (1 - a_[N]) * B[N - 1]) / (2 + (1 - a_[N]) * A[N - 1]);
  for i := N - 1 downto 1 do
  begin
    m[i] := A[i] * m[i + 1] + B[i];
  end;

  i := 1;
  while xi > X[i] do
  begin
    inc(i);
  end;

  i := i - 1;

  //------构造( X[i],X[i+1] )区间上的三次样条插值函数--------------------------------------------------
  s := (1 + 2 * (xi - X[i]) / (X[i + 1] - X[i])) * ((xi - X[i + 1]) / (X[i] - X[i
    + 1])) * (xi - X[i + 1]) / (X[i] - X[i + 1]) * Y[i]
    + (1 + 2 * (xi - X[i + 1]) / (X[i] - X[i + 1])) * ((xi - X[i]) / (X[i + 1] -
    X[i])) * (xi - X[i]) / (X[i + 1] - X[i]) * Y[i + 1]
    + (xi - X[i]) * ((xi - X[i + 1]) / (X[i] - X[i + 1])) * ((xi - X[i + 1]) /
    (X[i] - X[i + 1])) * m[i]
    + (xi - X[i + 1]) * ((xi - X[i]) / (X[i + 1] - X[i])) * ((xi - X[i]) / (X[i
    + 1] - X[i])) * m[i + 1];
  Split := s;
end;
//------------------------------------------------------------------------------

function TForm_main.determine_start_point(x_i: double; dimen: integer): integer;
var
  i, start_point    : integer;
begin                                   //判断插值节点, 待优化改进
  for i := 1 to N do
  begin
    if ((X[i] <= x_i) and (X[i + 1] >= x_i)) then
    begin
      if ((dimen mod 2) = 0) then
      begin
        start_point := i - Trunc(((dimen - 2) / 2));
      end
      else
      begin                             //判断离xi最近的节点
        if (abs((x_i - X[i])) < abs((x_i - X[i + 1]))) then //离i节点更近
          start_point := i - Trunc((dimen - 2) / 2) - 1
        else                            //离i+1节点更近
          start_point := i - Trunc((dimen - 2) / 2);
      end;
      break;
    end;
  end;
  if (start_point <= 0) then
    start_point := 1;
  determine_start_point := start_point;
end;
//------------------------------------------------------------------------------

//--------------------------------NewTon插值函数--------------------------------

function TForm_main.Newton(xi: double; dim: integer): double;
var
  i_start, i_start_point, i, j, k: integer;
  op_matrix         : array of array of double; //N行,N+1列数组,x,fx,df1,df2,df3...
  sum_i, muti       : double;
begin
  i_start := self.determine_start_point(xi, dim); //自动判断插值节点
  i_start_point := i_start;
  SetLength(op_matrix, dim);            //设置行数:N
  for i := 0 to dim - 1 do
  begin
    SetLength(op_matrix[i], dim + 1);   //设置列数:N+1
    op_matrix[i, 0] := X[i_start];      //填充第一列为Xi
    op_matrix[i, 1] := Y[i_start];
    inc(i_start);                       //填充第二列为Yi
  end;
  //------------------------对N*(N+1)维矩进行运算,求各阶均差---------------
  for i := 2 to dim do                  //i表示列号
  begin
    k := 0;
    for j := i - 1 to dim - 1 do        //j表示行号
    begin
      op_matrix[j, i] := (op_matrix[j, i - 1] - op_matrix[j - 1, i - 1]) /
        (op_matrix[j, 0] - op_matrix[K, 0]); //计算均差
      inc(K);
    end;
  end;
  sum_i := 0;
  for j := 1 to dim - 1 do
  begin
    muti := 1;
    for k := 0 to j - 1 do
    begin
      muti := muti * (xi - X[i_start_point + k]);
      //X数组下标是从1开始的,第一个数是X[1]
    end;
    sum_i := sum_i + muti * op_matrix[j, j + 1];
  end;

  sum_i := sum_i + Y[i_start_point];

  Newton := sum_i;
end;
//------------------------------------------------------------------------------

//-------------------------------Lagrange插值函数-------------------------------

function TForm_main.Lagrange(xi: double; dim: integer): double;
var
  i, j, nearest_point, i_start, i_over: integer;
  fx, fxi, t        : double;
begin
  fx := 0;
  i_start := self.determine_start_point(xi, dim); //自动选择插上值节点
  i_over := i_start + dim - 1;
  for i := i_start to i_over do         //构造Lagrange插值多项式
  begin
    t := 1;
    for j := i_start to i_over do
    begin
      if (j <> i) then t := t * (xi - X[j]) / (X[i] - X[j]);
    end;
    fxi := Y[i] * t;
    fx := fx + fxi;
  end;
  Lagrange := fx;                       //返回函数值
end;

//-------------------------------显示数据---------------------------------------//

procedure TForm_main.FormCreate(Sender: TObject);
var
  i                 : integer;
begin
  paint_num := 0;
  self.StringGrid1.RowCount := N + 1;
  StringGrid1.Cells[0, 0] := 'X';
  StringGrid1.Cells[1, 0] := 'Y';
  for i := 1 to N do
  begin
    self.StringGrid1.Cells[0, i] := floattostr(X[i]);
    self.StringGrid1.Cells[1, i] := floattostr(Y[i]);
  end;
   // fun100(1.1, 2.2, 3);                //计算数据

end;
//-------------------------------------------------------------------------------

procedure TForm_main.Button1Click(Sender: TObject);
var
  xi, yi            : double;
  times             : integer;
begin
//{
  times := Trunc(strtofloat(self.Edit_num.Text));
  if self.Edit_x.Text = '' then
  begin
    MessageDlg('请输入插值点的值!', mtInformation, [mbOk], 0);
    self.Edit_x.SetFocus;
  end
  else
    if (X[1] > strtofloat(self.Edit_x.Text)) or (strtofloat(self.Edit_x.Text) >
      X[N]) then
    begin
      MessageDlg('插值节点无效,请重新输入!', mtInformation, [mbOk], 0);
      self.Edit_x.SetFocus;
      Self.Edit_x.Text := '';
    end
    else
      if times + 1 > N then
      begin
        MessageDlg('插值次数大于节点数,不能运算!', mtInformation, [mbOk], 0);
        self.Edit_x.Text := '';
        self.Edit_x.SetFocus;
      end
      else
      begin
        xi := strtofloat(self.Edit_x.text);
        if self.RadioButton1.Checked then
          yi := self.Lagrange(xi, times + 1) //调用Lagrange插值函数
        else
          if self.RadioButton2.Checked then
            yi := self.Newton(xi, times + 1) //调用Newton插值函数
          else
            if self.RadioButton3.Checked then
              yi := self.Split(xi)      //调用三次样条插值
            else
              yi := self.Cut_linear(xi); //调用分段线性插值

        self.Edit_Y.Text := floattostr(yi);
      end;
     // }
end;

//--------------------------------------------------------------------------------

function TForm_main.coordinate_pos(): integer;
var
  i, j, flag_x, flag_y: integer;
begin
  //----------------------------判断X,Y数组中是否有同号,方便轴------------------------------
  flag_x := 0;                          //flag=0时表示X或Y数组中的值均为同号, 否则flag=1表示为异号
  flag_y := 0;
  for i := 2 to N do
    if (X[1] * X[i] < 0) then
      flag_x := 1;
  if (Y[1] * Y[i] < 0) then
    flag_y := 1;

end;

//----------------------------------------------------------------------------------

procedure TForm_main.PaintFrame();
var
  i, j, x_df, y_df, xi, yi, i_height, co_margin: integer;
  yh, xw            : integer;
  sw                : integer;
  margin, x0, y0    : integer;
  dy, dx            : double;
  t, max, min, X_range, Y_range, unit_x_range, unit_y_range: double;
  x_string_width, y_string_width, y_string_height: integer;
  x_string, y_string: string;
  tempy             : real;             //对Y的值进行处理
begin
  Y_range := 0;
  self.PaintBox1.Canvas.Pen.Color := clWhite;

  margin := 50;                         //设置边距

  yh := self.PaintBox1.Height - 2 * margin;
  xw := self.PaintBox1.Width - 2 * margin;
  x0 := margin;
  y0 := margin + yh;

  //画边框
  self.PaintBox1.Canvas.MoveTo(margin, margin);
  self.PaintBox1.Canvas.LineTo(margin, yh + margin);
  self.PaintBox1.Canvas.LineTo(xw + margin, yh + margin);
  self.PaintBox1.Canvas.LineTo(xw + margin, margin);
  self.PaintBox1.Canvas.LineTo(margin, margin);

  //画坐标
  X_range := X[N] - X[1];
  max := Y[1];
  min := Y[1];
  for i := 2 to N do
  begin
    if Y[i] >= max then
      max := Y[i];
    if Y[i] <= min then
      min := Y[i];
  end;
  Y_range := 0;
  Y_range := max - min;

  x_df := N - 1;                        //X坐标等分单元格数