XProfan

Window Title "Interpolation zw. Stützwerten through Polynomkoeffizienten-Anpassung"
Window Style 24:randomize:CLS rnd(8^8):font 2:set("decimals",18)
'{********************************************************
'*       Polynomial Interpolation or Extrapolation       *
'*              of a Discreet Function F(x)              *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'* (Example: Function sin(x) - 2*cos(x) is given by 12   *
'*          points from x=0 to x=1.1.                    *
'*          Extrapolate for x=1.255).                    *
'*                                                       *
'*  For X             =  1.255                           *
'*  Estimated Y value =  .3294023272245815               *
'*  Estimated Error   = -8.273064603451457E-11           *
'*  Exact Y value     =  .3294023272200048               *
'*                                                       *
'* ----------------------------------------------------- *
'* REFERENCE: "Numerical Recipes, The manner of Scientific  *
'*             Computing By W.H. Press, B.P. Flannery,   *
'*             s.A. Teukolsky and W.T. Vetterling,       *
'*             cambridge University Press, 1986"         *
'*                                                       *
'*                  Basic Release By J-P Moreau, Paris.  *
'*                           (www.jpmoreau.fr)           *
'*********************************************************
'*                                                       *
'*      XProfan-Version  2014-10 by P.woodpecker, Wien       *
'*                                                       *
'*********************************************************
'}
' PROGRAM TEST_POLINT
Var n&=12' Number of points
Declare X![N&],Y![N&],C![N&],D![N&]
Declare i&,x1!,xx!,fct!,yy!,DY!

REPEAT

    ' The Stützwerte müssten NOT absolutely in equal intervals lying!
    ' define tables X and Y 'ACHTUNG: ARRAY IS MIT BASISINDEX 1 GEFÜHRT!
    X![1] = 0.0
    X![2] = 0.1
    X![3] = 0.2
    X![4] = 0.3
    X![5] = 0.4
    X![6] = 0.5
    X![7] = 0.6
    X![8] = 0.7
    X![9] = 0.8
    X![10]= 0.9
    X![11]= 1.0
    X![12]= 1.1

    Whileloop n&:i&=&Loop

        X1! = X![I&]
        FCT!=FCT(X1!)
        Y![I&] = FCT!

    Endwhile

    proc FCT :parameters x1!

        ' FUNCTION FCT(X1) ' but not Tabelleneingabe the Y-Stützwerte
        ' becomes here a known function herangezogen.
        ' the allows a check the accuracy the Interpolation
        FCT! = SIN(X1!) - 2.0 * COS(X1!)
        RETURN FCT!

    endproc

    ' ANWENDUNG DES GEFUNDENEN POLYNOMS
    ' default one X-Wertes and query the intern found Interpolationsformel
    print "\n EINGABE:  X-worth, for Y To interpolate is "
    print " (with X=0 becomes eingebauter Testwert 1.255 uses) X = ";
    input xx! : case xx!=0 : XX! = 1.255
    ' INTERPOLATION
    yy!=POLINT(X1!,N&,XX!,YY!)
    ' AUSGABE
    case %csrlin>20:cls rnd(8^8)
    PRINT
    PRINT "     for the sought X = ";stature$("%g",XX!)
    PRINT "  Interpolierter Y-worth = ";stature$("%g",YY!)
    PRINT "       latest Korrektur = ";stature$("%g",DY!)
    X1! = XX! : FCT!=FCT(X1!)
    PRINT " Exakter Vergleichswert = ";stature$("%g",FCT!)
    PRINT "--------------------------------------------------\n"

UNTIL 0

proc STOP :sound 2000,100: waitinput:END

endproc

Proc POLINT :parameters X!,N&,XX!,YY!

    '*****************************************************
    '  Origianl-Subroutine: POLINT(X,Y,n,XX,YY,DY)       *
    '*****************************************************
    '*     Polynomial Interpolation or Extrapolation     *
    '*            of a Discreet Function                 *
    '* ------------------------------------------------- *
    '* INPUTS:                                           *
    '*    X:    Table of abscissas (n)                   *
    '*    Y:    Table of ordinates (n)                   *
    '*    n:    Number of points                         *
    '*   XX:    Interpolation abscissa value             *
    '* OUTPUT:                                           *
    '*   YY:    Returned estimation of function for X    *
    '*   DY:    Estimated error for YY                   *
    '*****************************************************
    Declare NS&,dif!,dift!,C![n&],D![n&],m&,ho!,hp!,w!,whom!
    NS& = 1
    DIF! = ABS(XX! - X![1])

    whileloop n&:i&=&Loop

        DIFT! = ABS(XX! - X![1])

        IF DIFT! < DIF!

            NS& = I&'index of closest table entry
            DIF! = DIFT!

        ENDIF

        C![I&] = Y![I&]'Initialize the C"s and D"s
        D![I&] = Y![I&]

    endwhile

    YY! = Y![NS&]'Initial approximation of Y
    NS& = NS& - 1

    whileloop n&-1:m&=&Loop

        whileloop n&-m&:i&=&Loop

            HO! = X![I&] - XX!
            HP! = X![I& + M&] - XX!
            W! = C![I& + 1] - D![I&]
            DEN! = HO! - HP!

            IF DEN! = 0

                PRINT
                PRINT " *** FEHLER: ZWEI STÜTZWERTE WIDERSPRECHEN SICH! *** "
                STOP

            ENDIF

            DEN! = W! / DEN!
            D![I&] = HP! * DEN!'Update the C's and D's
            C![I&] = HO! * DEN!

        endwhile

        IF (2*NS&) < (N&-M&)' After each column in the tableau XA is completed,

            DY! = C![NS&+1]' we decide which correction, C or D, we want to

        ELSE' add to our accumulating value of Y, i.e. which

            DY! = D![NS&]' path to take through the tableau, forking up or
            NS& = NS& - 1' down. We do this in such a way as to take the

        ENDIF' most "straight line" route through the tableau to

        ' its apex, updating NS accordingly to keep track
        YY! = YY! + DY!' of where we are. Diese route keeps the partial

    endwhile' approximations centered (insofar as possible) on

    ' the target X.The last DY added is thus the error
    RETURN YY!' indication.

endproc