aboutsummaryrefslogtreecommitdiffstats
path: root/buch/chapters/110-elliptisch/images
diff options
context:
space:
mode:
authorAndreas Müller <andreas.mueller@ost.ch>2021-10-16 21:44:20 +0200
committerAndreas Müller <andreas.mueller@ost.ch>2021-10-16 21:44:20 +0200
commiteeb4c0186284aca82d85c0e45e7471a3be1b9cd7 (patch)
tree7fa4ceecc2724f886a237652c579950b149bb914 /buch/chapters/110-elliptisch/images
parentnew files (diff)
downloadSeminarSpezielleFunktionen-eeb4c0186284aca82d85c0e45e7471a3be1b9cd7.tar.gz
SeminarSpezielleFunktionen-eeb4c0186284aca82d85c0e45e7471a3be1b9cd7.zip
wertebereich elliptischer integrale
Diffstat (limited to '')
-rw-r--r--buch/chapters/110-elliptisch/images/Makefile14
-rw-r--r--buch/chapters/110-elliptisch/images/rechteck.cpp369
-rw-r--r--buch/chapters/110-elliptisch/images/rechteck.pdfbin0 -> 70950 bytes
-rw-r--r--buch/chapters/110-elliptisch/images/rechteck.tex80
4 files changed, 431 insertions, 32 deletions
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
new file mode 100644
index 0000000..5fa07cb
--- /dev/null
+++ b/buch/chapters/110-elliptisch/images/rechteck.pdf
Binary files differ
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}
+