diff options
author | Nicolas Tobler <nicolas.tobler@ost.ch> | 2022-08-03 20:37:12 +0200 |
---|---|---|
committer | Nicolas Tobler <nicolas.tobler@ost.ch> | 2022-08-03 20:37:12 +0200 |
commit | e08392d4bacb9a54c3ab755fa6445514749b608f (patch) | |
tree | 67af5a4a6ed541b1b425de89fd05c2a74a265571 /buch/chapters/110-elliptisch/agm | |
parent | improved Einleitung (diff) | |
parent | Merge pull request #39 from NaoPross/master (diff) | |
download | SeminarSpezielleFunktionen-e08392d4bacb9a54c3ab755fa6445514749b608f.tar.gz SeminarSpezielleFunktionen-e08392d4bacb9a54c3ab755fa6445514749b608f.zip |
Merge branch 'master' of https://github.com/AndreasFMueller/SeminarSpezielleFunktionen
Diffstat (limited to '')
-rw-r--r-- | buch/chapters/110-elliptisch/agm/Makefile | 15 | ||||
-rw-r--r-- | buch/chapters/110-elliptisch/agm/agm.cpp | 75 | ||||
-rw-r--r-- | buch/chapters/110-elliptisch/agm/agm.m | 20 | ||||
-rw-r--r-- | buch/chapters/110-elliptisch/agm/agm.maxima | 26 | ||||
-rw-r--r-- | buch/chapters/110-elliptisch/agm/sn.cpp | 52 |
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; +} |