Das sieht auf den ersten Blick auch gar nicht mal schlecht aus… wenn ich nur wüsste, was O(m^5) bedeutet.
Kovergieren tut da, soweit ich das sehe, überhaupt nichts, zumal das (möglicherweise bis auf besagtes O()) alles in einer Zeile geht. Wie genau ist denn das?
Anscheinend brauch ich Wolfram Pro um das Ganze dort zu finden. Aber Registrieren oder irgend etwas erwerben wollte ich für dieses kleine Problem eigentlich nicht.
Ja, das hilft schon mal weiter. Die Genauigkeit hängt dann wohl davon ab, wieviele Glieder man addiert. Für die Faktoren im Nenner braucht man auf OEIS nur nach der Sequenz 8, 128, 512 und 32768 suchen. Ich werde da mal mit rumexperimentieren.
Danke
Tja schade eigentlich…
Nach 20 Gliedern unterscheidet sich ein Grenzfall (0.9999) von dem ursprünglichen Wert 1.0005145000840605 noch unheimlich von 1.0124863641844215
Für mehr als 16 Glieder reicht ein Double auch schon nicht mehr aus.
In der Hoffnung, dass 105 Glieder ausreichen, würde ich als nächsten Schritt aπ/b über BigDecimals zusammengefasst berechnen, damit ich ein einziges Array an Koeffizienten bekomme, die ich nur noch mit m^n (für 0<=n<104) multiplizieren muss. Leider kann ich Dank binomial() mit der Formel
a(n) = 2^(-2 w[n])binomial[2n, n]^2 (-1)^(2n)/(1-2n)
mal so gar nichts anfangen.
Edit: Die 105 Koeffizienten für die Nenner waren schon mal kein Problem.
Edit2: Die 105 Koeffizienten für die Zähler funktionieren nun auch, aber leider leider konvergiert auch dieses immer langsamer, je größer n wird. Nach den 105 Gliedern ergibt k=0,9996 1.0035864106410735 obwohl 1.0017808456277624 heraus kommen sollte.
Optimal wäre ein Algo, der auf den Bereich 0,9996<k<1.0 optimiert ist. k=1.0 ergibt 1.0 und k>1.0 ergibt -Infinity.
Das ist nur für das normale Integral. Beim vollständigen bleibt der Parameter x unangetastet auf PI/2.
Ich denke mal, es läuft darauf hinaus, dass ich bei Werten zwischen 0,9996 und 1 auf die Annäherung durch ein Polynom zurückgreifen muss.
public static double completeEllipticIntegralSecondKind(double k) {
if (k > 1.0) {
return Double.NEGATIVE_INFINITY;
}
if (k == 1.0) {
return 1.0;
}
double k1 = k;
k *= k;
double k2 = sqrt(1 - k);
if(k1 > k2) {
double kk1 = completeEllipticIntegralFirstKind(k1);
double kk2 = completeEllipticIntegralFirstKind(k2);
double ek1 = completeEllipticIntegralSecondKind(k2);
double ek2 = (PI_HALFED + kk1 * kk2 - kk1 * ek1) / kk2;
return ek2;
}
double sum = 1.0;
int n = 1;
double term = 1.0;
int above = 1;
int below = 2;
double sumOld = 0.0;
while (sum != sumOld && n > 0) {
sumOld = sum;
term *= above / (double) below;
sum -= pow(k, n) * pow(term, 2) / above;
above += 2;
below += 2;
n++;
}
sum *= PI_HALFED;
return sum;
}
public static double completeEllipticIntegralFirstKind(double k) {
return PI_HALFED / agm(1, sqrt(1 - k * k));
}
public static double agm(double a, double g) {
if(a == g) {
return a;
}
if(g == 0) {
return 0;
}
double a1 = 0.5 * (a + g);
double g1 = sqrt(a * g);
if(a == a1 && g == g1) {
return (a + g) / 2.0;
}
return agm(a1, g1);
}
Das Intergal 1. Art geht dank dem Arithmetisch-geometrisches Mittel immer recht schnell und dank der Legendres Relation, braucht man nur noch abzuwägen, wann k > √(1-k²) ist. Dann berechnet man das kleinere Intergral 2. Art und berechnet daraus das Größere. Das ist zwar auf den letzten beiden Nachkommastellen ungenau, aber was solls. Ich will gar nicht wissen, wie ungenau die polynimiale Annäherung gewesen wär.
