diff options
Diffstat (limited to '')
-rw-r--r-- | buch/chapters/060-integral/gq/Makefile | 22 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/gq.cpp | 64 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/gq.pdf | bin | 0 -> 25439 bytes | |||
-rw-r--r-- | buch/chapters/060-integral/gq/gq.tex | 57 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/kreis.cpp | 66 |
5 files changed, 209 insertions, 0 deletions
diff --git a/buch/chapters/060-integral/gq/Makefile b/buch/chapters/060-integral/gq/Makefile new file mode 100644 index 0000000..2ce859c --- /dev/null +++ b/buch/chapters/060-integral/gq/Makefile @@ -0,0 +1,22 @@ +# +# Makefile -- Gauss quadrature experiments +# +# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +# +all: kreis gq.pdf + +gq.o: gq.cpp + g++ `pkg-config --cflags gsl` -o gq.o -g -Wall -O2 -c gq.cpp + +gq: gq.o + g++ -g -Wall -O2 gq.o `pkg-config --libs gsl` -o gq + +kreis: kreis.cpp + g++ `pkg-config --cflags gsl` -g -Wall -O2 kreis.cpp \ + `pkg-config --libs gsl` -o kreis + +gqpaths.tex: gq + ./gq > gqpaths.tex + +gq.pdf: gq.tex gqpaths.tex + pdflatex gq.tex diff --git a/buch/chapters/060-integral/gq/gq.cpp b/buch/chapters/060-integral/gq/gq.cpp new file mode 100644 index 0000000..3eb3ac8 --- /dev/null +++ b/buch/chapters/060-integral/gq/gq.cpp @@ -0,0 +1,64 @@ +/* + * gq.cpp -- Gauss Quadrature Example + * + * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include <iostream> +#include <iomanip> +#include <cstdlib> +#include <cstdio> +#include <cmath> +#include <gsl/gsl_integration.h> + +int counter; + +double f(double x, void *) { + counter++; + return sqrt(1-x*x); +} + +void do_integration(size_t n, const char *name) { + double scaling = 1e6; + // set up integration gauss quadrature + counter = 0; + gsl_function func; + func.function = f; + func.params = NULL; + + // prepare plotting of results + int steps = 200; + double h = 1. / steps; + + printf("\\def\\%s{\n\t(0,0)", name); + for (int i = 1; i <= steps; i++) { + double a = i * h; + a = 1 - (1 - a) * (1 - a); + double expected = (asin(a) + a * sqrt(1 - a*a)) / 2; + gsl_integration_fixed_workspace *ws + = gsl_integration_fixed_alloc( + gsl_integration_fixed_legendre, + n, 0, a, 0, 0); + double result = 0; + gsl_integration_fixed(&func, &result, ws); + double y = fabs(scaling * (result - expected)); + if (y > 10) { + y = 10; + } + printf("\n\t--({\\dx*%.4f},{\\dy*%.6f})", a, y); + gsl_integration_fixed_free(ws); + } + printf("\n}\n"); +} + +int main(int argc, char *argv[]) { + do_integration(2, "two"); + do_integration(4, "four"); + do_integration(6, "six"); + do_integration(8, "eight"); + do_integration(10, "ten"); + do_integration(12, "twelve"); + do_integration(14, "fourteen"); + do_integration(16, "sixteen"); + + return EXIT_SUCCESS; +} diff --git a/buch/chapters/060-integral/gq/gq.pdf b/buch/chapters/060-integral/gq/gq.pdf Binary files differnew file mode 100644 index 0000000..7a0b3c4 --- /dev/null +++ b/buch/chapters/060-integral/gq/gq.pdf diff --git a/buch/chapters/060-integral/gq/gq.tex b/buch/chapters/060-integral/gq/gq.tex new file mode 100644 index 0000000..e5a9c4d --- /dev/null +++ b/buch/chapters/060-integral/gq/gq.tex @@ -0,0 +1,57 @@ +% +% gq.tex -- Fehler +% +% (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} +\def\dx{12} +\def\dy{8} +\input{gqpaths.tex} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\begin{scope} +\clip (0,0) rectangle ({\dx*1},{\dy*1}); +\draw[color=red,line width=1.4pt] \two; +\draw[color=red,line width=1.4pt] \four; +\draw[color=red,line width=1.4pt] \six; +\draw[color=red,line width=1.4pt] \eight; +\draw[color=red,line width=1.4pt] \ten; +\draw[color=red,line width=1.4pt] \twelve; +\draw[color=red,line width=1.4pt] \fourteen; +\draw[color=red,line width=1.4pt] \sixteen; +\end{scope} + +% add image content here +\draw[->] (0,-0.1) -- (0,{\dy+0.5}) coordinate[label={right:$e$}]; +\draw[->] (-0.1,0) -- ({\dx+0.5},0) coordinate[label={$a$}]; + +\draw ({\dx},-0.1) -- ({\dx},0.1); +\draw[line width=0.3pt] ({\dx},0) -- ({\dx},{\dy}); +\node at ({\dx},0) [below] {$1$}; +\node at (0,0) [below] {$0$}; + +\draw (-0.1,\dy) -- (0.1,\dy); +\node at (0,\dy) [left] {$10^{-6}$}; + +\node at ({0.227*\dx},{0.7*\dy}) [rotate=83.5] {$n=2$}; +\node at ({0.678*\dx},{0.7*\dy}) [rotate=83] {$n=4$}; +\node at ({0.850*\dx},{0.7*\dy}) [rotate=84] {$n=6$}; +\node at ({0.915*\dx},{0.7*\dy}) [rotate=85] {$n=8$}; +\node at ({1.015*\dx},{0.7*\dy}) [rotate=90] {$n=16$}; + +\foreach \x in {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}{ + \draw ({\x*\dx},-0.1) -- ({\x*\dx},0.1); + \node at ({\x*\dx},0) [below] {$\x$}; +} + +\end{tikzpicture} +\end{document} + diff --git a/buch/chapters/060-integral/gq/kreis.cpp b/buch/chapters/060-integral/gq/kreis.cpp new file mode 100644 index 0000000..c9723cc --- /dev/null +++ b/buch/chapters/060-integral/gq/kreis.cpp @@ -0,0 +1,66 @@ +/* + * kreis.cpp -- Gauss Quadrature Example + * + * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include <iostream> +#include <iomanip> +#include <cstdlib> +#include <cstdio> +#include <cmath> +#include <gsl/gsl_integration.h> + +double xmax = 0.990; + +double f(double x, void *) { + return sqrt(1-x*x); +} + +double do_integration(size_t n) { + // set up integration gauss quadrature + gsl_function func; + func.function = f; + func.params = NULL; + + gsl_integration_fixed_workspace *ws + = gsl_integration_fixed_alloc( + gsl_integration_fixed_legendre, + n, -xmax, xmax, 0, 0); + double result = 0; + gsl_integration_fixed(&func, &result, ws); + gsl_integration_fixed_free(ws); + return result; +} + +double do_trapez(size_t n) { + double h = 2 * xmax / n; + double s = 0; + for (int i = 0; i <= n; i++) { + double x = -xmax + i * h; + double a = f(x, NULL); + if ((i == 0) || (i == n)) { + a = a * 0.5; + } + s += a; + } + return h * s; +} + +void do_table(double limit) { + xmax = limit; + for (int n = 2; n <= 20; n += 2) { + printf("%d & %.16f & %.16f \\\\\n", n, do_integration(n), + do_trapez(10 * n)); + } + double amax = acos(xmax); + double expected = (M_PI / 2 - amax) + sin(amax) * cos(amax); + printf("\\infty & %.16f & %.16f \\\\\n", expected, expected); +} + +int main(int argc, char *argv[]) { + do_table(0.5); + do_table(0.9); + do_table(0.99); + do_table(0.999); + return EXIT_SUCCESS; +} |