XProfan

Window Title "ROMBERG-INTEGRATION"
'{Numerische Integration of/ one Funktionskurve
' to the Rundungsfehler-resistenten Romberg-take action
' (See https://de.wikipedia.org/wiki/Werner_Romberg)
' Konvergiert with all stetigen functions!
' demonstration-Translation to XProfan 11.2a 2012-05 by P.woodpecker,
' fountain: https://de.wikipedia.org/wiki/Romberg-Integration,
' plus two 30 year old EDV-Zeitschriften. without Gewähr!
' legal situation ungeprüft! only for Demonstrationszwecke!
'}
'{ Initialisierung
Font 2:randomize:cls rnd(8^8):set("decimals",18)
AppendMenuBar 10,"Numerische Integration to Rhomberg"
Declare xu!,xo!,n!,g!,tmp!,Integral!
'}

proc Fnk :parameters x!

    var a!=100
    var b!=50
    declare y!,aa!,it%
    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

'{ One/Ausgabeteil
Begin:
print "\n The function is in Programmteil A program!"
print "\n downstairs Integrationsgrenze Xu= ";:input xu!
print "\n  Obere Integrationsgrenze Xo= ";:input xo!
print "\n discontinue-accuracy [to put]: ";:input g!
g!=1/10^g!
Integral! = Romberg(xu!,xo!,g!)
print "\n Result:\n"
print " the Flächen-Integral between "
print " ";xu!;"And ";xo!
print " totals ";stature$("%e",Integral!)
print " with Fehlerschranke ";stature$("%e",tmp!)
WaitInput
Goto "Begin"
'}

proc Romberg : parameters xu!,xo!

    var anz&=10' Streifenzahl, recommended 10, 12, max.16 (Laufzeit quadratic!)
    Declare i&,j&,k&,n&[anz&+1],H![anz&+1],L![anz&,anz&],Q!
    n&[0]=2
    H![0]=(xo!-xu!)/n&[0]
    ' using 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! contains actually Fehlergrenze
    return L![k&,0]

endproc