unit2.pas

《Delphi常用数值算法集》的配书源码

unit Unit2;

interface
uses
  unit1;

procedure CHEBFT(A, B:real; var C:array of real; N:integer);
procedure CHEBPC(C:array of real;var D:array of real; N:integer);
procedure PCSHFT(A, B:real;var D:array of real; N:integer);

implementation
procedure CHEBFT(A, B:real; var C:array of real; N:integer);
var
     K,J:integer;  SUM,Y,BMA,BPA,FAC:real;
    F:array[0..50] of real;
const
    NMAX = 50;    PI = 3.14159265358979;
begin
    BMA:=0.5 * (B - A);
    BPA:=0.5 * (B + A);
    For K:=1 To N do
    begin
        Y:=Cos(PI * (K - 0.5) / N);
        F[K]:=FUNC(Y * BMA + BPA);
    end;
    FAC:=2/ N;
    For J:=1 To N do
    begin
        Sum:=0 ;
        For K:=1 To N do
        begin
            Sum:=Sum + F[K] * Cos((PI * (J - 1)) * ((K - 0.5) / N));
        end;
        C[J]:=FAC * Sum;
    end;
end;

procedure CHEBPC(C:array of real;var D:array of real; N:integer);
var
    DD:array[0..50] of real;  J,K:integer;
    SV:real;
begin
    For J:=1 To N do
    begin
        D[J]:=0 ;
        DD[J]:=0; 
    end;
    D[1]:=C[N];
    For J:=N - 1 DownTo 2 do
    begin
        For K:=N - J + 1 DownTo 2 do
        begin
            SV:=D[K];
            D[K]:=2* D[K - 1] - DD[K];
            DD[K]:=SV;
        end;
        SV:=D[1];
        D[1]:=-DD[1] + C[J];
        DD[1]:=SV;
    end;
    For J:=N DownTo 2 do
        D[J]:=D[J - 1] - DD[J];
    D[1]:=-DD[1] + 0.5 * C[1];
end;

procedure PCSHFT(A, B:real;var D:array of real; N:integer);
var
    FAC,CONST1:real;  J,K:integer;
begin
    CONST1:=2/ (B - A);
    FAC:=CONST1;
    For J:=2 To N do
    begin
        D[J]:=D[J] * FAC;
        FAC:=FAC * CONST1;
    end;
    CONST1:=0.5 * (A + B);
    For J:=1 To N - 1 do
        For K:=N - 1 DownTo J do
            D[K]:=D[K] - CONST1 * D[K + 1];
end;
end.