From b4ded112e21249774920681e00456b2a5c2b39de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Mon, 7 Jun 2021 09:03:26 +0200 Subject: add numerical computation of lambert W function --- buch/chapters/020-exponential/code/Makefile | 8 ++ buch/chapters/020-exponential/code/lambertw.c | 102 +++++++++++++++++++++ buch/chapters/020-exponential/lambertw.tex | 42 +++++++++ .../chapters/020-exponential/uebungsaufgaben/1.tex | 1 + buch/chapters/references.bib | 11 +++ 5 files changed, 164 insertions(+) create mode 100644 buch/chapters/020-exponential/code/lambertw.c 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 +#include +#include +#include + +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 +} + + + + -- cgit v1.2.1