XProfan

WindowTitle "Tschebyschow-Approximation einer benutzerdefinierten Funktion"
' Quelle: https://jean-pierre.moreau.pagesperso-orange.fr/Basic/tchebysh_bas.txt
' Transponiert nach XProfan 11.2a (D) Demo by P.Specht, Vienna/Austria
' No warranty whatsoever! Ohne jegliche Gewähr! Nutzung auf eigene Gefahr!
'********************************************************
'* Chebyshev Approximation of a user defined real       *
'* function FUNC(X) in double precision.                *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'* (Approximate sin(x) from x=0 to x=PI).               *
'*                                                      *
'* Chebyshev coefficients (N= 10 ):                     *
'*  0.944002431536470                                   *
'* -0.000000000000000                                   *
'* -0.499403258270407                                   *
'* -0.000000000000000                                   *
'*  0.027992079617546                                   *
'* -0.000000000000000                                   *
'* -0.000596695195801                                   *
'*  0.000000000000000                                   *
'*  0.000006704175524                                   *
'*  0.000000000000000                                   *
'*      X      Chebyshev Eval.   SIN(X)                 *
'*-----------------------------------------             *
'*  0.00000000   0.00000005   0.00000000                *
'*  0.34906585   0.34202018   0.34202014                *
'*  0.69813170   0.64278757   0.64278761                *
'*  1.04719755   0.86602545   0.86602540                *
'*  1.39626340   0.98480773   0.98480775                *
'*  1.74532925   0.98480773   0.98480775                *
'*  2.09439510   0.86602545   0.86602540                *
'*  2.44346095   0.64278757   0.64278761                *
'*  2.79252680   0.34202018   0.34202014                *
'*  3.14159265   0.00000005   0.00000000                *
'*                                                      *
'*                 Basic Release By J-P Moreau, Paris.  *
'*                          (www.jpmoreau.fr)           *
'*  XProfan 11-release by P.Specht, Vienna / Austria    *
'* ---------------------------------------------------- *
'* 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"        *
'*             [BIBLI 08].                              *
'********************************************************
'PROGRAM PART 1: TESTCHEBY
'DEFDBL A-H, O-Z
'DEFINT I-N
WindowStyle 24:Window 0,0-%maxx,%maxy-40:Font 2:Set("Decimals",17)
Declare N&,Zero!,Pi!,x0!,x1!, half!,two!,A!,B!,BmA!,BpA!
Declare i&,j&,k&,y!,xx!, func!
Declare fac!,sum! ,  f0$,f1$ ,dx!,x!
Declare m&,d!,dd!,y2!,sv!, chebev!
ZERO! = 0 : Pi! = 4*Arctan(1)':Print Pi!:WaitInput
N& = 10
Declare C![N&], F![N&]
X0! = ZERO!: X1! = Pi!
Gosub "ES1000"'call CHEBFT(X0,X1,C,N)
Print "\n\n Tschebyschow-Koeffizienten per Grad N = "; N&; \
":\n ----------------------------------------------------"
F0$ = "###0.###############;###0.###############;###0.###############"

WhileLoop n&:i&=&Loop

    Print "   C(";int(i&-1);") = ";Format$(F0$,C![I&])

EndWhile

WaitInput 5000
DX! = (X1!-X0!) / (N&-1)
X! = X0! - DX!
Font 0
PRINT "\n\n   x:    Tschebyschow-Entwicklung:   vgl.Orig SIN(X)     "
PRINT " ------------------------------------------------------------"
F1$ = "###0.########"

WhileLoop n&:i&=&Loop

    X! = X! + DX!
    ES2000'call CHEBEV(X0,X1,C,N,X)
    Print "   ";Format$(F1$, X!),
    Print Tab(20);Format$(F1$, CHEBEV!),
    Print Tab(40);Format$(F1$, Sin(X!) )

EndWhile

Beep: WaitInput
End'of main program
'user defined function FUNC(XX)

Proc ES500

    FUNC! = Sin(XX!)

EndProc

ES1000:
'subroutine CHEBFT(A,B,C,N)
'********************************************************
'* Chebyshev fit: Given a real function FUNC(X), lower  *
'* and upper limits of the interval [A,B] for X, and a  *
'* maximum degree N, this routine computes the N Cheby- *
'* shev coefficients Ck, such that FUNC(X) is approxima-*
'* ted by:  N                                           *
'*         [Sum Ck Tk-1(Y)] - C1/2, where X and Y are   *
'*         k=1                                          *
'* related by:     Y = (X - 1/2(A+B)) / (1/2(B-A))      *
'* This routine is to be used with moderately large N   *
'* (e.g. 30 or 50), the array of C's subsequently to be *
'* truncated at the smaller value m such that Cm+1 and  *
'* subsequent elements are negligible.                  *
'********************************************************
HALF! = 0.5: TWO! = 2
A! = X0!: B! = X1!
BMA! = HALF! * (B! - A!): BPA! = HALF! * (B! + A!)

WhileLoop n&:k&=&Loop

    Y! = Cos(Pi! * (K&-HALF!)/N&)
    XX! = Y! * BMA! + BpA!
    ES500
    F![K&] = FUNC!

EndWhile

FAC! = TWO! / N&

WhileLoop n&:j&=&Loop

    SUM! = ZERO!

    WhileLoop n&:k&=&Loop

        SUM! = SUM! + F![K&] * Cos((Pi!*(J&-1)) * ((K&-HALF!) / N&))

    EndWhile

    C![J&] = FAC! * SUM!

EndWhile

RETURN

Proc ES2000'function CHEBEV(A,B,C,M,X)

    '**********************************************************
    '* Chebyshev evaluation: All arguments are input. C is an *
    '* array of Chebyshev coefficients, of length M, the first*
    '* M elements of Coutput from subroutine CHEBFT (which    *
    '* must have been called with the same A and B). The Che- *
    '* byshev polynomial is evaluated at a point Y determined *
    '* from X, A and B, and the result FUNC(X) is returned as *
    '* the function value.                                    *
    '**********************************************************

    If ((X!-A!) * (X!-B!)) > ZERO!

        Print " X liegt nicht im Bereich!"
        Return

    EndIF

    M& = N&
    D! = ZERO!: Dd! = ZERO!
    Y! = (TWO!*X!-A!-B!) / (B!-A!)'change of variable
    Y2! = TWO! * Y!

    WhileLoop m&,2,-1:j&=&Loop

        SV! = D!
        D! = Y2! * D! - Dd! + C![J&]
        Dd! = SV!

    EndWhile

    CHEBEV! = Y! * D! - Dd! + HALF! * C![1]

EndProc