From eeb4c0186284aca82d85c0e45e7471a3be1b9cd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Sat, 16 Oct 2021 21:44:20 +0200 Subject: wertebereich elliptischer integrale --- buch/chapters/110-elliptisch/images/rechteck.cpp | 369 +++++++++++++++++++++-- 1 file changed, 341 insertions(+), 28 deletions(-) (limited to 'buch/chapters/110-elliptisch/images/rechteck.cpp') 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 #include #include +#include +#include +#include +#include double ast = 1; -std::complex integrand(double k, const std::complex& z) { - std::complex s = z * z; - std::complex eins(1); - std::complex 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 operator()(const std::complex& z) const = 0; + static double kprime(double k); +}; + +double integrand::kprime(double k) { + return sqrt(1 - k*k); } -std::complex segment(double k, const std::complex& z1, - const std::complex& z2, int n = 100) { +/** + * \brief Elliptic integral of the first kind + */ +class integrand1 : public integrand { +public: + integrand1(double k) : integrand(k) { } + virtual std::complex operator()( + const std::complex& z) const { + std::complex square = z * z; + std::complex eins(1.0); + std::complex 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 > curve_t; + curvetracer(const integrand& f) : _f(f) { } + + std::complex startpoint(const std::complex& z, + int n = 100) const; + + std::complex segment(const std::complex& z1, + const std::complex& z2, int n) const; + + std::pair,std::complex > segmentl( + const std::complex& start, + const std::complex& dir, + double stepsize, int n = 10) const; + + curve_t trace(const std::complex& startz, + const std::complex& dir, + const std::complex& 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 curvetracer::segment(const std::complex& z1, + const std::complex& z2, int n = 100) const { std::complex dz = z2 - z1; std::complex 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 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 z(0, y); - std::complex w = segment(k, std::complex(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 > curvetracer::segmentl( + const std::complex& start, + const std::complex& dir, + double stepsize, + int n) const { + std::complex s(0.); + std::complex z = start; int counter = 100; - while (counter-- > 0) { - std::complex znext = z + std::complex(xstep); - std::complex 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& startz, + const std::complex& dir, + const std::complex& startw, + int maxsteps) const { + curve_t result; + std::complex z = startz; + std::complex 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 curvetracer::startpoint(const std::complex& z, + int n) const { + std::cerr << "start at " << z.real() << "," << z.imag() << std::endl; + return segment(std::complex(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 dir(0.01, 0); + for (double im = 0.2; im < 3; im += 0.2) { + std::complex startz(0, im); + std::complex 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(0, 0.01); + for (double re = 0.2; re < 1; re += 0.2) { + std::complex startz(re, 0.001); + std::complex startw = ct.startpoint(startz); + startw = std::complex(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 startz(1.); + std::complex 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 startz(1/k); + std::complex 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; } -- cgit v1.2.1