From 335c3a23f09759be380291ec89b0f2c43c2d3db6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Sat, 25 Jun 2022 22:46:16 +0200 Subject: fix agm --- buch/chapters/110-elliptisch/agm/Makefile | 10 ++++++++ buch/chapters/110-elliptisch/agm/agm.cpp | 42 +++++++++++++++++++++++++++++++ buch/chapters/110-elliptisch/agm/agm.m | 20 +++++++++++++++ 3 files changed, 72 insertions(+) create mode 100644 buch/chapters/110-elliptisch/agm/Makefile create mode 100644 buch/chapters/110-elliptisch/agm/agm.cpp create mode 100644 buch/chapters/110-elliptisch/agm/agm.m (limited to 'buch/chapters/110-elliptisch/agm') diff --git a/buch/chapters/110-elliptisch/agm/Makefile b/buch/chapters/110-elliptisch/agm/Makefile new file mode 100644 index 0000000..e7975e1 --- /dev/null +++ b/buch/chapters/110-elliptisch/agm/Makefile @@ -0,0 +1,10 @@ +# +# Makefile +# +# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +# + +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..fdb0441 --- /dev/null +++ b/buch/chapters/110-elliptisch/agm/agm.cpp @@ -0,0 +1,42 @@ +/* + * agm.cpp + * + * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include +#include +#include +#include +#include + + + +int main(int argc, char *argv[]) { + long double a = 1; + long double b = sqrtl(2.)/2; + if (argc >= 3) { + a = std::stod(argv[1]); + b = std::stod(argv[2]); + } + + { + long double an = a; + long double bn = b; + for (int i = 0; i < 10; i++) { + printf("%d %24.18Lf %24.18Lf %24.18Lf\n", + i, an, bn, a * M_PI / (2 * an)); + long double A = (an + bn) / 2; + 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); + } + + 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) + -- cgit v1.2.1 From 05d75b0f467b2535db538ecaee461cf0c8b637d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Mon, 27 Jun 2022 20:17:16 +0200 Subject: add stuff for elliptic filters --- buch/chapters/110-elliptisch/agm/Makefile | 5 +++ buch/chapters/110-elliptisch/agm/agm.cpp | 37 ++++++++++++++++++-- buch/chapters/110-elliptisch/agm/agm.maxima | 26 +++++++++++++++ buch/chapters/110-elliptisch/agm/sn.cpp | 52 +++++++++++++++++++++++++++++ 4 files changed, 118 insertions(+), 2 deletions(-) create mode 100644 buch/chapters/110-elliptisch/agm/agm.maxima create mode 100644 buch/chapters/110-elliptisch/agm/sn.cpp (limited to 'buch/chapters/110-elliptisch/agm') diff --git a/buch/chapters/110-elliptisch/agm/Makefile b/buch/chapters/110-elliptisch/agm/Makefile index e7975e1..8dab511 100644 --- a/buch/chapters/110-elliptisch/agm/Makefile +++ b/buch/chapters/110-elliptisch/agm/Makefile @@ -3,6 +3,11 @@ # # (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` diff --git a/buch/chapters/110-elliptisch/agm/agm.cpp b/buch/chapters/110-elliptisch/agm/agm.cpp index fdb0441..8abb4b2 100644 --- a/buch/chapters/110-elliptisch/agm/agm.cpp +++ b/buch/chapters/110-elliptisch/agm/agm.cpp @@ -9,23 +9,54 @@ #include #include +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\n", - i, an, bn, a * M_PI / (2 * an)); + 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; } @@ -36,6 +67,8 @@ int main(int argc, char *argv[]) { 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.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..ff2ed17 --- /dev/null +++ b/buch/chapters/110-elliptisch/agm/sn.cpp @@ -0,0 +1,52 @@ +/* + * ns.cpp + * + * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include +#include +#include +#include +#include +#include + +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; +} -- cgit v1.2.1 From 6c3d3784c6d03fd89450bcf05b60cdf888d23333 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Tue, 28 Jun 2022 18:14:35 +0200 Subject: typos --- buch/chapters/110-elliptisch/agm/sn.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'buch/chapters/110-elliptisch/agm') diff --git a/buch/chapters/110-elliptisch/agm/sn.cpp b/buch/chapters/110-elliptisch/agm/sn.cpp index ff2ed17..9e1b047 100644 --- a/buch/chapters/110-elliptisch/agm/sn.cpp +++ b/buch/chapters/110-elliptisch/agm/sn.cpp @@ -1,5 +1,5 @@ /* - * ns.cpp + * sn.cpp * * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule */ -- cgit v1.2.1