XProfan

Italia

p.specht

Wenn´s darum geht, die Länge eines Bogenstückes der Ellipse zwischen zwei durch Zentriwinkel definierten Punkten zu ermitteln, gibt es mangels einfacher analytischer Formel (beim Kreis z.B. L = r * alpha) nur iterative Ansätze. Weil numerische Integration deutlich schneller geht als die sonst übliche MacLaurin-Reihenentwicklung mit zusätzlicher Integralrekursion, bemühte ich Prof. Werner Rombergs Integrationsalgorithmus.

Kompilieren Markieren Separieren

WindowTitle "Ellipsenbogenlänge via Romberg-Wegintegral"
'{ Initialisierung
' (D) Demoware 2012-06 by P. Specht. Ohne jedwede Gewähr!
Font 2:randomize:set("decimals",18)
Declare a!,b!,xu!,xo!,n!,g!,tmp!,Integral!,lambda2!,epsi!,epsi2!
var f!=pi()/180
'}

proc Fnk :parameters x!

    declare y!,er%,sq!
    sq!=sin(x!):sq!=sq!*sq!

    if (epsi2!*sq!)<=1

        ' ======= PROGRAMMTEIL A =============
        y! = a!*sqrt(1 - epsi2!*sq!)
        ' ====================================

    else

        y!=-1*10^-35

    endif

    return y!

endproc

'{ Ein/Ausgabeteil
Begin:
cls rnd(8^8):
print "\n Ellipsenbogenlänge: "
print "\n Grosse HALBachse a = "; :input a!
case a!=0 : goto "Begin"
print "\n Kleine Halbachse b = "; :input b!
case b!=0 : goto "Begin"
case a!=b!:print " Offenbar ein Kreis als Test! "
epsi2!=1-(b!*b!)/(a!*a!)
epsi!=sqrt(epsi2!)
print "\n Numer. Exzentrizität Epsilon = "; epsi!
print "\n Von Zentriwinkel alpha_1[° ] = ";:input xu!
print "\n   bis Zentriwinkel alpha_2[° ] =";:input xo!
print "\n Abbruch-Genauigkeit [Stellen]: ";:input g!
case g!=0:g!=10
g!=1/10^g!
Integral! = Rhomberg(xu!*f!,xo!*f!,g!)
print "\n Ergebnis:\n"
print " Die Bogenlänge zwischen "
print " ";xu!;" und ";xo!; " [°Grad]"
print " beträgt ";Integral!;"  bzw.  "; format$("%e",Integral!)
print " mit einem Fehler <= ";format$("%e",tmp!)

if (xu!=0) and (xo!=90)

    print "\n Kontrolle per Winkel 0° - 90° mit 1/4 der Ramanujan-UmfangsFormel: "
    lambda2!=(a!-b!)/(a!+b!):lambda2!=lambda2!*lambda2!
    print  " ";1/4*Pi()*(a!+b!)*(1+3*lambda2! / (10+sqrt(4-3*lambda2!) ) )
    print " mit einem rel. Fehler kleiner ";
    case (epsi! >= 0) and (epsi! < 0.8820): print "10^-9"
    case (epsi! >= 0.8820) and (epsi! < 0.9242): print "10^-8"
    case (epsi! >= 0.9242) and (epsi! < 0.9577): print "10^-7"
    case (epsi! >= 0.9577) and (epsi! < 0,9812): print "10^-6"
    case (epsi! >= 0.9812) and (epsi! < 0.9944): print "10^-5"
    case (epsi! >= 0.9944) and (epsi! < 0.9995): print "10^-4"
    case (epsi! >= 0.9995) and (epsi! < 0.9999): print "< 0.000403"
    case (epsi! >= 0.99999) and (epsi! <=1): print "< -0.5%"

endif

WaitInput
Goto "Begin"
'}

proc Rhomberg : parameters xu!,xo!

    var anz&=10' GERADE ZAHL!
    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