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.m | 20 ------------- buch/chapters/110-elliptisch/agm/Makefile | 10 +++++++ buch/chapters/110-elliptisch/agm/agm.cpp | 42 ++++++++++++++++++++++++++++ buch/chapters/110-elliptisch/agm/agm.m | 20 +++++++++++++ buch/chapters/110-elliptisch/ellintegral.tex | 34 ++++++++++++---------- 5 files changed, 91 insertions(+), 35 deletions(-) delete mode 100644 buch/chapters/110-elliptisch/agm.m 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') diff --git a/buch/chapters/110-elliptisch/agm.m b/buch/chapters/110-elliptisch/agm.m deleted file mode 100644 index 2f0a1ea..0000000 --- a/buch/chapters/110-elliptisch/agm.m +++ /dev/null @@ -1,20 +0,0 @@ -# -# 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 - -E = 2 / (pi * a) - 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) + diff --git a/buch/chapters/110-elliptisch/ellintegral.tex b/buch/chapters/110-elliptisch/ellintegral.tex index 970d8fa..79ed91e 100644 --- a/buch/chapters/110-elliptisch/ellintegral.tex +++ b/buch/chapters/110-elliptisch/ellintegral.tex @@ -578,9 +578,9 @@ Gauss hat gefunden, dass die Substitution \end{equation} zu \begin{equation} -\frac{dt}{a^2\cos^2 t + b^2 \sin^2 t} +\frac{dt}{\sqrt{a^2_{\phantom{1}}\cos^2 t + b^2_{\phantom{1}} \sin^2 t}} = -\frac{dt_1}{a_1^2\cos^2 t_1 + b_1^2 \sin^2 t_1} +\frac{dt_1}{\sqrt{a_1^2\cos^2 t_1 + b_1^2 \sin^2 t_1}} \label{buch:elliptisch:agm:dtdt1} \end{equation} führt. @@ -608,7 +608,7 @@ ableiten, es ist \frac{dt}{dt_1} \biggr)^2 = -\frac{a^2 \cos^2 t + b^2 \sin^2 t}{a_1^2 \cos^2 t_1 + b_1^2 \sin^2 t_1}. +\frac{a^2_{\phantom{1}} \cos^2 t + b^2_{\phantom{1}} \sin^2 t}{a_1^2 \cos^2 t_1 + b_1^2 \sin^2 t_1}. \] Man muss also nachprüfen, dass \begin{equation} @@ -618,7 +618,7 @@ Man muss also nachprüfen, dass \frac{a^2 \cos^2 t + b^2 \sin^2 t}{a_1^2 \cos^2 t_1 + b_1^2 \sin^2 t_1}. \label{buch:elliptisch:agm:deq} \end{equation} -Dazu muss man zunächst $a_1=(a+b)/2$, $b_1=\sqrt{ab}$ setzen. +Dazu muss man zunächst $a_1=(a+b)/2$, $b_1=\!\sqrt{ab}$ setzen. Ausserdem muss man $\cos^2 t$ durch $1-\sin^2t$ ersetzen und $\sin t$ durch \eqref{buch:elliptisch:agm:subst}. Auch $\cos^2 t_1$ muss man durch $1-\sin^2t_1$ ersetzt werden. @@ -724,15 +724,19 @@ K(k) = I(1,\sqrt{1-k^2}) = \frac{\pi}{2M(1,\sqrt{1-k^2})} \subsubsection{Numerisches Beispiel} \begin{table} \centering -\begin{tabular}{|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|} +\begin{tabular}{|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|} \hline -n& a_n & b_n \\ +n& a_n & b_n & \pi/2a_n \mathstrut +\text{\vrule height12pt depth6pt width0pt}\\ \hline -0 & 1.0000000000000000 & 0.7071067811865476\\ -1 & 0.\underline{8}535533905932737 & 0.\underline{84}08964152537146\\ -2 & 0.\underline{8472}249029234942 & 0.\underline{8472}012667468916\\ -3 & 0.\underline{847213084}8351929 & 0.\underline{8472130847}527654\\ -4 & 0.\underline{847213084793979}2 & 0.\underline{847213084793979}1\\ +\text{\vrule height12pt depth0pt width0pt} + 0 & 1.0000000000000000000 & 0.7071067811865475243 & 1.5707963267948965579 \\ + 1 & 0.8535533905932737621 & 0.8408964152537145430 & 1.\underline{8}403023690212201581 \\ + 2 & 0.8472249029234941526 & 0.8472012667468914603 & 1.\underline{8540}488143993356315 \\ + 3 & 0.8472130848351928064 & 0.8472130847527653666 & 1.\underline{854074677}2111781089 \\ + 4 & 0.8472130847939790865 & 0.8472130847939790865 & 1.\underline{854074677301371}8463 \\ +\infty& & & 1.8540746773013719184  +\text{\vrule height12pt depth6pt width0pt}\\ \hline \end{tabular} \caption{Die Berechnung des arithmetisch-geometrischen Mittels für @@ -747,11 +751,11 @@ Konvergenz der Berechnung des arithmetisch-geometrischen Mittels von $1$ und $\sqrt{2}/2$. Mit Satz~\ref{buch:elliptisch:agm:satz:Ek} folgt jetzt \[ -K(\sqrt{2}/2) +K(\!\sqrt{2}/2) = -\frac{\pi}{2M(1,\sqrt{2}/2)} +\frac{\pi}{2M(1,\!\sqrt{2}/2)} = -0.751428163461842. +1.854074677301372. \] Die Berechnung hat nur 4 Mittelwerte, 4 Produkte, 4 Quadratwurzeln und eine Division erfordert. @@ -848,7 +852,7 @@ Integrale vorkommen, haben Nullstellen bei $\pm1$, $\pm1/k$ und $\pm 1/\sqrt{n}$ % XXX Additionstheoreme \\ -XXX Parameterkonventionen \\ +% XXX Parameterkonventionen \\ % % Wertebereich -- cgit v1.2.1