add numerical computation of lambert W function

4 files changed, 153 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}