testratinterpolation.pas

maths lib with source

unit testratinterpolation;
interface
uses Math, Ap, Sysutils, polinterpolation, ratinterpolation;

function TestRI(Silent : Boolean):Boolean;
function testratinterpolation_test_silent():Boolean;
function testratinterpolation_test():Boolean;

implementation

function FTask(T : Double; Task : Integer):Double;forward;
procedure BuildFloaterHormannTestTable();forward;


function TestRI(Silent : Boolean):Boolean;
var
    WasErrors : Boolean;
    NPErrors : Boolean;
    X : TReal1DArray;
    Y : TReal1DArray;
    W : TReal1DArray;
    W2 : TReal1DArray;
    N : Integer;
    M : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    D : Integer;
    Pass : Integer;
    PassCount : Integer;
    Err : Double;
    MaxErr : Double;
    T : Double;
    A : Double;
    B : Double;
    S : Double;
    Info : Integer;
begin
    WasErrors := False;
    PassCount := 20;
    
    //
    // Testing "No Poles" interpolation
    //
    NPErrors := False;
    MaxErr := 0;
    N:=2;
    while N<=10 do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        SetLength(W, N-1+1);
        SetLength(W2, N-1+1);
        
        //
        // D=1, non-equidistant nodes
        //
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Initialize X, Y, W
            //
            A := -1-1*RandomReal;
            B := +1+1*RandomReal;
            I:=0;
            while I<=N-1 do
            begin
                X[I] := ArcTan((B-A)*I/(N-1)+A);
                Inc(I);
            end;
            W[0] := -1/(X[1]-X[0]);
            S := 1;
            I:=1;
            while I<=N-2 do
            begin
                W[I] := S*(1/(X[I]-X[I-1])+1/(X[I+1]-X[I]));
                S := -S;
                Inc(I);
            end;
            W[N-1] := S/(X[N-1]-X[N-2]);
            
            //
            // Mix
            //
            I:=0;
            while I<=N-1 do
            begin
                K := RandomInteger(N);
                if K<>I then
                begin
                    T := X[I];
                    X[I] := X[K];
                    X[K] := T;
                    T := Y[I];
                    Y[I] := Y[K];
                    Y[K] := T;
                    T := W[I];
                    W[I] := W[K];
                    W[K] := T;
                end;
                Inc(I);
            end;
            
            //
            // Build and test
            //
            BuildFloaterHormannRationalInterpolant(X, N, 1, W2);
            S := W[0]/W2[0];
            APVMul(@W2[0], 0, N-1, S);
            I:=0;
            while I<=N-1 do
            begin
                MaxErr := Max(MaxErr, AbsReal((W[I]-W2[I])/W[I]));
                Inc(I);
            end;
            Inc(Pass);
        end;
        
        //
        // D = 0, 1, 2. Equidistant nodes.
        //
        D:=0;
        while D<=2 do
        begin
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Skip incorrect (N,D) pairs
                //
                if N<2*D then
                begin
                    Inc(Pass);
                    Continue;
                end;
                
                //
                // Initialize X, Y, W
                //
                A := -1-1*RandomReal;
                B := +1+1*RandomReal;
                I:=0;
                while I<=N-1 do
                begin
                    X[I] := (B-A)*I/(N-1)+A;
                    Inc(I);
                end;
                S := 1;
                if D=0 then
                begin
                    I:=0;
                    while I<=N-1 do
                    begin
                        W[I] := S;
                        S := -S;
                        Inc(I);
                    end;
                end;
                if D=1 then
                begin
                    W[0] := -S;
                    I:=1;
                    while I<=N-2 do
                    begin
                        W[I] := 2*S;
                        S := -S;
                        Inc(I);
                    end;
                    W[N-1] := S;
                end;
                if D=2 then
                begin
                    W[0] := S;
                    W[1] := -3*S;
                    I:=2;
                    while I<=N-3 do
                    begin
                        W[I] := 4*S;
                        S := -S;
                        Inc(I);
                    end;
                    W[N-2] := 3*S;
                    W[N-1] := -S;
                end;
                
                //
                // Mix
                //
                I:=0;
                while I<=N-1 do
                begin
                    K := RandomInteger(N);
                    if K<>I then
                    begin
                        T := X[I];
                        X[I] := X[K];
                        X[K] := T;
                        T := Y[I];
                        Y[I] := Y[K];
                        Y[K] := T;
                        T := W[I];
                        W[I] := W[K];
                        W[K] := T;
                    end;
                    Inc(I);
                end;
                
                //
                // Build and test
                //
                BuildFloaterHormannRationalInterpolant(X, N, D, W2);
                S := W[0]/W2[0];
                APVMul(@W2[0], 0, N-1, S);
                I:=0;
                while I<=N-1 do
                begin
                    MaxErr := Max(MaxErr, AbsReal((W[I]-W2[I])/W[I]));
                    Inc(I);
                end;
                Inc(Pass);
            end;
            Inc(D);
        end;
        Inc(N);
    end;
    if MaxErr>100*MachineEpsilon then
    begin
        NPErrors := True;
    end;
    
    //
    // report
    //
    WasErrors := NPErrors;
    if  not Silent then
    begin
        Write(Format('TESTING RATIONAL INTERPOLATION'#13#10'',[]));
        Write(Format('FLOATER-HORMANN TEST                     ',[]));
        if NPErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        if WasErrors then
        begin
            Write(Format('TEST FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST PASSED'#13#10'',[]));
        end;
        Write(Format(''#13#10''#13#10'',[]));
    end;
    
    //
    // end
    //
    Result :=  not WasErrors;
end;


function FTask(T : Double; Task : Integer):Double;
begin
    if (Task=0) or (Task=1) then
    begin
        Result := Sin(T);
    end
    else
    begin
        Result := 1/(1+Sqr(T));
    end;
end;


procedure BuildFloaterHormannTestTable();
var
    I : Integer;
    J : Integer;
    N : Integer;
    D : Integer;
    EMatrix : TReal1DArray;
    EMin : Double;
    DMin : Integer;
    Err : Double;
    T : Double;
    X : TReal1DArray;
    Y : TReal1DArray;
    W : TReal1DArray;
    A : Double;
    B : Double;
    NTbl : TInteger1DArray;
    NCnt : Integer;
    Task : Integer;
begin
    NCnt := 4;
    SetLength(NTbl, NCnt-1+1);
    NTbl[0] := 11;
    NTbl[1] := 21;
    NTbl[2] := 41;
    NTbl[3] := 81;
    A := -5;
    B := 5;
    
    //
    // Output
    //
    Task:=0;
    while Task<=3 do
    begin
        if (Task=0) or (Task=1) then
        begin
            Write(Format('TESTING f=sin(x)'#13#10'',[]));
        end
        else
        begin
            Write(Format(''#13#10''#13#10'TESTING f=1/(1+x^2)'#13#10'',[]));
        end;
        if (Task=0) or (Task=2) then
        begin
            Write(Format('EQUIDISTANT NODES'#13#10'',[]));
        end
        else
        begin
            Write(Format('CHEBYSHEV NODES'#13#10'',[]));
        end;
        Write(Format('N     D=3       D=5       D=8       D=N       Optimal D'#13#10'',[]));
        J:=0;
        while J<=NCnt-1 do
        begin
            N := NTbl[J];
            SetLength(X, N-1+1);
            SetLength(Y, N-1+1);
            SetLength(W, N-1+1);
            SetLength(EMatrix, N-1+1);
            
            //
            // Fill X, Y
            //
            I:=0;
            while I<=N-1 do
            begin
                if (Task=0) or (Task=2) then
                begin
                    X[I] := A+(B-A)*I/(N-1);
                end
                else
                begin
                    X[I] := 0.5*(A+B)+0.5*(B-A)*Cos(Pi*(2*I+1)/(2*(N-1)+2));
                end;
                Y[I] := FTask(X[I], Task);
                Inc(I);
            end;
            
            //
            // Test
            //
            EMin := MaxRealNumber;
            DMin := 0;
            D:=0;
            while D<=N-1 do
            begin
                BuildFloaterHormannRationalInterpolant(X, N, D, W);
                Err := 0;
                T := A;
                while T<=B*1.000001 do
                begin
                    Err := Max(Err, AbsReal(BarycentricInterpolation(X, Y, W, N, T)-FTask(T, Task)));
                    T := T+0.0005*(B-A);
                end;
                if Err<EMin then
                begin
                    EMin := Err;
                    DMin := D;
                end;
                EMatrix[D] := Err;
                Inc(D);
            end;
            
            //
            // Output
            //
            Write(Format('%3d   %8.2e  %8.2e  %8.2e  %8.2e  %8.2e (D=%0d)'#13#10'',[
                N-1,
                EMatrix[3],
                EMatrix[5],
                EMatrix[8],
                EMatrix[N-1],
                EMin,
                DMin]));
            Inc(J);
        end;
        Inc(Task);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testratinterpolation_test_silent():Boolean;
begin
    Result := TestRI(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testratinterpolation_test():Boolean;
begin
    Result := TestRI(False);
end;


end.