XProfan

WindowTitle "Interpolation zw. Stützwerten durch Polynomkoeffizienten-Anpassung"
WindowStyle 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 Art 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.Specht, 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

    ' Die Stützwerte müssten NICHT unbedingt in gleichen Abständen liegen!
    ' define tables X and Y 'ACHTUNG: ARRAY WIRD 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) ' Statt Tabelleneingabe der Y-Stützwerte
        ' wird hier eine bekannte Funktion herangezogen.
        ' Das erlaubt eine Prüfung der Genauigkeit der Interpolation
        FCT! = SIN(X1!) - 2.0 * COS(X1!)
        RETURN FCT!

    endproc

    ' ANWENDUNG DES GEFUNDENEN POLYNOMS
    ' Vorgabe eines X-Wertes und Abfrage der intern gefundenen Interpolationsformel
    print "\n EINGABE:  X-Wert, per den Y zu interpolieren ist "
    print " (Bei X=0 wird eingebauter Testwert 1.255 verwendet) 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 "     Für das gesuchte X = ";format$("%g",XX!)
    PRINT "  Interpolierter Y-Wert = ";format$("%g",YY!)
    PRINT "       Letzte Korrektur = ";format$("%g",DY!)
    X1! = XX! : FCT!=FCT(X1!)
    PRINT " Exakter Vergleichswert = ";format$("%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!,den!
    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. This 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