XProfan

Deutsch

p.specht

Integration - die zugrundeliegende Frage ist stets: Wie groß ist die Fläche unter einer bestimmten Funktionskurve y=f(x) zwischen den x-Grenzen a und b ? In der Praxis hat man allerdings nur y-Meßwerte zu bestimmten x-Achsenwerten, kennt also die zugrundeliegende mathematische Funktion nicht wirklich. Fast genial erscheint daher die Idee, das Gummilineal der Mathematik, die Polynomfunktion in die Meßpunkte einzupassen und anschließend zu Integrieren. Das geht, weil nämlich das Integral eines Polynoms formelmäßig sehr gut bekannt ist. Bloß die entsprechenen Koeffizienten fehlen noch - doch da können Kapazunder wie Kepler, Euler, Gauss, Lagrange, Cotes oder eben auch der Schottisch-Britische Mathematiker, Astronom und Universitätsrektor Sir David Gregory (1659-1708), ein Zeitgenosse Newtons, helfen - und zwar gründlich:

Die von ihm verwendeten Formeln sind vom Grad her einstellbar, sodaß er die Verfahren aller vorgenannter Herren beliebig emulieren kann. Weiters sind seine Koeffizientenformeln als Gewichte eines iterativen Verfahrens leicht implementierbar.

P.S.: Mathematische interessanter mögen Lagrange- und Tschebeyschow-Polynome sein, aber nicht praxisgerechter. Lediglich das Romberg-Verfahren erscheint noch etwas eleganter!

Kompilieren Markieren Separieren

WindowTitle " GREGORY-INTEGRATION einer eingebauten Testfunktion x^n "
WindowStyle 24:Window 0,0-%maxx,%maxy-40:font 0:set("decimals",17)
'---------------------------------------------------------------------
' (C) ACM-TOMS ALGORITHM 280, COLLECTED ALGORITHMS FROM ACM.
' THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM,
' VOL. 9, NO. 4, April, 1966, P.  271
'---------------------------------------------------------------------
' (D) Demo-Migration von C++ nach XProfan11.2a 2014-11
' by P.Specht, Wien (Austria). OHNE JEDE GEWÄHR!
'---------------------------------------------------------------------
'------------------------------------------------------------------------
' Programmanpassungen für Arraygröße und zu integrierende Funktion
'------------------------------------------------------------------------
Var nmax&=100' Verfahren ausgelegt für diese maximalen Arraygrößen

Proc Fn_PowN :parameters a!,n&

    ' y=x^n
    case n&=0:return 1

    if n&>0 : var x!=a!

        whileloop n&-1

            x!=x!*a!

        endwhile

        return x!

    else

        Print " Negativer Exponent in Fn_PowN nicht erlaubt!"
        beep:waitinput 20000:return 0

    endif

EndProc

'------------------------------------------------------------------------
'------------------------------------------------------------------------

Proc Gregory :parameters n&,r&,t![],w![]'findet Gewichtungskoeffizienten

    '  Computes the abscissas and weights of the Gregory quadrature rule
    '  with r differences: Given h = (t_n - t_0)/n, then holds:
    '  Integral_from t_0 to t_n of: f(t) dt =~
    '  	  h*(f_0/2 + f_1 + ... +f_n-1 + f_n/2 )
    '  	- h*12*( nabla(f_n) - delta(f_0) )
    '  	- h*24*( nabla^2(f_n) + delta^2(f_0))
    '   - ...
    '   - h*c[r+1] *( nabla^r*f(n) + delta^r*(f_0) )
    '=  Sum[j=0..n]of w[j]*f(t[j])
    'where h = (t_n - t_0)/n, and the c_j' are given in Henrici (1964, Schweiz & USA)
    'The number r must be an integer from 0 to n, the number of subdivisions.
    'The left and right endpoints must be in t(0) and t(n) respectively.
    'The abscissas are returned in t(0) to t(n) and the corresponding weights
    'in w(0) to w(n).
    '
    'If r=0 the Gregory rule is the same as the repeated trapezoid rule, and
    'if r=n the same as the Newton-Cotes rule (closed type). The order p of the
    'quadrature rule is p = r+1 for r odd and p = r+2 for r even.
    'For n >= 9 and large r some of the weights can be negative.
    '
    'For n <= 32 and r<= 24, the numerical integration of powers (less than r)
    'of x on the interval [0,1] gave 9 significant digits correct in an 11
    'digit mantissa.
    '
    '2 References:
    'Hildebrand, F. B. Introduction to Numerical Analysis.
    'McGraw-Hill, New York 1956, p. 155 ;
    'Henrici, Peter. Elements of Numerical Analysis.
    'Wiley, New York 1964, p. 252.
    'Person: https://en.wikipedia.org/wiki/Peter_Henrici_%28mathematician%28  : The Book:
    'https://ia700700.us.archive.org/23/items/ElementsOfNumericalAnalysis/Henrici-ElementsOfNumericalAnalysis.pdf
    '--------------------------------------------------------------------
    Declare i&,j&,h!,cj!,c![nmax&+1],b![nmax&]
    b![0]=1:b![n&]=0
    c![0]=1:c![1]=-0.5
    w![0]=0.5:w![n&]=w![0]
    h!=(t![n&]-t![0])/n&

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

        w![i&]=1
        b![i&]=0
        t![i&]=i&*h!+t![0]

    endwhile

    case r&>n&:r&=n&

    Whileloop r&:j&=&Loop

        cj!=0.5*c![j&]

        whileloop j&,1,-1:i&=&Loop

            b![i&]=b![i&]-b![i&-1]

        endwhile

        whileloop 3,j&+2:i&=&Loop

            cj!=cj!+c![j&+2-i&]/i&

        endwhile

        c![j&+1]=-cj!

        whileloop 0,n&:i&=&Loop

            w![i&]=w![i&]-cj!*(b![n&-i&]+b![i&])

        endwhile

    Endwhile

    WhileLoop 0,n&:i&=&Loop

        w![i&]=w![i&]*h!

    endWhile

EndProc'Gregory

'--------------------------------------------------------------------
'  Testteil
'--------------------------------------------------------------------
declare t![nmax&],w![nmax&]:declare I!,In!,i&,n&,r&,p&:declare a!,b!
a!=-1.0'Integration von a
b!= 1.0'bis b
n&= 7'Grad des impliziten Polynoms
r&=n&
t![0]=a!:t![n&]=b!' Integrationsgrenzen belegen
Gregory(n&,r&,t![],w![])' Polynom-Koeffizienten (Gregory-Gewichtungen) berechnen
' Kontrollausgabe:
font 2:Print "\n       Gregory-Gewichtungen:\n "+mkstr$("-",40)

WhileLoop 0,n&:i&=&Loop

    print tab(2);if(w![i&]<0,""," ");format$("%g",w![i&]),
    print tab(27);if(t![i&]<0,""," ");format$("%g",t![i&])
    endwhile:print
    '-----------------------------------------------------------------------
    ' Testanwendung: Integriere die Funktion x^p und gib die Ergebnisse aus:
    '-----------------------------------------------------------------------

    whileloop 0,r&+4:p&=&Loop

        I!= ( Fn_PowN(b!,p&+1) - Fn_PowN(a!,p&+1) ) / (p&+1)
        In!=0

        whileloop 0,n&:i&=&Loop

            In!=In!+w![i&]*Fn_PowN(t![i&],p&)

        endwhile

        print tab(2);n&,tab(5);r&,tab(8);p&,
        print tab(13);if(I!<0,""," ");format$("%g",I!),
        print tab(43);if(In!<0,""," ");format$("%g",In!),
        print tab(73);if((I!-In!)<0,""," ");format$("%g",I!-In!)

    Endwhile

    BEEP
    Print "\n Es folgen die Vergleichsdaten lt. Buch [Taste]\n"+mkstr$("-",80)
    waitinput
    var wdata$=\
    "  0.086921 -1.000000 \n  0.414005 -0.714286 \n"+\
    "  0.153125 -0.428571 \n  0.345949 -0.142857 \n"+\
    "  0.345949  0.142857 \n  0.153125  0.428571 \n"+\
    "  0.414005  0.714286 \n  0.086921  1.000000 \n"
    var outdata$= "\n "+\
    "7 7  0  2.000000  2.000000  4.440892e-16 \n "+\
    "7 7  1  0.000000 -0.000000  1.110223e-16 \n "+\
    "7 7  2  0.666667  0.666667  2.220446e-16 \n "+\
    "7 7  3  0.000000 -0.000000  8.326673e-17 \n "+\
    "7 7  4  0.400000  0.400000  1.665335e-16 \n "+\
    "7 7  5  0.000000 -0.000000  8.326673e-17 \n "+\
    "7 7  6  0.285714  0.285714  0.000000e+00 \n "+\
    "7 7  7  0.000000 -0.000000  4.163336e-17 \n "+\
    "7 7  8  0.222222  0.230297 -8.075245e-03 \n "+\
    "7 7  9  0.000000 -0.000000  2.775558e-17 \n "+\
    "7 7 10  0.181818  0.202532 -2.071405e-02 \n "+\
    "7 7 11  0.000000 -0.000000  1.387779e-17"
    font 2:print Wdata$,outdata$
    waitinput
    End