testpolynomialinterpolationunit.pas

maths lib with source

unit testpolynomialinterpolationunit;
interface
uses Math, Ap, Sysutils, polinterpolation;

function TestPolynomialInterpolation(Silent : Boolean):Boolean;
function testpolynomialinterpolationunit_test_silent():Boolean;
function testpolynomialinterpolationunit_test():Boolean;

implementation

function PolInt(const x : TReal1DArray;
     F : TReal1DArray;
     n : Integer;
     t : Double):Double;forward;


function TestPolynomialInterpolation(Silent : Boolean):Boolean;
var
    WasErrors : Boolean;
    GuardedAccErrors : Boolean;
    GenIntErrors : Boolean;
    NevIntErrors : Boolean;
    NevDIntErrors : Boolean;
    EquidistantIntErrors : Boolean;
    Chebyshev1IntErrors : Boolean;
    Chebyshev2IntErrors : Boolean;
    EpsTol : Double;
    N : Integer;
    MaxN : Integer;
    I : Integer;
    J : Integer;
    Pass : Integer;
    PassCount : Integer;
    A : Double;
    B : Double;
    P : Double;
    DP : Double;
    D2P : Double;
    Q : Double;
    DQ : Double;
    D2Q : Double;
    X : TReal1DArray;
    Y : TReal1DArray;
    W : TReal1DArray;
    V1 : Double;
    V2 : Double;
    V3 : Double;
    H : Double;
    MaxErr : Double;
begin
    WasErrors := False;
    MaxN := 5;
    PassCount := 20;
    EpsTol := 1.0E+4;
    
    //
    // Testing general barycentric formula
    //
    
    //
    // Crash test: three possible types of overflow
    //
    SetLength(X, 2+1);
    SetLength(Y, 2+1);
    SetLength(W, 2+1);
    X[0] := -1;
    Y[0] := 1;
    W[0] := 1;
    X[1] := 0;
    Y[1] := 0.002;
    W[1] := 1;
    X[2] := 1;
    Y[2] := 1;
    W[2] := 1;
    V2 := BarycentricInterpolation(X, Y, W, 3, 0);
    V2 := EquidistantPolynomialInterpolation(-1.0, 1.0, Y, 3, 0);
    V1 := 0.5;
    while V1>0 do
    begin
        V2 := BarycentricInterpolation(X, Y, W, 3, V1);
        V2 := EquidistantPolynomialInterpolation(-1.0, 1.0, Y, 3, V1);
        V1 := 0.5*V1;
    end;
    Y[1] := 200;
    V1 := 0.5;
    while V1>0 do
    begin
        V2 := BarycentricInterpolation(X, Y, W, 3, V1);
        V2 := EquidistantPolynomialInterpolation(-1.0, 1.0, Y, 3, V1);
        V1 := 0.5*V1;
    end;
    
    //
    // Accuracy test of guarded formula (near-overflow task)
    //
    SetLength(X, 2+1);
    SetLength(Y, 2+1);
    SetLength(W, 2+1);
    X[0] := -1;
    Y[0] := -1;
    W[0] := 1;
    X[1] := 0;
    Y[1] := 0;
    W[1] := -2;
    X[2] := 1;
    Y[2] := 1;
    W[2] := 1;
    V1 := 0.5;
    MaxErr := 0;
    while V1>0 do
    begin
        V2 := BarycentricInterpolation(X, Y, W, 3, V1);
        MaxErr := Max(MaxErr, AbsReal(V2-V1)/AbsReal(V1));
        V2 := EquidistantPolynomialInterpolation(-1.0, 1.0, Y, 3, V1);
        MaxErr := Max(MaxErr, AbsReal(V2-V1)/AbsReal(V1));
        V1 := 0.5*V1;
    end;
    GuardedAccErrors := MaxErr>EpsTol*MachineEpsilon;
    
    //
    // Testing polynomial interpolation using a
    // general barycentric formula. Test on Chebyshev nodes
    // (to avoid big condition numbers).
    //
    MaxErr := 0;
    H := 0.0001;
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        SetLength(W, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            V1 := 0.7*RandomReal+0.3;
            I:=0;
            while I<=N-1 do
            begin
                X[I] := Cos(I*Pi/N);
                Y[I] := Cos(V1*X[I])+2.3*X[I]-2;
                Inc(I);
            end;
            I:=0;
            while I<=N-1 do
            begin
                W[I] := 1.0;
                J:=0;
                while J<=N-1 do
                begin
                    if J<>I then
                    begin
                        W[I] := W[I]/(X[I]-X[J]);
                    end;
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test at intermediate points
            //
            V1 := -1;
            while V1<=1.0001 do
            begin
                V2 := PolInt(X, Y, N, V1);
                V3 := BarycentricInterpolation(X, Y, W, N, V1);
                MaxErr := Max(MaxErr, AbsReal(V2-V3));
                V1 := V1+0.01;
            end;
            
            //
            // Test at nodes
            //
            I:=0;
            while I<=N-1 do
            begin
                MaxErr := Max(MaxErr, AbsReal(BarycentricInterpolation(X, Y, W, N, X[I])-Y[I]));
                Inc(I);
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    GenIntErrors := MaxErr>EpsTol*MachineEpsilon;
    
    //
    // Testing polynomial interpolation using a Neville's algorithm.
    // Test on Chebyshev nodes (don't want to work with big condition numbers).
    //
    MaxErr := 0;
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            V1 := 0.7*RandomReal+0.3;
            I:=0;
            while I<=N-1 do
            begin
                X[I] := Cos(I*Pi/N);
                Y[I] := Cos(V1*X[I])+2.3*X[I]-2;
                Inc(I);
            end;
            
            //
            // Test
            //
            V1 := -1;
            while V1<=1.0001 do
            begin
                V2 := PolInt(X, Y, N, V1);
                V3 := NevilleInterpolation(X, Y, N, V1);
                MaxErr := Max(MaxErr, AbsReal(V2-V3));
                V1 := V1+0.01;
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    NevIntErrors := MaxErr>EpsTol*MachineEpsilon;
    
    //
    // Testing polynomial interpolation using an
    // equidistant barycentric algorithm.
    //
    MaxErr := 0;
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            V1 := 0.7*RandomReal+0.3;
            A := -(0.5+0.5*RandomReal);
            B := +(0.5+0.5*RandomReal);
            I:=0;
            while I<=N-1 do
            begin
                X[I] := A+(B-A)*I/Max(N-1, 1);
                Y[I] := Cos(V1*X[I])+2.3*X[I]-2;
                Inc(I);
            end;
            
            //
            // Test
            //
            V1 := -1;
            while V1<=1.0001 do
            begin
                V2 := PolInt(X, Y, N, V1);
                V3 := EquidistantPolynomialInterpolation(A, B, Y, N, V1);
                MaxErr := Max(MaxErr, AbsReal(V2-V3));
                V1 := V1+0.01;
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    EquidistantIntErrors := MaxErr>EpsTol*MachineEpsilon;
    
    //
    // Testing polynomial interpolation using
    // Chebyshev-1 barycentric algorithm.
    //
    MaxErr := 0;
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            V1 := 0.7*RandomReal+0.3;
            A := -(0.5+0.5*RandomReal);
            B := +(0.5+0.5*RandomReal);
            I:=0;
            while I<=N-1 do
            begin
                X[I] := 0.5*(A+B)+0.5*(B-A)*Cos(Pi*(2*I+1)/(2*(N-1)+2));
                Y[I] := Cos(V1*X[I])+2.3*X[I]-2;
                Inc(I);
            end;
            
            //
            // Test
            //
            V1 := A;
            while V1<=B do
            begin
                V2 := PolInt(X, Y, N, V1);
                V3 := Chebyshev1Interpolation(A, B, Y, N, V1);
                MaxErr := Max(MaxErr, AbsReal(V2-V3));
                V1 := V1+0.01;
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    Chebyshev1IntErrors := MaxErr>EpsTol*MachineEpsilon;
    
    //
    // Testing polynomial interpolation using
    // Chebyshev-2 barycentric algorithm.
    //
    MaxErr := 0;
    N:=2;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            V1 := 0.7*RandomReal+0.3;
            A := -(0.5+0.5*RandomReal);
            B := +(0.5+0.5*RandomReal);
            I:=0;
            while I<=N-1 do
            begin
                X[I] := 0.5*(A+B)+0.5*(B-A)*Cos(Pi*I/(N-1));
                Y[I] := Cos(V1*X[I])+2.3*X[I]-2;
                Inc(I);
            end;
            
            //
            // Test
            //
            V1 := A;
            while V1<=B do
            begin
                V2 := PolInt(X, Y, N, V1);
                V3 := Chebyshev2Interpolation(A, B, Y, N, V1);
                MaxErr := Max(MaxErr, AbsReal(V2-V3));
                V1 := V1+0.01;
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    Chebyshev2IntErrors := MaxErr>EpsTol*MachineEpsilon;
    
    //
    // Testing polynomial interpolation using a Neville's algorithm.
    // with derivatives. Test on Chebyshev nodes (don't want to work
    // with big condition numbers).
    //
    MaxErr := 0;
    H := 0.0001;
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            V2 := 0.7*RandomReal+0.3;
            I:=0;
            while I<=N-1 do
            begin
                X[I] := Cos(I*Pi/N);
                Y[I] := Cos(V2*X[I]);
                Inc(I);
            end;
            
            //
            // Test
            //
            V1 := -1;
            while V1<=1.0001 do
            begin
                P := PolInt(X, Y, N, V1);
                DP := (PolInt(X, Y, N, V1+H)-PolInt(X, Y, N, V1-H))/(2*H);
                D2P := (PolInt(X, Y, N, V1-H)-2*PolInt(X, Y, N, V1)+PolInt(X, Y, N, V1+H))/Sqr(H);
                NevilleDifferentiation(X, Y, N, V1, Q, DQ, D2Q);
                MaxErr := Max(MaxErr, AbsReal(P-Q));
                MaxErr := Max(MaxErr, AbsReal(DP-DQ));
                MaxErr := Max(MaxErr, AbsReal(D2P-D2Q));
                V1 := V1+0.01;
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    NevDIntErrors := MaxErr>0.001;
    
    //
    // report
    //
    WasErrors := GuardedAccErrors or GenIntErrors or NevIntErrors or EquidistantIntErrors or Chebyshev1IntErrors or Chebyshev2IntErrors or NevDIntErrors;
    if  not Silent then
    begin
        Write(Format('TESTING POLYNOMIAL INTERPOLATION'#13#10'',[]));
        
        //
        // Normal tests
        //
        Write(Format('BARYCENTRIC INTERPOLATION TEST:          ',[]));
        if GenIntErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('NEVILLE INTERPOLATON TEST:               ',[]));
        if NevIntErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('EQUIDISTANT INTERPOLATON TEST:           ',[]));
        if EquidistantIntErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('CHEBYSHEV-1 INTERPOLATON TEST:           ',[]));
        if Chebyshev1IntErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('CHEBYSHEV-2 INTERPOLATON TEST:           ',[]));
        if Chebyshev2IntErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('NEVILLE DIFFERENTIATION TEST:            ',[]));
        if NevDIntErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('GUARDED FORMULA ACCURACY TEST:           ',[]));
        if GuardedAccErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        
        //
        // if we are here then crash test is passed :)
        //
        Write(Format('CRASH TEST:                              OK'#13#10'',[]));
        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 PolInt(const x : TReal1DArray;
     F : TReal1DArray;
     n : Integer;
     t : Double):Double;
var
    I : Integer;
    J : Integer;
begin
    F := DynamicArrayCopy(F);
    N := N-1;
    j:=0;
    while j<=n-1 do
    begin
        i:=j+1;
        while i<=n do
        begin
            F[i] := ((t-x[j])*F[i]-(t-x[i])*F[j])/(x[i]-x[j]);
            Inc(i);
        end;
        Inc(j);
    end;
    Result := F[n];
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testpolynomialinterpolationunit_test_silent():Boolean;
begin
    Result := TestPolynomialInterpolation(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testpolynomialinterpolationunit_test():Boolean;
begin
    Result := TestPolynomialInterpolation(False);
end;


end.