diff options
Diffstat (limited to '')
| -rw-r--r-- | buch/chapters/110-elliptisch/ellintegral.tex | 10 | ||||
| -rw-r--r-- | buch/chapters/110-elliptisch/images/Makefile | 14 | ||||
| -rw-r--r-- | buch/chapters/110-elliptisch/images/rechteck.cpp | 369 | ||||
| -rw-r--r-- | buch/chapters/110-elliptisch/images/rechteck.pdf | bin | 0 -> 70950 bytes | |||
| -rw-r--r-- | buch/chapters/110-elliptisch/images/rechteck.tex | 80 |
5 files changed, 441 insertions, 32 deletions
diff --git a/buch/chapters/110-elliptisch/ellintegral.tex b/buch/chapters/110-elliptisch/ellintegral.tex index 7ac09ca..b6f35fc 100644 --- a/buch/chapters/110-elliptisch/ellintegral.tex +++ b/buch/chapters/110-elliptisch/ellintegral.tex @@ -329,6 +329,16 @@ XXX Parameterkonventionen \\ XXX Wertebereich (Rechtecke) \\ XXX Komplementäre Integrale \\ +\begin{figure} +\centering +\includegraphics{chapters/110-elliptisch/images/rechteck.pdf} +\caption{Der Wertebereich der Funktion $F(k,z)$ ist ein Rechteck +der Breite $2K(k)$ und $2K(k')$. +Die obere Halbebene wird in das rote Rechteck abgebildet, die unter +in das blaue. +\label{buch:elliptisch:fig:rechteck}} +\end{figure} + \subsection{Potenzreihe} XXX Potenzreihen \\ XXX Als hypergeometrische Funktionen \\ diff --git a/buch/chapters/110-elliptisch/images/Makefile b/buch/chapters/110-elliptisch/images/Makefile index 4912cca..6e3a5ce 100644 --- a/buch/chapters/110-elliptisch/images/Makefile +++ b/buch/chapters/110-elliptisch/images/Makefile @@ -3,7 +3,7 @@ # # (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule # -all: lemniskate.pdf ellipsenumfang.pdf unvollstaendig.pdf +all: lemniskate.pdf ellipsenumfang.pdf unvollstaendig.pdf rechteck.pdf lemniskate.pdf: lemniskate.tex pdflatex lemniskate.tex @@ -15,13 +15,19 @@ ekplot.tex: ellipsenumfang.m octave ellipsenumfang.m rechteck: rechteck.cpp - g++ -O -Wall -g rechteck.cpp -o rechteck + g++ `pkg-config --cflags gsl` `pkg-config --libs gsl` -O -Wall -g -std=c++11 rechteck.cpp -o rechteck -test: rechteck - ./rechteck +rechteckpfade.tex: rechteck + ./rechteck --outfile rechteckpfade.tex + +rechteck.pdf: rechteck.tex rechteckpfade.tex + pdflatex rechteck.tex unvollstaendig.pdf: unvollstaendig.tex unvollpath.tex pdflatex unvollstaendig.tex unvollpath.tex: unvollstaendig.m octave unvollstaendig.m + + + diff --git a/buch/chapters/110-elliptisch/images/rechteck.cpp b/buch/chapters/110-elliptisch/images/rechteck.cpp index 6bd6a06..3ece11c 100644 --- a/buch/chapters/110-elliptisch/images/rechteck.cpp +++ b/buch/chapters/110-elliptisch/images/rechteck.cpp @@ -8,47 +8,360 @@ #include <cstdio> #include <complex> #include <iostream> +#include <fstream> +#include <list> +#include <getopt.h> +#include <gsl/gsl_sf_ellint.h> double ast = 1; -std::complex<double> integrand(double k, const std::complex<double>& z) { - std::complex<double> s = z * z; - std::complex<double> eins(1); - std::complex<double> result = eins / sqrt((eins - s) * (eins - k * k * s)); - return ast * ((result.imag() < 0) ? -result : result); +/** + * \brief Base class for integrands + */ +class integrand { +protected: + double _k; +public: + double k() const { return _k; } + integrand(double k) : _k(k) { } + virtual std::complex<double> operator()(const std::complex<double>& z) const = 0; + static double kprime(double k); +}; + +double integrand::kprime(double k) { + return sqrt(1 - k*k); } -std::complex<double> segment(double k, const std::complex<double>& z1, - const std::complex<double>& z2, int n = 100) { +/** + * \brief Elliptic integral of the first kind + */ +class integrand1 : public integrand { +public: + integrand1(double k) : integrand(k) { } + virtual std::complex<double> operator()( + const std::complex<double>& z) const { + std::complex<double> square = z * z; + std::complex<double> eins(1.0); + std::complex<double> result = eins + / sqrt((eins - square) * (eins - k() * k() * square)); + return ast * ((result.imag() < 0) ? -result : result); + } +}; + +/** + * \brief A class to trace curves + */ +class curvetracer { + const integrand& _f; +public: + typedef std::list<std::complex<double> > curve_t; + curvetracer(const integrand& f) : _f(f) { } + + std::complex<double> startpoint(const std::complex<double>& z, + int n = 100) const; + + std::complex<double> segment(const std::complex<double>& z1, + const std::complex<double>& z2, int n) const; + + std::pair<std::complex<double>,std::complex<double> > segmentl( + const std::complex<double>& start, + const std::complex<double>& dir, + double stepsize, int n = 10) const; + + curve_t trace(const std::complex<double>& startz, + const std::complex<double>& dir, + const std::complex<double>& startw, + int maxsteps) const; + + static curve_t mirrorx(const curve_t& c); + static curve_t mirrory(const curve_t& c); + static curve_t mirror(const curve_t& c); +}; + +/** + * \brief Perform integration for a + * + * \param z1 the start point of the segment in the domain + * \param z2 the end point of the segment in the domain + * \param n the number of integration steps to use + * \return the increment along that segment + */ +std::complex<double> curvetracer::segment(const std::complex<double>& z1, + const std::complex<double>& z2, int n = 100) const { std::complex<double> dz = z2 - z1; std::complex<double> summe(0); - double h = 1. / n; - summe = integrand(k, z1); - summe += integrand(k, z2); - for (int i = 1; i < n; i++) { + double h = 1. / (2 * n); + for (int i = 1; i < 2 * n; i += 2) { double t = i * h; std::complex<double> z = (1 - t) * z1 + t * z2; - summe += 2. * integrand(k, z); + summe += _f(z); } - return dz * h * summe / 2.; + return dz * h * summe * 2.; } -int main(int argc, char *argv[]) { - double k = 0.5; - double y = -0.0001; - double xstep = -0.1; - ast = 1; - std::complex<double> z(0, y); - std::complex<double> w = segment(k, std::complex<double>(0), z); - std::cout << z << std::endl; +/** + * \brief Exception thrown when the number of iterations is exceeded + */ +class toomanyiterations : public std::runtime_error { +public: + toomanyiterations(const std::string& cause) + : std::runtime_error(cause) { + } +}; + +/** + * \brief Perform integration with a given target length + * + * \param start the start point + * \param dir the direction in which to integrate + * \param stepsize the required length of the step + * \param n the number of integration steps + * \return the increment in the range by this segment + */ +std::pair<std::complex<double>, std::complex<double> > curvetracer::segmentl( + const std::complex<double>& start, + const std::complex<double>& dir, + double stepsize, + int n) const { + std::complex<double> s(0.); + std::complex<double> z = start; int counter = 100; - while (counter-- > 0) { - std::complex<double> znext = z + std::complex<double>(xstep); - std::complex<double> incr = segment(k, z, znext); - w += incr; - std::cout << znext << " -> " << w << ", "; - std::cout << integrand(k, znext) << std::endl; - z = znext; + while (abs(s) < stepsize) { + s = s + segment(z, z + dir, n); + z = z + dir; + if (counter-- == 0) { + throw toomanyiterations("too many iterations"); + } + } + return std::make_pair(z, s); +} + +/** + * \brief Trace a curve from a starting point in the domain + * + * \param startz the start point of the curve in the domain + * \param dir the direction of the curve in the domain + * \param startw the initial function value where the curve starts in + * the range + * \param maxsteps the maximum number of dir-steps before giving up + * \return the curve as a list of points + */ +curvetracer::curve_t curvetracer::trace(const std::complex<double>& startz, + const std::complex<double>& dir, + const std::complex<double>& startw, + int maxsteps) const { + curve_t result; + std::complex<double> z = startz; + std::complex<double> w = startw; + result.push_back(w); + while (maxsteps-- > 0) { + try { + auto seg = segmentl(z, dir, abs(dir), 40); + z = seg.first; + w = w + seg.second; + result.push_back(w); + } catch (const toomanyiterations& x) { + std::cerr << "iterations exceeded after "; + std::cerr << result.size(); + std::cerr << " points"; + maxsteps = 0; + } + } + return result; +} + +/** + * \brief Find the initial point for a coordinate curve + * + * \param k the elliptic integral parameter + * \param z the desired starting point argument + */ +std::complex<double> curvetracer::startpoint(const std::complex<double>& z, + int n) const { + std::cerr << "start at " << z.real() << "," << z.imag() << std::endl; + return segment(std::complex<double>(0.), z, n); +} + + +curvetracer::curve_t curvetracer::mirrorx(const curve_t& c) { + curve_t result; + for (auto z : c) { + result.push_back(-std::conj(z)); + } + return result; +} + +curvetracer::curve_t curvetracer::mirrory(const curve_t& c) { + curve_t result; + for (auto z : c) { + result.push_back(std::conj(z)); + } + return result; +} + +curvetracer::curve_t curvetracer::mirror(const curve_t& c) { + curve_t result; + for (auto z : c) { + result.push_back(-z); + } + return result; +} + +/** + * \brief Class to draw the curves to a file + */ +class curvedrawer { + std::ostream *_out; + std::string _color; +public: + curvedrawer(std::ostream *out) : _out(out), _color("red") { } + const std::string& color() const { return _color; } + void color(const std::string& c) { _color = c; } + void operator()(const curvetracer::curve_t& curve); + std::ostream *out() { return _out; } +}; + +/** + * \brief Operator to draw a curve + * + * \param curve the curve to draw + */ +void curvedrawer::operator()(const curvetracer::curve_t& curve) { + double first = true; + for (auto z : curve) { + if (first) { + *_out << "\\draw[color=" << _color << "] "; + first = false; + } else { + *_out << std::endl << " -- "; + } + *_out << "({" << z.real() << "*\\dx},{" << z.imag() << "*\\dy})"; } + *_out << ";" << std::endl; +} + +static struct option longopts[] = { +{ "outfile", required_argument, NULL, 'o' }, +{ "k", required_argument, NULL, 'k' }, +{ NULL, 0, NULL, 0 } +}; + +/** + * \brief Main function + */ +int main(int argc, char *argv[]) { + double k = 0.6; + + int c; + int longindex; + std::string outfilename; + while (EOF != (c = getopt_long(argc, argv, "o:k:", longopts, &longindex))) + switch (c) { + case 'o': + outfilename = std::string(optarg); + break; + case 'k': + k = std::stod(optarg); + break; + } + + double kprime = integrand::kprime(k); + + double xmax = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); + double ymax = gsl_sf_ellint_Kcomp(kprime, GSL_PREC_DOUBLE); + std::cout << "xmax = " << xmax << std::endl; + std::cout << "ymax = " << ymax << std::endl; + + curvedrawer *cdp; + std::ofstream *outfile = NULL; + if (outfilename.c_str()) { + outfile = new std::ofstream(outfilename.c_str()); + } + if (outfile) { + cdp = new curvedrawer(outfile); + } else { + cdp = new curvedrawer(&std::cout); + } + + integrand1 f(k); + curvetracer ct(f); + + // fill + (*cdp->out()) << "\\fill[color=red!10] ({" << (-xmax) << "*\\dx},0) " + << "rectangle ({" << xmax << "*\\dx},{" << ymax << "*\\dy});" + << std::endl; + (*cdp->out()) << "\\fill[color=blue!10] ({" << (-xmax) << "*\\dx},{" + << (-ymax) << "*\\dy}) rectangle ({" << xmax << "*\\dx},0);" + << std::endl; + + // "circles" + std::complex<double> dir(0.01, 0); + for (double im = 0.2; im < 3; im += 0.2) { + std::complex<double> startz(0, im); + std::complex<double> startw = ct.startpoint(startz); + curvetracer::curve_t curve = ct.trace(startz, dir, + startw, 1000); + cdp->color("red"); + (*cdp)(curve); + (*cdp)(curvetracer::mirrorx(curve)); + cdp->color("blue"); + (*cdp)(curvetracer::mirrory(curve)); + (*cdp)(curvetracer::mirror(curve)); + } + + // imaginary axis + (*cdp->out()) << "\\draw[color=red] (0,0) -- (0,{" << ymax + << "*\\dy});" << std::endl; + (*cdp->out()) << "\\draw[color=blue] (0,0) -- (0,{" << (-ymax) + << "*\\dy});" << std::endl; + + // arguments between 0 and 1 + dir = std::complex<double>(0, 0.01); + for (double re = 0.2; re < 1; re += 0.2) { + std::complex<double> startz(re, 0.001); + std::complex<double> startw = ct.startpoint(startz); + startw = std::complex<double>(startw.real(), 0); + curvetracer::curve_t curve = ct.trace(startz, dir, + startw, 1000); + cdp->color("red"); + (*cdp)(curve); + (*cdp)(curvetracer::mirrorx(curve)); + cdp->color("blue"); + (*cdp)(curvetracer::mirrory(curve)); + (*cdp)(curvetracer::mirror(curve)); + } + + // argument 1 (singularity) + { + std::complex<double> startz(1.); + std::complex<double> startw(xmax); + 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 + { + std::complex<double> startz(1/k); + std::complex<double> startw(xmax, 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; + (*cdp->out()) << "\\def\\ymax{" << ymax << "}" << std::endl; + return EXIT_SUCCESS; } diff --git a/buch/chapters/110-elliptisch/images/rechteck.pdf b/buch/chapters/110-elliptisch/images/rechteck.pdf Binary files differnew file mode 100644 index 0000000..5fa07cb --- /dev/null +++ b/buch/chapters/110-elliptisch/images/rechteck.pdf diff --git a/buch/chapters/110-elliptisch/images/rechteck.tex b/buch/chapters/110-elliptisch/images/rechteck.tex new file mode 100644 index 0000000..622a9e9 --- /dev/null +++ b/buch/chapters/110-elliptisch/images/rechteck.tex @@ -0,0 +1,80 @@ +% +% rechteck.tex -- rechteck für Wertebereich der ell. Integrale +% +% (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{3} +\def\dy{3} + +\input{rechteckpfade.tex} + +\begin{scope} + \clip ({-\xmax*\dx},{-\ymax*\dy}) rectangle ({\xmax*\dx},{\ymax*\dy}); + \fill[color=red!10] (0,{\ymax*\dy}) circle[radius=0.45]; + \fill[color=blue!10] (0,{-\ymax*\dy}) circle[radius=0.45]; + \draw[color=gray!80] (0,{\ymax*\dy}) -- (0,{\ymax*\dy-0.45}); + \draw[color=gray!80] (0,{-\ymax*\dy}) -- (0,{-\ymax*\dy+0.45}); +\end{scope} + +\draw ({-\xmax*\dx},{-\ymax*\dy}) rectangle ({\xmax*\dx},{\ymax*\dy}); + +\draw[->] ({-\xmax*\dx-0.3},0) -- ({\xmax*\dx+0.9},0) + coordinate[label={$\operatorname{Re}F(k,z)$}]; + +\draw[<->,line width=0.7pt] + ({-\xmax*\dx},{\ymax*\dy+0.2}) + -- + ({\xmax*\dx},{\ymax*\dy+0.2}); +\draw[line width=0.2pt] + ({-\xmax*\dx},{\ymax*\dy}) + -- + ({-\xmax*\dx},{\ymax*\dy+0.3}); +\draw[line width=0.2pt] + ({\xmax*\dx},{\ymax*\dy}) + -- + ({\xmax*\dx},{\ymax*\dy+0.3}); + +\node at ({-0.5*\xmax*\dx},{\ymax*\dy+0.2}) [above] {$2K(k)$}; +\draw[->] (0,{\ymax*\dy}) -- (0,{\ymax*\dy+0.7}) + coordinate[label={right:$\operatorname{Im}F(k,z)$}]; +\draw (0,{-\ymax*\dy}) -- (0,{-\ymax*\dy-0.2}); + +\draw[line width=0.2pt] + ({-\xmax*\dx-0.3},{-\ymax*\dy}) + -- + ({-\xmax*\dx},{-\ymax*\dy}); +\draw[line width=0.2pt] + ({-\xmax*\dx-0.3},{\ymax*\dy}) + -- + ({-\xmax*\dx},{\ymax*\dy}); +\draw[<->,line width=0.7pt] + ({-\xmax*\dx-0.2},{\ymax*\dy}) + -- + ({-\xmax*\dx-0.2},0); + +\node at ({-\xmax*\dx-0.2},{-0.5*\ymax*\dy}) [rotate=90,above] {$K(k')$}; +\node at ({-\xmax*\dx-0.2},{0.5*\ymax*\dy}) [rotate=90,above] {$K(k')$}; + +\draw[<->,line width=0.7pt] + ({-\xmax*\dx-0.2},{-\ymax*\dy}) + -- + ({-\xmax*\dx-0.2},0); + +\node at ({\xmax*\dx},{\ymax*\dy}) [right] {$K(k)+iK(k')$}; +\fill[color=white] ({\xmax*\dx},{\ymax*\dy}) circle[radius=0.05]; +\draw ({\xmax*\dx},{\ymax*\dy}) circle[radius=0.05]; + +\end{tikzpicture} +\end{document} + |