aboutsummaryrefslogtreecommitdiffstats
path: root/buch/chapters/110-elliptisch/agm
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--buch/chapters/110-elliptisch/agm/Makefile15
-rw-r--r--buch/chapters/110-elliptisch/agm/agm.cpp75
-rw-r--r--buch/chapters/110-elliptisch/agm/agm.m20
-rw-r--r--buch/chapters/110-elliptisch/agm/agm.maxima26
-rw-r--r--buch/chapters/110-elliptisch/agm/sn.cpp52
5 files changed, 188 insertions, 0 deletions
diff --git a/buch/chapters/110-elliptisch/agm/Makefile b/buch/chapters/110-elliptisch/agm/Makefile
new file mode 100644
index 0000000..8dab511
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/Makefile
@@ -0,0 +1,15 @@
+#
+# Makefile
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+all: sn
+
+sn: sn.cpp
+ g++ -O -Wall -g -std=c++11 sn.cpp -o sn `pkg-config --cflags gsl` `pkg-config --libs gsl`
+
+
+agm: agm.cpp
+ g++ -O -Wall -g -std=c++11 agm.cpp -o agm `pkg-config --cflags gsl` `pkg-config --libs gsl`
+ ./agm
+
diff --git a/buch/chapters/110-elliptisch/agm/agm.cpp b/buch/chapters/110-elliptisch/agm/agm.cpp
new file mode 100644
index 0000000..8abb4b2
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/agm.cpp
@@ -0,0 +1,75 @@
+/*
+ * agm.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+#include <iostream>
+#include <gsl/gsl_sf_ellint.h>
+
+inline long double sqrl(long double x) {
+ return x * x;
+}
+
+long double Xn(long double a, long double b, long double x) {
+ double long epsilon = fabsl(a - b);
+ if (epsilon > 0.001) {
+ return (a - sqrtl(sqrl(a) - sqrl(x) * (a + b) * (a - b)))
+ / (x * (a - b));
+ }
+ long double d = a + b;
+ long double x1 = 0;
+ long double y2 = sqrl(x/a);
+ long double c = 1;
+ long double s = 0;
+ int k = 1;
+ while (c > 0.0000000000001) {
+ c *= (0.5 - (k - 1)) / k;
+ c *= (d - epsilon) * y2;
+ s += c;
+ c *= epsilon;
+ c = -c;
+ k++;
+ }
+ return s * a / x;
+}
+
+int main(int argc, char *argv[]) {
+ long double a = 1;
+ long double b = sqrtl(2.)/2;
+ long double x = 0.7;
+ if (argc >= 3) {
+ a = std::stod(argv[1]);
+ b = std::stod(argv[2]);
+ }
+ if (argc >= 4) {
+ x = std::stod(argv[3]);
+ }
+
+ {
+ long double an = a;
+ long double bn = b;
+ long double xn = x;
+ for (int i = 0; i < 10; i++) {
+ printf("%d %24.18Lf %24.18Lf %24.18Lf %24.18Lf\n",
+ i, an, bn, xn, a * asin(xn) / an);
+ long double A = (an + bn) / 2;
+ xn = Xn(an, bn, xn);
+ bn = sqrtl(an * bn);
+ an = A;
+ }
+ }
+
+ {
+ double k = b/a;
+ k = sqrt(1 - k*k);
+ double K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+ printf(" %24.18f %24.18f\n", k, K);
+ double F = gsl_sf_ellint_F(asinl(x), k, GSL_PREC_DOUBLE);
+ printf(" %24.18f %24.18f\n", k, F);
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/buch/chapters/110-elliptisch/agm/agm.m b/buch/chapters/110-elliptisch/agm/agm.m
new file mode 100644
index 0000000..dcb3ad8
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/agm.m
@@ -0,0 +1,20 @@
+#
+# agm.m
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+format long
+
+n = 10;
+a = 1;
+b = sqrt(0.5);
+
+for i = (1:n)
+ printf("%20.16f %20.16f\n", a, b);
+ A = (a+b)/2;
+ b = sqrt(a*b);
+ a = A;
+end
+
+K = pi / (2 * a)
+
diff --git a/buch/chapters/110-elliptisch/agm/agm.maxima b/buch/chapters/110-elliptisch/agm/agm.maxima
new file mode 100644
index 0000000..c7facd4
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/agm.maxima
@@ -0,0 +1,26 @@
+/*
+ * agm.maxima
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+
+S: 2*a*sin(theta1) / (a+b+(a-b)*sin(theta1)^2);
+
+C2: ratsimp(diff(S, theta1)^2 / (1 - S^2));
+C2: ratsimp(subst(sqrt(1-sin(theta1)^2), cos(theta1), C2));
+C2: ratsimp(subst(S, sin(theta), C2));
+C2: ratsimp(subst(sqrt(1-S^2), cos(theta), C2));
+
+D2: (a^2 * cos(theta)^2 + b^2 * sin(theta)^2)
+ /
+ (a1^2 * cos(theta1)^2 + b1^2 * sin(theta1)^2);
+D2: subst((a+b)/2, a1, D2);
+D2: subst(sqrt(a*b), b1, D2);
+D2: ratsimp(subst(1-S^2, cos(theta)^2, D2));
+D2: ratsimp(subst(S, sin(theta), D2));
+D2: ratsimp(subst(1-sin(theta1)^2, cos(theta1)^2, D2));
+
+Q: D2/C2;
+Q: ratsimp(subst(x, sin(theta1), Q));
+
+Q: ratsimp(expand(ratsimp(Q)));
diff --git a/buch/chapters/110-elliptisch/agm/sn.cpp b/buch/chapters/110-elliptisch/agm/sn.cpp
new file mode 100644
index 0000000..9e1b047
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/sn.cpp
@@ -0,0 +1,52 @@
+/*
+ * sn.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+#include <iostream>
+#include <gsl/gsl_sf_ellint.h>
+#include <gsl/gsl_sf_elljac.h>
+
+static const int N = 10;
+
+inline long double sqrl(long double x) {
+ return x * x;
+}
+
+int main(int argc, char *argv[]) {
+ long double u = 0.6;
+ long double k = 0.9;
+ long double kprime = sqrt(1 - sqrl(k));
+
+ long double a[N], b[N], x[N+1];
+ a[0] = 1;
+ b[0] = kprime;
+
+ for (int n = 0; n < N-1; n++) {
+ printf("a[%d] = %22.18Lf b[%d] = %22.18Lf\n", n, a[n], n, b[n]);
+ a[n+1] = (a[n] + b[n]) / 2;
+ b[n+1] = sqrtl(a[n] * b[n]);
+ }
+
+ x[N] = sinl(u * a[N-1]);
+ printf("x[%d] = %22.18Lf\n", N, x[N]);
+
+ for (int n = N - 1; n >= 0; n--) {
+ x[n] = 2 * a[n] * x[n+1] / (a[n] + b[n] + (a[n] - b[n]) * sqrl(x[n+1]));
+ printf("x[%2d] = %22.18Lf\n", n, x[n]);
+ }
+
+ printf("sn(%7.4Lf, %7.4Lf) = %20.24Lf\n", u, k, x[0]);
+
+ double sn, cn, dn;
+ double m = sqrl(k);
+ gsl_sf_elljac_e((double)u, m, &sn, &cn, &dn);
+ printf("sn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, sn);
+ printf("cn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, cn);
+ printf("dn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, dn);
+
+ return EXIT_SUCCESS;
+}