From 24e0ebb5475b275af76900be4cbe5d7a8df0ef06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Mon, 18 Oct 2021 08:51:09 +0200 Subject: improve elliptic integrals --- buch/chapters/110-elliptisch/ellintegral.tex | 34 +++++++++++---- buch/chapters/110-elliptisch/images/rechteck.cpp | 53 +++++++++++++++++++++-- buch/chapters/110-elliptisch/images/rechteck.pdf | Bin 70950 -> 91630 bytes 3 files changed, 76 insertions(+), 11 deletions(-) (limited to 'buch/chapters/110-elliptisch') diff --git a/buch/chapters/110-elliptisch/ellintegral.tex b/buch/chapters/110-elliptisch/ellintegral.tex index 1bec096..86b6431 100644 --- a/buch/chapters/110-elliptisch/ellintegral.tex +++ b/buch/chapters/110-elliptisch/ellintegral.tex @@ -414,7 +414,7 @@ Die Ecken auf der reellen Achse liegen bei den reellen Koordinaten \] Für die Höhe muss das Integral \begin{equation} -l=\int_1^{\frac1{k}} +l(\frac{1}{k})=\int_1^{\frac1{k}} \frac{dt}{\sqrt{(t^2-1)(1-k^2t^2)}} \label{buch:elliptisch:eqn:hoeheintegral} \end{equation} @@ -451,11 +451,20 @@ Mit der Ableitung \frac{k'^2 y}{(1-k'^2y^2)^{\frac32}} \] der Substitution -wird das Integral~\eqref{buch:elliptisch:eqn:hoeheintegral} jetzt zu +wird das Integral~\eqref{buch:elliptisch:eqn:hoeheintegral} mit der +oberen Grenze $x$ zu einem Integral mit oberer Grenze +\[ +x^2 = \frac{1}{1-k'^2y_0^2} +\quad\Rightarrow\quad +y_0^2 = \frac{1}{k'^2}\biggl(1-\frac{1}{x^2}\biggr) +\quad\Rightarrow\quad +y_0=\frac{1}{k'}\sqrt{1-\frac{1}{x^2}} +\] +jetzt zu \begin{align*} -l +l(x) &= -\int_0^1 +\int_0^{y_0} \frac{1}{\sqrt{\frac{1}{1-k'^2y^2}-1}} \cdot \frac{1}{\sqrt{1-\frac{k^2}{1-k'^2y^2}}} @@ -466,7 +475,7 @@ l \,dy \\ &= -\int_0^1 +\int_0^{y_0} \frac{\sqrt{1-k'^2y^2}}{\sqrt{k'^2y^2}} \cdot \frac{1}{\sqrt{1-k^2 -k'^2y^2}} @@ -475,7 +484,7 @@ l \,dy \\ &= -\int_0^1 +\int_0^{y_0} \sqrt{1-k'^2y^2} \cdot \frac{1}{k'\sqrt{1-y^2}} @@ -484,10 +493,19 @@ l \,dy \\ &= -\int_0^1 \frac{dy}{\sqrt{(1-y^2)(1-k'^2y^2)}} +\int_0^{y_0} \frac{dy}{\sqrt{(1-y^2)(1-k'^2y^2)}} = -K(k'). +F(y_0,k') \end{align*} +Die gesuchte Höhe des Rechtecks ergibt sich für die obere Grenze $\frac1k$. +In diesem Fall ist +\[ +y_0 += +\frac{1}{k'}\sqrt{1-k^2} = 1 +\] +und das unvollständig elliptische Integral wird zum vollständigen +elliptischen Integral $K(k')$. Die Höhe des Rechtecks des Wertebereichs der oberen Halbebene ist als der Wert des vollständigen elliptischen Integrals erster Art für den Komplementärmodul. diff --git a/buch/chapters/110-elliptisch/images/rechteck.cpp b/buch/chapters/110-elliptisch/images/rechteck.cpp index 3ece11c..c3f6f1b 100644 --- a/buch/chapters/110-elliptisch/images/rechteck.cpp +++ b/buch/chapters/110-elliptisch/images/rechteck.cpp @@ -243,6 +243,7 @@ void curvedrawer::operator()(const curvetracer::curve_t& curve) { static struct option longopts[] = { { "outfile", required_argument, NULL, 'o' }, { "k", required_argument, NULL, 'k' }, +{ "deltax", required_argument, NULL, 'd' }, { NULL, 0, NULL, 0 } }; @@ -250,13 +251,18 @@ static struct option longopts[] = { * \brief Main function */ int main(int argc, char *argv[]) { - double k = 0.6; + double k = 0.625; + double deltax = 0.2; int c; int longindex; std::string outfilename; - while (EOF != (c = getopt_long(argc, argv, "o:k:", longopts, &longindex))) + while (EOF != (c = getopt_long(argc, argv, "o:k:d:", longopts, + &longindex))) switch (c) { + case 'd': + deltax = std::stod(optarg); + break; case 'o': outfilename = std::string(optarg); break; @@ -296,7 +302,7 @@ int main(int argc, char *argv[]) { // "circles" std::complex dir(0.01, 0); - for (double im = 0.2; im < 3; im += 0.2) { + for (double im = deltax; im < 3; im += deltax) { std::complex startz(0, im); std::complex startw = ct.startpoint(startz); curvetracer::curve_t curve = ct.trace(startz, dir, @@ -345,7 +351,28 @@ int main(int argc, char *argv[]) { (*cdp)(curvetracer::mirrory(curve)); } + // arguments between 1 and 1/k + { + for (double x0 = 1 + deltax; x0 < 1/k; x0 += deltax) { + double y0 = sqrt(1-1/(x0*x0))/kprime; + //std::cout << "y0 = " << y0 << std::endl; + double y = gsl_sf_ellint_F(asin(y0), kprime, + GSL_PREC_DOUBLE); + std::complex startz(x0); + std::complex startw(xmax, y); + curvetracer::curve_t curve = ct.trace(startz, dir, + startw, 1000); + cdp->color("red"); + (*cdp)(curve); + (*cdp)(curvetracer::mirrorx(curve)); + cdp->color("blue"); + (*cdp)(curvetracer::mirror(curve)); + (*cdp)(curvetracer::mirrory(curve)); + } + } + // argument 1/k +#if 0 { std::complex startz(1/k); std::complex startw(xmax, ymax); @@ -358,6 +385,26 @@ int main(int argc, char *argv[]) { (*cdp)(curvetracer::mirror(curve)); (*cdp)(curvetracer::mirrory(curve)); } +#endif + + // arguments larger than 1/k + { + double x0 = 1; + while (x0 <= 1/k + 0.0001) { x0 += deltax; } + for (; x0 < 4; x0 += deltax) { + std::complex startz(x0); + std::complex startw(gsl_sf_ellint_F( + asin(1/(k*x0)), k, GSL_PREC_DOUBLE), ymax); + curvetracer::curve_t curve = ct.trace(startz, dir, + startw, 1000); + cdp->color("red"); + (*cdp)(curve); + (*cdp)(curvetracer::mirrorx(curve)); + cdp->color("blue"); + (*cdp)(curvetracer::mirror(curve)); + (*cdp)(curvetracer::mirrory(curve)); + } + } // border (*cdp->out()) << "\\def\\xmax{" << xmax << "}" << std::endl; diff --git a/buch/chapters/110-elliptisch/images/rechteck.pdf b/buch/chapters/110-elliptisch/images/rechteck.pdf index 5fa07cb..845a550 100644 Binary files a/buch/chapters/110-elliptisch/images/rechteck.pdf and b/buch/chapters/110-elliptisch/images/rechteck.pdf differ -- cgit v1.2.1