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/Makefile.inc | 3 + 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 ++++++ buch/chapters/110-elliptisch/chapter.tex | 7 +- buch/chapters/110-elliptisch/dglsol.tex | 89 +++++++++++ buch/chapters/110-elliptisch/ellintegral.tex | 141 +++++++++++++++- buch/chapters/110-elliptisch/elltrigo.tex | 2 +- buch/chapters/110-elliptisch/experiments/KK.pdf | Bin 0 -> 23248 bytes buch/chapters/110-elliptisch/experiments/KK.tex | 66 ++++++++ buch/chapters/110-elliptisch/experiments/KN.cpp | 177 +++++++++++++++++++++ buch/chapters/110-elliptisch/experiments/Makefile | 15 ++ .../chapters/110-elliptisch/experiments/agm.maxima | 26 --- buch/chapters/110-elliptisch/uebungsaufgaben/2.tex | 61 +++++++ buch/chapters/110-elliptisch/uebungsaufgaben/3.tex | 135 ++++++++++++++++ .../110-elliptisch/uebungsaufgaben/landen.m | 60 +++++++ buch/chapters/references.bib | 9 ++ 18 files changed, 877 insertions(+), 34 deletions(-) create mode 100644 buch/chapters/110-elliptisch/agm/agm.maxima create mode 100644 buch/chapters/110-elliptisch/agm/sn.cpp create mode 100644 buch/chapters/110-elliptisch/experiments/KK.pdf create mode 100644 buch/chapters/110-elliptisch/experiments/KK.tex create mode 100644 buch/chapters/110-elliptisch/experiments/KN.cpp create mode 100644 buch/chapters/110-elliptisch/experiments/Makefile delete mode 100644 buch/chapters/110-elliptisch/experiments/agm.maxima create mode 100644 buch/chapters/110-elliptisch/uebungsaufgaben/2.tex create mode 100644 buch/chapters/110-elliptisch/uebungsaufgaben/3.tex create mode 100644 buch/chapters/110-elliptisch/uebungsaufgaben/landen.m (limited to 'buch/chapters') diff --git a/buch/chapters/110-elliptisch/Makefile.inc b/buch/chapters/110-elliptisch/Makefile.inc index 639cb8f..ef6ea51 100644 --- a/buch/chapters/110-elliptisch/Makefile.inc +++ b/buch/chapters/110-elliptisch/Makefile.inc @@ -12,4 +12,7 @@ CHAPTERFILES += \ chapters/110-elliptisch/mathpendel.tex \ chapters/110-elliptisch/lemniskate.tex \ chapters/110-elliptisch/uebungsaufgaben/1.tex \ + chapters/110-elliptisch/uebungsaufgaben/2.tex \ + chapters/110-elliptisch/uebungsaufgaben/3.tex \ + chapters/110-elliptisch/uebungsaufgaben/4.tex \ chapters/110-elliptisch/chapter.tex 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; +} diff --git a/buch/chapters/110-elliptisch/chapter.tex b/buch/chapters/110-elliptisch/chapter.tex index e05f3bd..d65570b 100644 --- a/buch/chapters/110-elliptisch/chapter.tex +++ b/buch/chapters/110-elliptisch/chapter.tex @@ -35,11 +35,14 @@ wieder hergestellt. \input{chapters/110-elliptisch/lemniskate.tex} -\section*{Übungsaufgabe} -\rhead{Übungsaufgabe} +\section*{Übungsaufgaben} +\rhead{Übungsaufgaben} \aufgabetoplevel{chapters/110-elliptisch/uebungsaufgaben} \begin{uebungsaufgaben} %\uebungsaufgabe{0} \uebungsaufgabe{1} +\uebungsaufgabe{2} +\uebungsaufgabe{3} +\uebungsaufgabe{4} \end{uebungsaufgaben} diff --git a/buch/chapters/110-elliptisch/dglsol.tex b/buch/chapters/110-elliptisch/dglsol.tex index 3303aee..3709300 100644 --- a/buch/chapters/110-elliptisch/dglsol.tex +++ b/buch/chapters/110-elliptisch/dglsol.tex @@ -340,7 +340,96 @@ Die Jacobischen elliptischen Funktionen sind daher inverse Funktionen der unvollständigen elliptischen Integrale. % +% Numerische Berechnung mit dem arithmetisch-geometrischen Mittel % +\subsubsection{Numerische Berechnung mit dem arithmetisch-geometrischen Mittel} +\begin{table} +\centering +\begin{tikzpicture}[>=latex,thick] + +\begin{scope}[xshift=-2.4cm,yshift=1.2cm] +\fill[color=red!20] + (-1.0,0) -- (-1.0,-2.1) -- (-1.8,-2.1) -- (0,-3.0) + -- (1.8,-2.1) -- (1.0,-2.1) -- (1.0,0) -- cycle; +\node[color=white] at (0,-1.2) [scale=7] {\sf 1}; +\end{scope} + +\begin{scope}[xshift=2.9cm,yshift=-1.8cm] +\fill[color=blue!20] + (0.8,0) -- (0.8,2.1) -- (1.4,2.1) -- (0,3.0) -- (-1.4,2.1) + -- (-0.8,2.1) -- (-0.8,0) -- cycle; +\node[color=white] at (0,1.2) [scale=7] {\sf 2}; +\end{scope} + +\node at (0,0) { +\begin{tabular}{|>{$}c<{$}|>{$}c<{$}>{$}c<{$}|>{$}c<{$}>{$}l<{$}|} +\hline +n & a_n & b_n & x_n & +\mathstrut\text{\vrule height12pt depth6pt width0pt}\\ +\hline +0 & 1.0000000000000000 & 0.4358898943540673 & 0.5422823228691580 & = \operatorname{sn}(u,k)% +\mathstrut\text{\vrule height12pt depth0pt width0pt}\\ +1 & 0.7179449471770336 & 0.6602195804079634 & 0.4157689781689663 & \mathstrut\\ +2 & 0.6890822637924985 & 0.6884775317911533 & 0.4017521410983242 & \mathstrut\\ +3 & 0.6887798977918259 & 0.6887798314243237 & 0.4016042867931862 & \mathstrut\\ +4 & 0.6887798646080748 & 0.6887798646080740 & 0.4016042705654757 & \mathstrut\\ +5 & 0.6887798646080744 & 0.6887798646080744 & 0.4016042705654755 & \mathstrut\\ +6 & & & 0.4016042705654755 & = \sin(a_5u) +\mathstrut\text{\vrule height0pt depth6pt width0pt}\\ +\hline +\end{tabular} +}; +\end{tikzpicture} +\caption{Berechnung von $\operatorname{sn}(u,k)$ für $u=0.6$ und $k=0.$2 +mit Hilfe des arithmetisch-geo\-me\-tri\-schen Mittels. +In der ersten Phase des Algorithmus (rot) wird die Folge der arithmetischen +und geometrischen Mittel berechnet, in der zweiten Phase werden die +Approximationen von $x_0=\operatorname{sn}(u,k)$. +Bei $n=5$ erreicht die Iteration des arithmetisch-geometrischen Mittels +Maschinengenauigkeit, was sich auch darin äussert, dass sich $x_5$ und +$x_6=\sin(a_5u)$ nicht unterscheiden. +\label{buch:elliptisch:agm:table:snberechnung}} +\end{table} +In Abschnitt~\ref{buch:elliptisch:subsection:agm} auf +Seite~\pageref{buch:elliptisch:subsubection:berechnung-fxk-agm} +wurde erklärt, wie das unvollständige elliptische Integral $F(x,k)$ mit +Hilfe des arithmetisch-geometrischen Mittels berechnet werden kann. +Da $\operatorname{sn}^{-1}(x,k) = F(x,k)$ die Umkehrfunktion ist, kann +man den Algorithmus auch zur Berechnung von $\operatorname{sn}(u,k)$ +verwenden. +Dazu geht man wie folgt vor: +\begin{enumerate} +\item +$k'=\sqrt{1-k^2}$. +\item +Berechne die Folgen des arithmetisch-geometrischen Mittels +$a_n$ und $b_n$ mit $a_0=1$ und $b_0=k'$, bis zum Folgenindex $N$, +bei dem ausreichende Konvergenz eintegreten ist. +\item +Setze $x_N = \sin(a_N \cdot u)$. +\item +Berechnet für absteigende $n=N-1,N-2,\dots$ die Folge $x_n$ mit Hilfe +der Rekursionsformel +\begin{equation} +x_{n} += +\frac{2a_nx_{n+1}}{a_n+b_n+(a_n-b_n)x_{n+1}^2}, +\label{buch:elliptisch:agm:xnrek} +\end{equation} +die aus \eqref{buch:elliptisch:agm:subst} +durch die Substitution $x_n = \sin t_n$ entsteht. +\item +Setze $\operatorname{sn}(u,k) = x_0$. +\end{enumerate} +Da die Formel \eqref{buch:elliptisch:agm:xnrek} nicht unter den +numerischen Stabilitätsproblemen leidet, die früher auf +Seite~\pageref{buch:elliptisch:agm:ellintegral-stabilitaet} +diskutiert wurden, ist die Berechnung stabil und sehr schnell. +Tabelle~\ref{buch:elliptisch:agm:table:snberechnung} +zeigt die Berechnung am Beispiel $u=0.6$ und $k=0.2$. + +% +% Pole und Nullstellen der Jacobischen elliptischen Funktionen % \subsubsection{Pole und Nullstellen der Jacobischen elliptischen Funktionen} \begin{figure} diff --git a/buch/chapters/110-elliptisch/ellintegral.tex b/buch/chapters/110-elliptisch/ellintegral.tex index 4589ffa..cc99218 100644 --- a/buch/chapters/110-elliptisch/ellintegral.tex +++ b/buch/chapters/110-elliptisch/ellintegral.tex @@ -459,7 +459,8 @@ Parameter $k$ mit der Ableitungsformel für die Funktion $\mathstrut_2F_1$. % % Berechnung mit dem arithmetisch-geometrischen Mittel % -\subsection{Berechnung mit dem arithmetisch-geometrischen Mittel} +\subsection{Berechnung mit dem arithmetisch-geometrischen Mittel +\label{buch:elliptisch:subsection:agm}} Die numerische Berechnung von elliptischer Integrale mit gewöhnlichen numerischen Integrationsroutinen ist nicht sehr effizient. Das in diesem Abschnitt vorgestellte arithmetisch-geometrische Mittel @@ -472,7 +473,11 @@ Sie ist ein Speziallfall der sogenannten Landen-Transformation, \index{Landen-Transformation}% welche ausser für die elliptischen Integrale auch für die Jacobischen elliptischen Funktionen formuliert werden kann und -für letztere ebenfalls sehr schnelle numerische Algorithmen liefert. +für letztere ebenfalls sehr schnelle numerische Algorithmen liefert +(siehe dazu auch die +Aufgaben~\ref{buch:elliptisch:aufgabe:2}--\ref{buch:elliptisch:aufgabe:4}). +Sie kann auch verwendet werden, um die Werte der Jacobischen elliptischen +Funktionen für komplexe Argument zu berechnen. % % Das arithmetisch-geometrische Mittel @@ -574,7 +579,7 @@ Gauss hat gefunden, dass die Substitution \begin{equation} \sin t = -\frac{2a\sin t_1}{a+b+(a-b)\sin t_1} +\frac{2a\sin t_1}{a+b+(a-b)\sin^2 t_1} \label{buch:elliptisch:agm:subst} \end{equation} zu @@ -1103,6 +1108,136 @@ F(x,k) = iK(k') - F\biggl(\frac1{kx},k\biggr) für die Werte des elliptischen Integrals erster Art für grosse Argumentwerte fest. +% +% AGM und Berechnung von F(x,k) +% +\subsubsection{Berechnung von $F(x,k)$ mit dem arithmetisch-geometrischen Mittel\label{buch:elliptisch:subsubection:berechnung-fxk-agm}} +Wie das vollständige elliptische Integral $K(k)$ kann auch das +unvollständige elliptische Integral +\begin{align*} +F(x,k) +&= +\int_0^x \frac{d\xi}{\sqrt{(1-\xi^2)(1-k^{\prime 2}\xi^2)}} += +\int_0^{\varphi} +\frac{dt}{\sqrt{1-k^2 \sin^2 t}} +\\ +&= +a +\int_0^{\varphi} \frac{dt}{a^2 \cos^2 t + b^2 \sin^2 t} +\qquad\text{mit $k=b/a$} +\end{align*} +mit dem arithmetisch-geometrischen Mittel berechnet werden. +Dazu muss die Substitution +\eqref{buch:elliptisch:agm:subst} +verwendet werden, um auch den Winkel $\varphi_1$ zu berechnen. +Dazu muss \eqref{buch:elliptisch:agm:subst} nach $x_1=\sin t_1$ +aufgelöst werden. +Durch Multiplikation mit dem Nenner erhält man mit der Abkürzung +$x=\sin t$ und $x_1=\sin t_1$ die quadratische Gleichung +\[ +(a-b)x x_1 +- +2ax_1 +(a+b)x += +0, +\] +mit der Lösung +\begin{equation} +x_1 += +\frac{a-\sqrt{a^2-(a^2-b^2)x^2}}{(a-b)x}. +\label{buch:elliptisch:unvollstagm:xrek} +\end{equation} +Der Algorithmus zur Berechnung des arithmetisch-geometrischen Mittels +muss daher verallgemeinert werden zu +\begin{equation} +\left. +\begin{aligned} +a_{n+1} &= \frac{a_n+b_n}2, &\qquad a_0 &= a +\\ +b_{n+1} &= \sqrt{a_nb_n}, & b_0 &= b +\\ +x_{n+1} &= \frac{a_n-\sqrt{a_n^2-(a_n^2-b_n^2)x_n^2}}{(a_n-b_n)x_n}, & x_0 &= x +\end{aligned} +\quad +\right\} +\label{buch:elliptisch:unvollstagm:rek} +\end{equation} +Die Folge $x_n$ konvergiert gegen einen Wert $x_{\infty} = \lim_{n\to\infty} x_n$. +Der Wert des unvollständigen elliptischen Integrals ist dann der Grenzwert +\[ +F(x,k) += +\lim_{n\to\infty} +\frac{\arcsin x_n}{M(a_n,b_n)} += +\frac{\arcsin x_{\infty}}{M(1,\sqrt{1-k^2})}. +\] + +In dieser Form ist die Berechnung allerdings nicht praktisch durchführbar. +Das Problem ist, dass die Differenz $a_n-b_n$, die in +\eqref{buch:elliptisch:unvollstagm:rek} +im Nenner vorkommt, sehr schnell gegen Null geht. +Ausserdem ist die Quadratwurzel im Zähler fast gleich gross wie +$a_n$, was zu Auslöschung und damit ungenauen Resultaten führt. +\label{buch:elliptisch:agm:ellintegral-stabilitaet} + +Eine Möglichkeit, das Problem zu entschärfen, ist, die Rekursionsformel +nach $\varepsilon = a-b$ zu entwickeln. +Mit $a+b=2a+\varepsilon$ kann man $b$ aus der Formel elimineren und erhält +mit Hilfe der binomischen Reihe +\begin{align*} +x_1 +&= +\frac{a}{x\varepsilon} +\left(1-\sqrt{1-\varepsilon(2a-\varepsilon)x^2/a^2}\right) +\\ +&= +\frac{a}{x\varepsilon} +\biggl( +1-\sum_{k=0}^\infty +(-1)^k +\frac{(\frac12)_k}{k!} \varepsilon^k(2a-\varepsilon)^k\frac{x^{2k}}{a^{2k}} +\biggr) +\\ +&= +\sum_{k=1}^\infty +(-1)^{k-1} +\frac{(\frac12)_k}{k!} \varepsilon^{k-1}(2a-\varepsilon)^k\frac{x^{2k-1}}{a^{2k-1}} +\\ +&= +\frac{\frac12}{1!}(2a-\varepsilon)\frac{x}{a} +- +\frac{\frac12\cdot(\frac12-1)}{2!}\varepsilon(2a-\varepsilon)^2\frac{x^3}{a^3} ++ +\frac{\frac12\cdot(\frac12-1)(\frac12-2)}{3!}\varepsilon^2(2a-\varepsilon)^3\frac{x^5}{a^5} +- +\dots +\\ +&= +x\biggl(1-\frac{\varepsilon}{2a}\biggr) +\biggl( +1 +- +\frac{\frac12-1}{2!}\varepsilon(2a-\varepsilon)\frac{x^2}{a^2} ++ +\frac{(\frac12-1)(\frac12-2)}{3!}\varepsilon^2(2a-\varepsilon)^2\frac{x^4}{a^4} +- +\dots +\biggr) +\\ +&= +x\biggl(1-\frac{\varepsilon}{2a}\biggr) +\cdot +\mathstrut_2F_1\biggl( +\begin{matrix}-\frac12,1\\2\end{matrix};-\varepsilon(2a-\varepsilon)\frac{x^2}{a^2} +\biggr). +\end{align*} +Diese Form ist wesentlich besser, aber leider kann es bei der numerischen +Rechnung passieren, dass $\varepsilon < 0$ wird. + %\subsection{Potenzreihe} %XXX Potenzreihen \\ %XXX Als hypergeometrische Funktionen \url{https://www.youtube.com/watch?v=j0t1yWrvKmE} \\ diff --git a/buch/chapters/110-elliptisch/elltrigo.tex b/buch/chapters/110-elliptisch/elltrigo.tex index 583e00a..c67870f 100644 --- a/buch/chapters/110-elliptisch/elltrigo.tex +++ b/buch/chapters/110-elliptisch/elltrigo.tex @@ -169,7 +169,7 @@ x^2(k^2-1) + y^2 = 1. an einer Ellipse mit Halbachsen $a$ und $1$. \label{buch:elliptisch:fig:jacobidef}} \end{figure} -\subsubsection{Definition der elliptischen Funktionen} +\subsubsection{Definition der Jacobischen elliptischen Funktionen} Die elliptischen Funktionen für einen Punkt $P$ auf der Ellipse mit Modulus $k$ können jetzt als Verhältnisse der Koordinaten des Punktes definieren. Es stellt sich aber die Frage, was man als Argument verwenden soll. diff --git a/buch/chapters/110-elliptisch/experiments/KK.pdf b/buch/chapters/110-elliptisch/experiments/KK.pdf new file mode 100644 index 0000000..13a2739 Binary files /dev/null and b/buch/chapters/110-elliptisch/experiments/KK.pdf differ diff --git a/buch/chapters/110-elliptisch/experiments/KK.tex b/buch/chapters/110-elliptisch/experiments/KK.tex new file mode 100644 index 0000000..a3ae425 --- /dev/null +++ b/buch/chapters/110-elliptisch/experiments/KK.tex @@ -0,0 +1,66 @@ +% +% KK.tex -- template for standalon tikz images +% +% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +% +\documentclass[tikz]{standalone} +\usepackage{amsmath} +\usepackage{times} +\usepackage{txfonts} +\usepackage{pgfplots} +\usepackage{csvsimple} +\usetikzlibrary{arrows,intersections,math} +\begin{document} +\def\skala{1} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\def\dx{10} +\def\dy{3} +\input{KKpath.tex} + +\draw[->] (-0.1,0) -- (10.3,0) coordinate[label={$k$}]; +\draw[->] (0,-0.1) -- (0,{2*\dy+0.3}) coordinate[label={right:$y$}]; + +\node at (3,{1.2*\dy}) {$\displaystyle y = \frac{K(k)}{K(\!\sqrt{1-k^2})}$}; + +\begin{scope} +\clip (0,0) rectangle (10,{2*\dy}); +\draw[color=red,line width=1.4pt] \KKpath; +\end{scope} + +\draw[line width=0.2pt] (10,0) -- (10,{2*\dy}); + +\foreach \y in {0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0}{ + \draw (-0.05,{\y*\dy}) -- (0.05,{\y*\dy}); + \node at (0,{\y*\dy}) [left] {$\y\mathstrut$}; +} + +\foreach \k in {1,...,9}{ + \draw ({\k*\dx/10},-0.05) -- ({\k*\dx/10},0.05); + \node at ({\k*\dx/10},0) [below] {$0.\k\mathstrut$}; +} +\node at (0,0) [below] {$0\mathstrut$}; +\node at (10,0) [below] {$1\mathstrut$}; + +\draw[color=blue] ({\knull*\dx},0) -- ({\knull*\dx},{\KKnull*\dy}); +\foreach \y in {1,2,3,4}{ + \draw[color=blue] + ({\knull*\dx-0.05},{\y*\KKnull*\dy/5}) + -- + ({\knull*\dx+0.05},{\y*\KKnull*\dy/5}); +} +\draw[color=black,line width=0.1pt] (0,{\KKnull*\dy}) -- ({\knull*\dx},{\KKnull*\dy}); +\draw[color=black,line width=0.1pt] (0,{\KKnull*\dy/5}) -- ({\kone*\dx},{\KKnull*\dy/5}); +\node at ({0.6*\dx},{\KKnull*\dy}) [above] {$y=1.7732$}; +\node at ({0.6*\dx},{\KKnull*\dy/5}) [above] {$y=0.3546$}; +\draw[color=blue] ({\kone*\dx},0) -- ({\kone*\dx},{\KKnull*\dy/5}); +\draw[color=blue] ({\kone*\dx},{\KKnull*\dy/5}) -- ({\knull*\dx},{\KKnull*\dy/5}); +\fill[color=blue] ({\kone*\dx},{\KKnull*\dy/5}) circle[radius=0.05]; +\fill[color=blue] ({\knull*\dx},{\KKnull*\dy/5}) circle[radius=0.05]; +\fill[color=blue] ({\knull*\dx},{\KKnull*\dy}) circle[radius=0.05]; +\node[color=blue] at ({\knull*\dx},0) [left,rotate=90] {$k=0.97\mathstrut$}; +\node[color=blue] at ({\kone*\dx},0) [left,rotate=90] {$k_1=0.0477$}; + +\end{tikzpicture} +\end{document} + diff --git a/buch/chapters/110-elliptisch/experiments/KN.cpp b/buch/chapters/110-elliptisch/experiments/KN.cpp new file mode 100644 index 0000000..1dcca9e --- /dev/null +++ b/buch/chapters/110-elliptisch/experiments/KN.cpp @@ -0,0 +1,177 @@ +/* + * KN.cpp + * + * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace KN { + +bool debug = false; + +static struct option longopts[] { +{ "debug", no_argument, NULL, 'd' }, +{ "N", required_argument, NULL, 'N' }, +{ "outfile", required_argument, NULL, 'o' }, +{ "min", required_argument, NULL, 'm' }, +{ NULL, 0, NULL, 0 } +}; + +double KprimeK(double k) { + double kprime = sqrt(1-k*k); + if (debug) + printf("%s:%d: k = %f, k' = %f\n", __FILE__, __LINE__, k, kprime); + double v + = + gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE) + / + gsl_sf_ellint_Kcomp(kprime, GSL_PREC_DOUBLE) + ; + if (debug) + printf("%s:%d: KprimeK(k = %f) = %f\n", __FILE__, __LINE__, k, v); + return v; +} + +static const int L = 100000000; +static const double h = 1. / L; + +double Kd(double k) { + double m = 0; + if (k < h) { + m = 2 * (KprimeK(k) - KprimeK(k / 2)) / k; + } else if (k > 1-h) { + m = 2 * (KprimeK((1 + k) / 2) - KprimeK(k)) / (1 - k); + + } else { + m = L * (KprimeK(k + h) - KprimeK(k)); + } + if (debug) + printf("%s:%d: Kd(%f) = %f\n", __FILE__, __LINE__, k, m); + return m; +} + +double k1(double y) { + if (debug) + printf("%s:%d: Newton for y = %f\n", __FILE__, __LINE__, y); + double kn = 0.5; + double delta = 1; + int n = 0; + while ((fabs(delta) > 0.000001) && (n < 10)) { + double yn = KprimeK(kn); + if (debug) + printf("%s:%d: k%d = %f, y%d = %f\n", __FILE__, __LINE__, n, kn, n, yn); + delta = (yn - y) / Kd(kn); + if (debug) + printf("%s:%d: delta = %f\n", __FILE__, __LINE__, delta); + double kneu = kn - delta; + if (kneu <= 0) { + kneu = kn / 4; + } + if (kneu >= 1) { + kneu = (3 + kn) / 4; + } + kn = kneu; + if (debug) + printf("%s:%d: kneu = %f, kn = %f\n", __FILE__, __LINE__, kneu, kn); + n++; + } + if (debug) + printf("%s:%d: Newton result: k = %f\n", __FILE__, __LINE__, kn); + return kn; +} + +double k1(int N, double k) { + return k1(KprimeK(k) / N); +} + +/** + * \brief Main function for the slcl program + */ +int main(int argc, char *argv[]) { + int longindex; + int c; + int N = 5; + double kmin = 0.01; + std::string outfilename; + while (EOF != (c = getopt_long(argc, argv, "d:N:o:m:", + longopts, &longindex))) + switch (c) { + case 'd': + debug = true; + break; + case 'N': + N = std::stoi(optarg); + break; + case 'o': + outfilename = std::string(optarg); + break; + case 'm': + kmin = std::stod(optarg); + break; + } + + double d = 0.01; + if (outfilename.size() > 0) { + FILE *fn = fopen(outfilename.c_str(), "w"); + fprintf(fn, "\\def\\KKpath{ "); + double k = d; + fprintf(fn, " (0,0)"); + double k0 = k/16; + while (k0 < k) { + fprintf(fn, "\n\t-- ({%.4f*\\dx},{%.4f*\\dy})", k0, KprimeK(k0)); + k0 *= 2; + } + while (k < 1-0.5*d) { + fprintf(fn, "\n\t-- ({%.4f*\\dx},{%.4f*\\dy})", k, KprimeK(k)); + k += d; + } + fprintf(fn, "}\n"); + + k0 = 0.97; + fprintf(fn, "\\def\\knull{%.4f}\n", k0); + double KK = KprimeK(k0); + fprintf(fn, "\\def\\KKnull{%.4f}\n", KK); + fprintf(fn, "\\def\\kone{%.4f}\n", k1(N, k0)); + + fclose(fn); + return EXIT_SUCCESS; + } + + for (double k = kmin; k < (1 - d/2); k += d) { + if (debug) + printf("%s:%d: k = %f\n", __FILE__, __LINE__, k); + double y = KprimeK(k); + double k0 = k1(y); + double kone = k1(N, k0); + printf("g(%4.2f) = %10.6f,", k, y); + printf(" g'(%.2f) = %10.6f,", k, Kd(k)); + printf(" g^{-1} = %10.6f,", k0); + printf(" k1 = %10.6f,", kone); + printf(" g(k1) = %10.6f\n", KprimeK(kone)); + } + + return EXIT_SUCCESS; +} + +} // namespace KN + +int main(int argc, char *argv[]) { + try { + return KN::main(argc, argv); + } catch (const std::exception& e) { + std::cerr << "terminated by exception: " << e.what(); + std::cerr << std::endl; + } catch (...) { + std::cerr << "terminated by unknown exception" << std::endl; + } + return EXIT_FAILURE; +} diff --git a/buch/chapters/110-elliptisch/experiments/Makefile b/buch/chapters/110-elliptisch/experiments/Makefile new file mode 100644 index 0000000..fac4fbc --- /dev/null +++ b/buch/chapters/110-elliptisch/experiments/Makefile @@ -0,0 +1,15 @@ +# +# Makefile +# +# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschlue +# +all: KK.pdf + +KN: KN.cpp + g++ -O -Wall -std=c++11 KN.cpp -o KN `pkg-config --cflags gsl` `pkg-config --libs gsl` + +KKpath.tex: KN + ./KN --outfile KKpath.tex + +KK.pdf: KK.tex KKpath.tex + pdflatex KK.tex diff --git a/buch/chapters/110-elliptisch/experiments/agm.maxima b/buch/chapters/110-elliptisch/experiments/agm.maxima deleted file mode 100644 index c7facd4..0000000 --- a/buch/chapters/110-elliptisch/experiments/agm.maxima +++ /dev/null @@ -1,26 +0,0 @@ -/* - * 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/uebungsaufgaben/2.tex b/buch/chapters/110-elliptisch/uebungsaufgaben/2.tex new file mode 100644 index 0000000..9a1cafc --- /dev/null +++ b/buch/chapters/110-elliptisch/uebungsaufgaben/2.tex @@ -0,0 +1,61 @@ +\label{buch:elliptisch:aufgabe:2}% +Die Landen-Transformation basiert auf der Iteration +\begin{equation} +\begin{aligned} +k_{n+1} +&= +\frac{1-k_n'}{1+k_n'} +& +&\text{und}& +k_{n+1}' +&= +\sqrt{1-k_{n+1}^2} +\end{aligned} +\label{buch:elliptisch:aufgabe:2:iteration} +\end{equation} +mit den Startwerten $k_0 = k$ und $k_0' = \sqrt{1-k_0^2}$. +Zeigen Sie, dass $k_n\to 0$ und $k_n'\to 1$ mit quadratischer Konvergenz. + +\begin{loesung} +\begin{table} +\centering +\begin{tabular}{|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|} +\hline +n & k & k' \\ +\hline +0 & 0.200000000000000 & 0.979795897113271 \\ +1 & 0.010205144336438 & 0.999947926158694 \\ +2 & 0.000026037598592 & 0.999999999661022 \\ +3 & 0.000000000169489 & 1.000000000000000 \\ +4 & 0.000000000000000 & 1.000000000000000 \\ +\hline +\end{tabular} +\caption{Numerisches Experiment zur Folge $(k_n,k_n')$ +gemäss \eqref{buch:elliptisch:aufgabe:2:iteration} +mit $k_0=0.2$ +\label{buch:ellptisch:aufgabe:2:numerisch}} +\end{table} +Es ist klar, dass $k'_n\to 1$ folgt, wenn man zeigen kann, dass +$k_n\to 0$ gilt. +Wir berechnen daher +\begin{align*} +k_{n+1} +&= +\frac{1-k_n'}{1+k_n'} += +\frac{1-\sqrt{1-k_n^2}}{1+\sqrt{1-k_n^2}} +\intertext{und erweitern mit dem Nenner $1+\sqrt{1-k_n^2}$ um} +&= +\frac{1-(1-k_n^2)}{(1+\sqrt{1-k_n^2})^2} += +\frac{ k_n^2 }{(1+\sqrt{1-k_n^2})^2} +\le +k_n^2 +\end{align*} +zu erhalten. +Daraus folgt jetzt sofort die quadratische Konvergenz von $k_n$ gegen $0$. + +Ein einfaches numerisches Experiment (siehe +Tabelle~\ref{buch:ellptisch:aufgabe:2:numerisch}) +bestätigt die quadratische Konvergenz der Folgen. +\end{loesung} diff --git a/buch/chapters/110-elliptisch/uebungsaufgaben/3.tex b/buch/chapters/110-elliptisch/uebungsaufgaben/3.tex new file mode 100644 index 0000000..a5d118f --- /dev/null +++ b/buch/chapters/110-elliptisch/uebungsaufgaben/3.tex @@ -0,0 +1,135 @@ +\label{buch:elliptisch:aufgabe:3}% +Aus der in Aufgabe~\ref{buch:elliptisch:aufgabe:2} konstruierten Folge +$k_n$ kann zu einem vorgegebenen $u$ ausserdem die Folge $u_n$ +mit der Rekursionsformel +\[ +u_{n+1} = \frac{u_n}{1+k_{n+1}} +\] +und Anfangswert $u_0=u$ konstruiert werden. +Die Landen-Transformation (siehe \cite[80]{buch:ellfun-applications}) +\index{Landen-Transformation}% +führt auf die folgenden Formeln für die Jacobischen elliptischen Funktionen: +\begin{equation} +\left.\qquad +\begin{aligned} +\operatorname{sn}(u_n,k_n) +&= +\frac{ +(1+k_{n+1})\operatorname{sn}(u_{n+1},k_{n+1}) +}{ +1 + k_{n+1} \operatorname{sn}(u_{n+1},k_{n+1})^2 +} +\\ +\operatorname{cn}(u_n,k_n) +&= +\frac{ +\operatorname{cn}(u_{n+1},k_{n+1}) +\operatorname{dn}(u_{n+1},k_{n+1}) +}{ +1 + k_{n+1} \operatorname{sn}(u_{n+1},k_{n+1})^2 +} +\\ +\operatorname{dn}(u_n,k_n) +&= +\frac{ +1 - k_{n+1} \operatorname{sn}(u_{n+1},k_{n+1})^2 +}{ +1 + k_{n+1} \operatorname{sn}(u_{n+1},k_{n+1})^2 +} +\end{aligned} +\qquad\right\} +\label{buch:elliptisch:aufgabe:3:gauss} +\end{equation} +Die Transformationsformeln +\eqref{buch:elliptisch:aufgabe:3:gauss} +sind auch als Gauss-Transformation bekannt. +\index{Gauss-Transformation}% +Konstruieren Sie daraus einen numerischen Algorithmus, mit dem sich +gleichzeitig die Werte aller drei Jacobischen elliptischen Funktionen +für vorgegebene Parameterwerte $u$ und $k$ berechnen lassen. + +\begin{loesung} +In der ersten Phase des Algorithmus werden die Folgen $k_n$ und $k_n'$ +sowie $u_n$ bis zum Folgenindex $N$ berechnet, bis $k_N\approx 0$ +angenommen werden darf. +Dann gilt +\begin{align*} +\operatorname{sn}(u_N, k_N) &= \operatorname{sn}(u_N,0) = \sin u_N +\\ +\operatorname{cn}(u_N, k_N) &= \operatorname{cn}(u_N,0) = \cos u_N +\\ +\operatorname{dn}(u_N, k_N) &= \operatorname{dn}(u_N,0) = 1. +\end{align*} +In der zweiten Phase des Algorithmus können für absteigende +$n$ jeweils die Formeln~\eqref{buch:elliptisch:aufgabe:3:gauss} +angewendet werden um nacheinander die Werte der Jacobischen +elliptischen Funktionen für Argument $u_n$ und Parameter $k_n$ +für $n=N-1,N-2,\dots,0$ zu bekommen. +\end{loesung} +\begin{table} +\centering +\begin{tikzpicture}[>=latex,thick] +\def\pfeil#1#2{ + \fill[color=#1!30] (-0.5,1) -- (-0.5,-1) -- (-0.8,-1) + -- (0,-1.5) -- (0.8,-1) -- (0.5,-1) -- (0.5,1) -- cycle; + \node[color=white] at (0,-0.2) [scale=5] {\sf #2\strut}; +} +\begin{scope}[xshift=-4.9cm,yshift=0.2cm] +\pfeil{red}{1} +\end{scope} + +\begin{scope}[xshift=-2.3cm,yshift=0.2cm] +\pfeil{red}{1} +\end{scope} + +\begin{scope}[xshift=0.35cm,yshift=-0.3cm,yscale=-1] +\pfeil{blue}{2} +\end{scope} + +\begin{scope}[xshift=2.92cm,yshift=-0.3cm,yscale=-1] +\pfeil{blue}{2} +\end{scope} + +\begin{scope}[xshift=5.60cm,yshift=-0.3cm,yscale=-1] +\pfeil{blue}{2} +\end{scope} + +\node at (0,0) { +\begin{tabular}{|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|>{$}c<{$}|} +\hline +n & k_n & u_n & \operatorname{sn}(u_n,k_n) & \operatorname{cn}(u_n,k_n) & \operatorname{dn}(u_n,k_n)% +\mathstrut\text{\vrule height12pt depth6pt width0pt} \\ +\hline +\mathstrut\text{\vrule height12pt depth0pt width0pt}% +%\small +0 & 0.90000000000 & 0.60000000000 & 0.54228232286 & 0.84019633556 & 0.87281338478 \\ +1 & 0.39286445838 & 0.43076696830 & 0.41576897816 & 0.90947026163 & 0.98656969610 \\ +2 & 0.04188568608 & 0.41344935827 & 0.40175214109 & 0.91574844642 & 0.99985840483 \\ +3 & 0.00043898784 & 0.41326793867 & 0.40160428679 & 0.91581329801 & 0.99999998445 \\ +4 & 0.00000004817 & 0.41326791876 & 0.40160427056 & 0.91581330513 & 1.00000000000 \\ +5 & 0.00000000000 & 0.41326791876 & 0.40160427056 & 0.91581330513 & 1.00000000000 \\ +%N & 0.00000000000 & 0.41326791876 & 0.40160427056 & 0.91581330513 & 1.00000000000% +N & & 0.41326791876 & \sin u_N & \cos u_N & 1% +%0 & 0.900000000000000 & 0.600000000000000 & 0.542282322869158 & 0.840196335569032 & 0.872813384788490 \\ +%1 & 0.392864458385019 & 0.430766968306220 & 0.415768978168966 & 0.909470261631645 & 0.986569696107075 \\ +%2 & 0.041885686080039 & 0.413449358275499 & 0.401752141098324 & 0.915748446421239 & 0.999858404836479 \\ +%3 & 0.000438987841605 & 0.413267938675096 & 0.401604286793186 & 0.915813298019491 & 0.999999984459261 \\ +%4 & 0.000000048177586 & 0.413267918764845 & 0.401604270565476 & 0.915813305135699 & 1.000000000000000 \\ +%5 & 0.000000000000001 & 0.413267918764845 & 0.401604270565476 & 0.915813305135699 & 1.000000000000000 \\ +%N & 0.000000000000000 & 0.413267918764845 & 0.401604270565476 & 0.915813305135699 & 1.000000000000000 \\ +\mathstrut\text{\vrule height12pt depth6pt width0pt} \\ +\hline +\end{tabular} +}; +\end{tikzpicture} +\caption{Durchführung des auf der Landen-Transformation basierenden +Algorithmus zur Berechnung der Jacobischen elliptischen Funktionen +für $u=0.6$ und $k=0.9$. +Die erste Phase (rot) berechnet die Folgen $k_n$ und $u_n$, die zweite +(blau) +transformiert die Wert der trigonometrischen Funktionen in die Werte +der Jacobischen elliptischen Funktionen. +\label{buch:elliptisch:aufgabe:3:resultate}} +\end{table} + + diff --git a/buch/chapters/110-elliptisch/uebungsaufgaben/landen.m b/buch/chapters/110-elliptisch/uebungsaufgaben/landen.m new file mode 100644 index 0000000..bba5549 --- /dev/null +++ b/buch/chapters/110-elliptisch/uebungsaufgaben/landen.m @@ -0,0 +1,60 @@ +# +# landen.m +# +# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +# +N = 10; + +function retval = M(a,b) + for i = (1:10) + A = (a+b)/2; + b = sqrt(a*b); + a = A; + endfor + retval = a; +endfunction; + +function retval = EllipticKk(k) + retval = pi / (2 * M(1, sqrt(1-k^2))); +endfunction + +k = 0.5; +kprime = sqrt(1-k^2); + +EK = EllipticKk(k); +EKprime = EllipticKk(kprime); + +u = EK + EKprime * i; + +K = zeros(N,3); +K(1,1) = k; +K(1,2) = kprime; +K(1,3) = u; + +format long + +for n = (2:N) + K(n,1) = (1-K(n-1,2)) / (1+K(n-1,2)); + K(n,2) = sqrt(1-K(n,1)^2); + K(n,3) = K(n-1,3) / (1 + K(n,1)); +end + +K(:,[1,3]) + +pi / 2 + +scd = zeros(N,3); +scd(N,1) = sin(K(N,3)); +scd(N,2) = cos(K(N,3)); +scd(N,3) = 1; + +for n = (N:-1:2) + nenner = 1 + K(n,1) * scd(n, 1)^2; + scd(n-1,1) = (1+K(n,1)) * scd(n, 1) / nenner; + scd(n-1,2) = scd(n, 2) * scd(n, 3) / nenner; + scd(n-1,3) = (1 - K(n,1) * scd(n,1)^2) / nenner; +end + +scd(:,1) + +cosh(2.009459377005286) diff --git a/buch/chapters/references.bib b/buch/chapters/references.bib index 571831a..fbbbf30 100644 --- a/buch/chapters/references.bib +++ b/buch/chapters/references.bib @@ -137,3 +137,12 @@ year = 2004 } +@book{buch:ellfun-applications, + author = { Derek F. Lawden }, + title = { Elliptic Functions and Applications }, + series = { Applied Mathematical Sciences }, + volume = { 80 }, + publisher = { Springer-Verlag }, + year = 2010, + ISBN = { 978-1-4419-3090-3 } +} -- cgit v1.2.1