testsplineinterpolationunit.cs

三样条插值函数算法

            //
            a = -1-AP.Math.RandomReal();
            b = +1+AP.Math.RandomReal();
            for(i=0; i<=n-1; i++)
            {
                x[i] = a+(b-a)*i/(n-1);
                y[i] = Math.Cos(1.3*Math.PI*x[i]+0.4);
            }
            spline3.buildcubicspline(x, y, n, 2, 0.0, 2, 0.0, ref c);
            
            //
            // Test diff
            //
            err = 0;
            for(pass=1; pass<=passcount; pass++)
            {
                t = a+(b-a)*AP.Math.RandomReal();
                spline3.splinedifferentiation(ref c, t, ref s, ref ds, ref d2s);
                vl = spline3.splineinterpolation(ref c, t-h);
                vm = spline3.splineinterpolation(ref c, t);
                vr = spline3.splineinterpolation(ref c, t+h);
                err = Math.Max(err, Math.Abs(s-vm));
                err = Math.Max(err, Math.Abs(ds-(vr-vl)/(2*h)));
                err = Math.Max(err, Math.Abs(d2s-(vr-2*vm+vl)/AP.Math.Sqr(h)));
            }
            dserrors = dserrors | err>0.001;
            
            //
            // Test copy
            //
            spline3.splinecopy(ref c, ref c2);
            err = 0;
            for(pass=1; pass<=passcount; pass++)
            {
                t = a+(b-a)*AP.Math.RandomReal();
                err = Math.Max(err, Math.Abs(spline3.splineinterpolation(ref c, t)-spline3.splineinterpolation(ref c2, t)));
            }
            cperrors = cperrors | err>100*AP.Math.MachineEpsilon;
            
            //
            // Test unpack
            //
            uperrors = uperrors | !testunpack(ref c, ref x);
            
            //
            // Test lin.trans.
            //
            err = 0;
            for(pass=1; pass<=passcount; pass++)
            {
                
                //
                // LinTransX, general A
                //
                sa = 4*AP.Math.RandomReal()-2;
                sb = 2*AP.Math.RandomReal()-1;
                t = a+(b-a)*AP.Math.RandomReal();
                spline3.splinecopy(ref c, ref c2);
                spline3.splinelintransx(ref c2, sa, sb);
                err = Math.Max(err, Math.Abs(spline3.splineinterpolation(ref c, t)-spline3.splineinterpolation(ref c2, (t-sb)/sa)));
                
                //
                // LinTransX, special case: A=0
                //
                sb = 2*AP.Math.RandomReal()-1;
                t = a+(b-a)*AP.Math.RandomReal();
                spline3.splinecopy(ref c, ref c2);
                spline3.splinelintransx(ref c2, 0, sb);
                err = Math.Max(err, Math.Abs(spline3.splineinterpolation(ref c, sb)-spline3.splineinterpolation(ref c2, t)));
                
                //
                // LinTransY
                //
                sa = 2*AP.Math.RandomReal()-1;
                sb = 2*AP.Math.RandomReal()-1;
                t = a+(b-a)*AP.Math.RandomReal();
                spline3.splinecopy(ref c, ref c2);
                spline3.splinelintransy(ref c2, sa, sb);
                err = Math.Max(err, Math.Abs(sa*spline3.splineinterpolation(ref c, t)+sb-spline3.splineinterpolation(ref c2, t)));
            }
            lterrors = lterrors | err>100*AP.Math.MachineEpsilon;
        }
        
        //
        // Testing integration
        //
        err = 0;
        for(n=20; n<=35; n++)
        {
            x = new double[n-1+1];
            y = new double[n-1+1];
            for(pass=1; pass<=passcount; pass++)
            {
                
                //
                // Prepare cubic spline
                //
                a = -1-0.2*AP.Math.RandomReal();
                b = +1+0.2*AP.Math.RandomReal();
                for(i=0; i<=n-1; i++)
                {
                    x[i] = a+(b-a)*i/(n-1);
                    y[i] = Math.Sin(Math.PI*x[i]+0.4)+Math.Exp(x[i]);
                }
                bl = Math.PI*Math.Cos(Math.PI*a+0.4)+Math.Exp(a);
                br = Math.PI*Math.Cos(Math.PI*b+0.4)+Math.Exp(b);
                spline3.buildcubicspline(x, y, n, 1, bl, 1, br, ref c);
                
                //
                // Test
                //
                t = a+(b-a)*AP.Math.RandomReal();
                v = -(Math.Cos(Math.PI*a+0.4)/Math.PI)+Math.Exp(a);
                v = -(Math.Cos(Math.PI*t+0.4)/Math.PI)+Math.Exp(t)-v;
                v = v-spline3.splineintegration(ref c, t);
                err = Math.Max(err, Math.Abs(v));
            }
        }
        ierrors = ierrors | err>0.001;
        
        //
        // report
        //
        waserrors = lserrors | cserrors | hserrors | aserrors | dserrors | cperrors | uperrors | lterrors | ierrors;
        if( !silent )
        {
            System.Console.Write("TESTING SPLINE INTERPOLATION");
            System.Console.WriteLine();
            
            //
            // Normal tests
            //
            System.Console.Write("LINEAR SPLINE TEST:                      ");
            if( lserrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("CUBIC SPLINE TEST:                       ");
            if( cserrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("HERMITE SPLINE TEST:                     ");
            if( hserrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("AKIMA SPLINE TEST:                       ");
            if( aserrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("DIFFERENTIATION TEST:                    ");
            if( dserrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("COPY TEST:                               ");
            if( cperrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("UNPACK TEST:                             ");
            if( uperrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("LIN.TRANS. TEST:                         ");
            if( lterrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            System.Console.Write("INTEGRATION TEST:                        ");
            if( ierrors )
            {
                System.Console.Write("FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("OK");
                System.Console.WriteLine();
            }
            if( waserrors )
            {
                System.Console.Write("TEST FAILED");
                System.Console.WriteLine();
            }
            else
            {
                System.Console.Write("TEST PASSED");
                System.Console.WriteLine();
            }
            System.Console.WriteLine();
            System.Console.WriteLine();
        }
        
        //
        // end
        //
        result = !waserrors;
        return result;
    }


    /*************************************************************************
    Lipschitz constants for spline inself, first and second derivatives.
    *************************************************************************/
    private static void lconst(double a,
        double b,
        ref double[] c,
        double lstep,
        ref double l0,
        ref double l1,
        ref double l2)
    {
        double t = 0;
        double vl = 0;
        double vm = 0;
        double vr = 0;
        double prevf = 0;
        double prevd = 0;
        double prevd2 = 0;
        double f = 0;
        double d = 0;
        double d2 = 0;

        l0 = 0;
        l1 = 0;
        l2 = 0;
        t = a-0.1;
        vl = spline3.splineinterpolation(ref c, t-2*lstep);
        vm = spline3.splineinterpolation(ref c, t-lstep);
        vr = spline3.splineinterpolation(ref c, t);
        f = vm;
        d = (vr-vl)/(2*lstep);
        d2 = (vr-2*vm+vl)/AP.Math.Sqr(lstep);
        while( t<=b+0.1 )
        {
            prevf = f;
            prevd = d;
            prevd2 = d2;
            vl = vm;
            vm = vr;
            vr = spline3.splineinterpolation(ref c, t+lstep);
            f = vm;
            d = (vr-vl)/(2*lstep);
            d2 = (vr-2*vm+vl)/AP.Math.Sqr(lstep);
            l0 = Math.Max(l0, Math.Abs((f-prevf)/lstep));
            l1 = Math.Max(l1, Math.Abs((d-prevd)/lstep));
            l2 = Math.Max(l2, Math.Abs((d2-prevd2)/lstep));
            t = t+lstep;
        }
    }


    /*************************************************************************
    Lipschitz constants for spline inself, first and second derivatives.
    *************************************************************************/
    private static bool testunpack(ref double[] c,
        ref double[] x)
    {
        bool result = new bool();
        int i = 0;
        int n = 0;
        double err = 0;
        double t = 0;
        double v1 = 0;
        double v2 = 0;
        int pass = 0;
        int passcount = 0;
        double[,] tbl = new double[0,0];

        passcount = 20;
        err = 0;
        spline3.splineunpack(ref c, ref n, ref tbl);
        for(i=0; i<=n-2; i++)
        {
            for(pass=1; pass<=passcount; pass++)
            {
                t = AP.Math.RandomReal()*(tbl[i,1]-tbl[i,0]);
                v1 = tbl[i,2]+t*tbl[i,3]+AP.Math.Sqr(t)*tbl[i,4]+t*AP.Math.Sqr(t)*tbl[i,5];
                v2 = spline3.splineinterpolation(ref c, tbl[i,0]+t);
                err = Math.Max(err, Math.Abs(v1-v2));
            }
        }
        for(i=0; i<=n-2; i++)
        {
            err = Math.Max(err, Math.Abs(x[i]-tbl[i,0]));
        }
        for(i=0; i<=n-2; i++)
        {
            err = Math.Max(err, Math.Abs(x[i+1]-tbl[i,1]));
        }
        result = err<100*AP.Math.MachineEpsilon;
        return result;
    }


    /*************************************************************************
    Silent unit test
    *************************************************************************/
    public static bool testsplineinterpolationunit_test_silent()
    {
        bool result = new bool();

        result = testsplineinterpolation(true);
        return result;
    }


    /*************************************************************************
    Unit test
    *************************************************************************/
    public static bool testsplineinterpolationunit_test()
    {
        bool result = new bool();

        result = testsplineinterpolation(false);
        return result;
    }
}