complexvector.~pas

以面向对象方法实现的数值算法类库

{-------------------------------------------------------------------------------
2003.08.11 定义向量类及相关运算
以基本类定义为基础，定义向量类。该类能够与传统的数据类型进行方便的转换。
--------------------------------------------------------------------------------
2003.08.20
  移动平均：开头或结尾的点，有几个就平均几个。
  求离散的差分：前差(i-(i-1))
--------------------------------------------------------------------------------
2003.08.14 实现求导、牛顿迭代
--------------------------------------------------------------------------------
2003.08.13 向量的移动:左移、右移;带参数，表示循环与否(缺省情况是循环，否则补0)。
-------------------------------------------------------------------------------}
unit ComplexVector;

interface

uses classes, MathTypes, DB;

Type

  TComplexVector = class (TPersistent)
  protected
    FData:          TComplexArray;             // 数据
    FCount:         integer;                   // 数组的大小
    FFormat:        string;                    // 显示格式
    FErrorLimit:    double;                    // 误差限
    FMaxIterations: Integer;                   // 最大迭代次数
    FRoot:          double;                    // 根
  public
    { 构造与析构方法-----------------------------------------------------------}
    constructor Create;
    destructor  Destroy; override;
    { 数据显示 ----------------------------------------------------------------}
    function    AsString: string;
    { 数据赋值 ----------------------------------------------------------------}
    procedure   Assign(Source: TPersistent); override;
    procedure   AssignArray(x: array of TComplex);
    procedure   LoadFromDatabase(ADataSet: TDataSet; ReFieldID: integer = 0;
        ImFieldID: integer = 1; Start: integer = 1; nCount: integer = 1024);
    { 数据维护 ----------------------------------------------------------------}
    function    LeftCopy(N: integer): TComplexVector;
    { 数据维护 ----------------------------------------------------------------}
    procedure   Resize(NCount: integer);
    { 常见的数值方法 ----------------------------------------------------------}
    // 平均值
    function    Average: double;
    // 多项式求值
    function    Value(x: double): double;
    // 多项式的导数
    function    Derivative: TComplexVector;
    // 牛顿迭代法求根
    function    NewtonIteration: double;
    // 移动平均
    function    MovingAverage(N: integer=3): TComplexVector;
    // 找极小值
    function    GetMinIndex: integer;
    // 找极大值
    function    GetMaxIndex: integer;
    // 找绝对极大值
    function    GetAbsMaxIndex: integer;
    // 找局部绝对极大值
    function    GetLocalAbsMaxIndex(FirstID, LastID: integer): integer;
    // 找局部极大值
    function    GetLocalMaxIndex(FirstID, LastID: integer): integer;
    // 找局部极小值
    function    GetLocalMinIndex(FirstID, LastID: integer): integer;

  published
    property Count: integer read FCount write Resize;
    property Values: TComplexArray read FData write FData;
    property DisplayFormat: string read FFormat write FFormat;
    property ErrorLimit: double read FErrorLimit write FErrorLimit;             // 误差限
    property MaxIterations: Integer read FMaxIterations write FMaxIterations;   // 最大迭代次数
    property Root: double read FRoot write FRoot;                               // 根

  end;

implementation

uses math, sysutils, dialogs;

{ TComplexVector }

procedure TComplexVector.AssignArray(x: array of TComplex);
Var
    i,N: integer;
begin
    N := Length(x);
    Resize(N);
    for i:=0 to N-1 do FData[i] := x[i];
end;

procedure TComplexVector.Assign(Source: TPersistent);
Var
    i: integer;
begin
    if Source is TComplexVector then begin
      Resize((Source as TComplexVector).FCount);
      for i:=0 to FCount-1 do FData[i] := (Source as TComplexVector).Values[i];
    end;
end;

function TComplexVector.AsString: string;
Var
    i: Integer;
begin
    Result := '';
    {
    for i:=0 to FCount-1 do
      if Result='' then
        Result := Format(FFormat,[FData[i].Re, FData[i].Im])
      else
        Result := Result + Format(' + ' + FFormat+' x^%d',[FData[i],i]);
    }
end;

function TComplexVector.Average: double;
Var
    r: double;
    i: integer;
begin
    r := 0;
    //for i:=0 to FCount-1 do r := r + FData[i];
    Result := r / FCount;
end;

constructor TComplexVector.Create;
begin
    inherited;
    FData := nil;
    FFormat := '%.4f';
    FErrorLimit := 1e-5;           // 误差限
    FMaxIterations := 1000;        // 最大迭代次数
    FRoot := 0;                    // 根
end;

// 求多项式的导数
function TComplexVector.Derivative: TComplexVector;
Var
    i: integer;
begin
    Result := TComplexVector.Create;
    //if FCount>1 then begin
    //  Result.Resize (FCount-1);     // 多项式的大小
    //  // 计算一阶导数
    //  for i:=0 to FCount-2 do Result.FData[i] := FData[i+1] * (i+1);
    //end else begin
    //  Result.Resize(1);
    //  Result.FData[0] := 0;
    //end;
end;

destructor TComplexVector.Destroy;
begin
    FData := nil;
    inherited;
end;

function TComplexVector.GetAbsMaxIndex: integer;
Var
    tmp: double;
    i,k: integer;
begin
    {
    tmp := abs(FData[0]);
    k := 0;
    for i:=1 to FCount-1 do begin
      if tmp<abs(FData[i]) then begin
        tmp := abs(FData[i]);
        k := i;
      end;
    end;
    }
    Result := k;
end;

function TComplexVector.GetLocalAbsMaxIndex(FirstID, LastID: integer): integer;
Var
    tmp: double;
    i,k: integer;
begin
    {
    tmp := abs(FData[FirstID]);
    k := FirstID;
    for i:=FirstID+1 to LastID do begin
      if tmp<abs(FData[i]) then begin
        tmp := abs(FData[i]);
        k := i;
      end;
    end;
    }
    Result := k;
end;

function TComplexVector.GetLocalMaxIndex(FirstID, LastID: integer): integer;
Var
    tmp: double;
    i,k: integer;
begin
    {
    tmp := FData[FirstID];
    k := FirstID;
    for i:=FirstID+1 to LastID do begin
      if tmp<FData[i] then begin
        tmp := FData[i];
        k := i;
      end;
    end;
    }
    Result := k;
end;

function TComplexVector.GetLocalMinIndex(FirstID, LastID: integer): integer;
Var
    tmp: double;
    i,k: integer;
begin
    {
    tmp := FData[FirstID];
    k := FirstID;
    for i:=FirstID+1 to LastID do begin
      if tmp>FData[i] then begin
        tmp := FData[i];
        k := i;
      end;
    end;
    }
    Result := k;
end;

function TComplexVector.GetMaxIndex: integer;
Var
    tmp: double;
    i,k: integer;
begin
    {
    tmp := FData[0];
    k := 0;
    for i:=1 to FCount-1 do begin
      if tmp<FData[i] then begin
        tmp := FData[i];
        k := i;
      end;
    end;
    }
    Result := k;
end;

function TComplexVector.GetMinIndex: integer;
Var
    tmp: double;
    i,k: integer;
begin
    {
    tmp := FData[0];
    k := 0;
    for i:=1 to FCount-1 do begin
      if tmp>FData[i] then begin
        tmp := FData[i];
        k := i;
      end;
    end;
    }
    Result := k;
end;

function TComplexVector.LeftCopy(N: integer): TComplexVector;
Var
    i: Integer;
begin
    Result := TComplexVector.Create;
    Result.Resize(N);
    for i:=0 to N-1 do Result.Values[i] := FData[i];
end;

procedure TComplexVector.LoadFromDatabase(ADataSet: TDataSet; ReFieldID: integer = 0;
        ImFieldID: integer = 1; Start: integer = 1; nCount: integer = 1024);
Var
    i,N: integer;
begin
    N := ADataSet.RecordCount;
    Resize(nCount);
    ADataSet.FindFirst;
    for i:=0 to Start-1 do ADataSet.Next;
    if nCount<0 then begin
      for i:=0 to N-1 do
        if NOT ADataSet.Eof then begin
          FData[i].Re := ADataSet.Fields[ReFieldId].AsFloat;
          FData[i].Im := ADataSet.Fields[ImFieldId].AsFloat;
        end  else break;
    end else begin
      for i:=0 to nCount-1 do
        if NOT ADataSet.Eof then begin
          FData[i].Re := ADataSet.Fields[ReFieldId].AsFloat;
          FData[i].Im := ADataSet.Fields[ImFieldId].AsFloat;
        end else break;
    end;
end;

function TComplexVector.MovingAverage(N: integer): TComplexVector;
Var
    L, R, i, Z, j: integer;
    y: TComplexVector;
    a: double;
begin
    if not odd(N) then N := N+1;
    Z := N div 2;
    y := TComplexVector.Create;
    y.Resize(FCount);
    for i:=0 to FCount-1 do begin
      a := 0;
      // 左端
      L := 0;
      if i>=Z then begin
        for j:=i-Z to i-1 do begin
          a := a + FData[j];
          inc(L);
        end;
      end else begin
        for j:=0 to i-1 do begin
          a := a + FData[j];
          inc(L);
        end;
      end;
      // 中心点
      a := a + FData[i];
      // 右端
      R := 0;
      if i<FCount-Z then begin
        for j:=i+1 to i+Z do begin
          a := a + FData[j];
          inc(R);
        end;
      end else begin
        for j:=i+1 to FCount-1 do begin
          a := a + FData[j];
          inc(R);
        end;
      end;
      y.FData[i] := a / (L + R + 1);
    end;
    Result := y;
end;

function TComplexVector.NewtonIteration: double;
Var
    i: Integer;
    xn, xn_1, e: double;
    dp: TComplexVector;
begin
    // 计算一阶导数
    dp := Derivative;
    // 开始迭代
    i := 0;
    xn_1 := FRoot;
    Repeat
      xn := xn_1 - poly(xn_1, FData) / poly(xn_1, dp.FData);
      e := abs(xn - xn_1);
      inc(i);
      if e >= FErrorLimit then xn_1 := xn;
    Until (i>FMaxIterations) or (e < FErrorLimit);
    if i>FMaxIterations then showmessage('牛顿迭代不收敛！');
    FRoot := (xn + xn_1) / 2;
    Result := FRoot;
end;

procedure TComplexVector.Resize(NCount: integer);
begin
    if NCount<>FCount then begin
      FData := nil;
      SetLength(FData, NCount);
      FCount := NCount;
    end;
end;

function TComplexVector.Value(x: double): double;
begin
    Result := Poly(x, FData);
end;


end.