aboutsummaryrefslogtreecommitdiffstats
path: root/buch
diff options
context:
space:
mode:
authorAndreas Müller <andreas.mueller@ost.ch>2021-06-07 09:03:26 +0200
committerAndreas Müller <andreas.mueller@ost.ch>2021-06-07 09:03:26 +0200
commitb4ded112e21249774920681e00456b2a5c2b39de (patch)
tree454c9541f2301f28a93141b067a71e4214f6b811 /buch
parentadd new problem (diff)
downloadSeminarSpezielleFunktionen-b4ded112e21249774920681e00456b2a5c2b39de.tar.gz
SeminarSpezielleFunktionen-b4ded112e21249774920681e00456b2a5c2b39de.zip
add numerical computation of lambert W function
Diffstat (limited to 'buch')
-rw-r--r--buch/chapters/020-exponential/code/Makefile8
-rw-r--r--buch/chapters/020-exponential/code/lambertw.c102
-rw-r--r--buch/chapters/020-exponential/lambertw.tex42
-rw-r--r--buch/chapters/020-exponential/uebungsaufgaben/1.tex1
-rw-r--r--buch/chapters/references.bib11
5 files changed, 164 insertions, 0 deletions
diff --git a/buch/chapters/020-exponential/code/Makefile b/buch/chapters/020-exponential/code/Makefile
index b189dae..f7c1e9c 100644
--- a/buch/chapters/020-exponential/code/Makefile
+++ b/buch/chapters/020-exponential/code/Makefile
@@ -4,6 +4,14 @@
#
# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
#
+all: lambertw-test
xxl: xxl.c
gcc -o xxl -W -O2 -lgsl xxl.c
+
+lambertw: lambertw.c
+ gcc -o lambertw -W -O2 lambertw.c
+
+lambertw-test: lambertw
+ ./lambertw -l l.txt -- -0.35 -0.367 -0.3678 -0.36787 -0.367879 -0.2 1 2 3 10 100 1000 -0.1 -0.01 -0.001 -0.0001
+
diff --git a/buch/chapters/020-exponential/code/lambertw.c b/buch/chapters/020-exponential/code/lambertw.c
new file mode 100644
index 0000000..14ee190
--- /dev/null
+++ b/buch/chapters/020-exponential/code/lambertw.c
@@ -0,0 +1,102 @@
+/*
+ * lambertw.c
+ *
+ * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <getopt.h>
+#include <math.h>
+
+static double precision = 1e-16;
+static FILE *logfile = NULL;
+
+double f(double x) {
+ return x * exp(x);
+}
+
+double lambertW0(double y) {
+ double x = log(1 + y);
+ double delta = 1;
+ int count = 0;
+ while ((count < 20) && (delta > precision)) {
+ double error = f(x) - y;
+ if (logfile) {
+ fprintf(logfile, "%20.16f %20.16f\n", x, error);
+ }
+ double xnew = x - error / ((1 + x) * exp(x));
+ delta = fabs(xnew - x);
+ x = xnew;
+ count++;
+ }
+ if (count == 0) {
+ fprintf(stderr, "count exceeded\n");
+ return -1;
+ }
+ if (logfile) {
+ fprintf(logfile, "%d\n", count);
+ }
+ return x;
+}
+
+double lambertW1(double y) {
+ double x = log(-y);
+ double delta = 1;
+ int count = 0;
+ while ((count < 20) && (delta > precision)) {
+ double error = f(x) - y;
+ if (logfile) {
+ fprintf(logfile, "%20.16f %20.16f\n", x, error);
+ }
+ double xnew = x - error / ((1 + x) * exp(x));
+ delta = fabs(xnew - x);
+ x = xnew;
+ count++;
+ }
+ if (count == 0) {
+ fprintf(stderr, "count exceeded\n");
+ return -1;
+ }
+ if (logfile) {
+ fprintf(logfile, "%d\n", count);
+ }
+ return x;
+}
+
+struct option options[] = {
+{ "precision", required_argument, NULL, 'p' },
+{ "logfile", required_argument, NULL, 'l' },
+{ "minus", no_argument, NULL, 'm' },
+{ NULL, 0, NULL, 0 }
+};
+
+int main(int argc, char *argv[]) {
+ int c;
+ int l;
+ int minusone = 0;
+ while (EOF != (c = getopt_long(argc, argv, "p:l:m", options, &l)))
+ switch (c) {
+ case 'l':
+ logfile = fopen(optarg, "w");
+ break;
+ case 'p':
+ precision = atof(optarg);
+ break;
+ case 'm':
+ minusone = 1;
+ break;
+ }
+
+ for (; optind < argc; optind++) {
+ double y = atof(argv[optind]);
+ double x = (y >= 0) ? lambertW0(y) :
+ ((minusone) ? lambertW1(y) : lambertW0(y));
+ printf("%20.16f -> %20.16f (%20.16f)\n", x, y, y - f(x));
+ }
+
+ if (logfile) {
+ fclose(logfile);
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/buch/chapters/020-exponential/lambertw.tex b/buch/chapters/020-exponential/lambertw.tex
index a7a882c..410d145 100644
--- a/buch/chapters/020-exponential/lambertw.tex
+++ b/buch/chapters/020-exponential/lambertw.tex
@@ -319,6 +319,48 @@ W(-cbe^{ac})
Die Gleichung hat eine Lösung wenn $-cbe^{ac} > -1/e$ ist.
\end{proof}
+\subsection{Numerische Berechnung
+\label{buch:subsection:lambertberechnung}}
+Die $W$-Funktionen sind nur dann nützlich, wenn man sie effizient
+berechnen kann.
+Leider ist sie nicht Teil der C- oder C++-Standardbibliothek,
+man muss sich also mit einer spezialisierten Bibliothek oder einer
+eigenen Implementation behelfen.
+
+\subsubsection{Berechnung mit dem Newton-Algorithmus}
+Für $x>-1$ ist die Funktion $W(x)$ ist die Umkehrfunktion der
+streng monoton wachsenden und konvexen Funktion $f(x)=xe^x$.
+In dieser Situation konvergiert der Newton-Algorithmus zur Bestimmung
+der Nullstelle $x=W_0(y)$ von $f(x)-y$ für alle Werte von $y>-1/e$.
+Für $W_{-1}(y)$ ist die Situation etwas komplizierter, da für
+$x<-1$ die Funktion $f(x)$ nicht konvex ist.
+
+Ausgehend vom Startwert $x_0$ ist die Iterationsfolge definiert
+durch
+\[
+x_{n+1}
+=
+x_n - \frac{f(x_n) - y}{f'(x_n)}
+=
+x_n - \frac{x_ne^{x_n}-y}{(1+x_n)e^{x_n}}.
+\]
+Die Theorie verspricht, dass die Folge quadratisch konvergiert, wenn
+der Startwert $x_0$ genügend genau ist.
+Für $W_0(y)$ scheint $x_0=\log(1+y)$ ein guter Startwert zu sein, für
+$W_{-1}(y)$ funktioniert $x_0=\log(-y)$.
+
+Die Steigung des Graphen der Funktion $f(x)$ ist für grosse positive
+Werte von $x$ sehr gross und für grosse negative Werte von $x$ sehr
+klein, was die Konvergenz stark beeinträchtigen kann.
+An der Stelle $x=-1$ mit dem Wert $f(-1)=-1/e$ ist die Steigung $0$
+und die Konvergenz des Newton-Algorithmus ist nur noch linear.
+Für mittelgrosse Werte von $y$ weg von $-1/e$ kann $W_0(y)$ oder $W_{-1}(y)$
+mit einer Genauigkeit $10^{-15}$ innert weniger als 10 Iterationen
+bestimmt werden.
+
+\subsubsection{GNU scientific library}
+Die Lambert $W$-Funktionen $W_0(x)$ und $W_{-1}(x)$ sind auch in der
+GNU scientific library \cite{buch:library:gsl} implementiert.
%
% Verfolgungskurven
diff --git a/buch/chapters/020-exponential/uebungsaufgaben/1.tex b/buch/chapters/020-exponential/uebungsaufgaben/1.tex
index c88bdde..3b71806 100644
--- a/buch/chapters/020-exponential/uebungsaufgaben/1.tex
+++ b/buch/chapters/020-exponential/uebungsaufgaben/1.tex
@@ -29,5 +29,6 @@ Für $W(\log 27)$ findet man
W(\log 27) = 1.098612
\qquad\Rightarrow\qquad
x=3.
+\qedhere
\]
\end{loesung}
diff --git a/buch/chapters/references.bib b/buch/chapters/references.bib
index a5d0201..38e56c8 100644
--- a/buch/chapters/references.bib
+++ b/buch/chapters/references.bib
@@ -139,3 +139,14 @@ abstract = "In this paper, we present Google, a prototype of a large-scale searc
language = {german},
}
+@online{buch:library:gsl,
+ title = {GNU scientific library},
+ version = {2.7},
+ year = 2021,
+ month = 6,
+ day = 1
+}
+
+
+
+