XProfan

WindowTitle "ROMBERG-INTEGRATION"
'{Numerische Integration einer Funktionskurve
' nach dem Rundungsfehler-resistenten Romberg-Verfahren
' (Siehe https://de.wikipedia.org/wiki/Werner_Romberg)
' Konvergiert bei allen stetigen Funktionen!
' Demo-Übersetzung nach XProfan 11.2a 2012-05 by P.Specht,
' Quelle: https://de.wikipedia.org/wiki/Romberg-Integration,
' plus zwei 30 Jahre alte EDV-Zeitschriften. Ohne Gewähr!
' Rechtslage ungeprüft! Nur per Demonstrationszwecke!
'}
'{ Initialisierung
Font 2:randomize:cls rnd(8^8):set("decimals",18)
AppendMenuBar 10,"Numerische Integration nach Rhomberg"
Declare xu!,xo!,n!,g!,tmp!,Integral!
'}

proc Fnk :parameters x!

    var a!=100
    var b!=50
    declare y!,aa!,er%
    aa!=a!*a!

    if (a!<>0) and (a!<>x!)

        ' ======= PROGRAMMTEIL A =============
        y! = b!*sqrt(1-x!*x!/aa!)' ELLIPSE a , b
        ' ====================================

    else

        y!=-1*10^-35

    endif

    return y!

endproc

'{ Ein/Ausgabeteil
Begin:
print "\n Die Funktion ist in Programmteil A programmiert!"
print "\n Untere Integrationsgrenze Xu= ";:input xu!
print "\n  Obere Integrationsgrenze Xo= ";:input xo!
print "\n Abbruch-Genauigkeit [Stellen]: ";:input g!
g!=1/10^g!
Integral! = Romberg(xu!,xo!,g!)
print "\n Ergebnis:\n"
print " Das Flächen-Integral zwischen "
print " ";xu!;" und ";xo!
print " beträgt ";format$("%e",Integral!)
print " mit Fehlerschranke ";format$("%e",tmp!)
WaitInput
Goto "Begin"
'}

proc Romberg : parameters xu!,xo!

    var anz&=10' Streifenzahl, empfohlen 10, 12, max.16 (Laufzeit quadratisch!)
    Declare i&,j&,k&,n&[anz&+1],H![anz&+1],L![anz&,anz&],Q!
    n&[0]=2
    H![0]=(xo!-xu!)/n&[0]
    ' benutzt Trapezregel:
    L![0,0]=H![0]/2*(Fnk(xu!)+Fnk(xo!)+2*Fnk(xu!+H![0]))

    WhileLoop Anz&:j&=&Loop

        H![j&]=H![0]/(2^j&)
        n&[j&]=n&[0]*(2^j&)
        Q!=0

        whileLoop 0,n&[j&-1]-1:i&=&Loop

            Q!=Q! + Fnk(xu!+(2*i&+1)*H![j&])

        endwhile

        L![0,j&]=L![0,j&-1]/2+H![j&]*Q!

    EndWhile

    WhileLoop Anz&:k&=&Loop

        whileloop 0,Anz&-1:j&=&Loop

            L![k&,j&]=1/(2^(2*k&)-1)*(2^(2*k&)*L![k&-1,j&+1]-L![k&-1,j&])

        endwhile

        tmp!=abs(L![k&,0]-L![k&-1,1])
        case tmp!<=g!:break

    Endwhile

    ' tmp! enthält aktuelle Fehlergrenze
    return L![k&,0]

endproc