[genius] Tue Oct 15 18:17:58 2013 Jiri (George) Lebl <jirka 5z com>

From: George Lebl <jirka src gnome org>
To: commits-list gnome org
Cc:
Subject: [genius] Tue Oct 15 18:17:58 2013 Jiri (George) Lebl <jirka 5z com>
Date: Tue, 15 Oct 2013 23:17:56 +0000 (UTC)
commit 7abdb25e15f2b95a013d296c8f742b38f887c39f
Author: Jiri (George) Lebl <jiri lebl gmail com>
Date:   Tue Oct 15 18:18:04 2013 -0500

    Tue Oct 15 18:17:58 2013  Jiri (George) Lebl <jirka 5z com>
    
        * lib/functions/lambert.gel, lib/equation_solving/newton.gel: Add
          LambertW, LambertWm1, NewtonsMethod
    
        * src/geniustests.txt: add tests for the above

 ChangeLog                        |    7 +++++
 lib/equation_solving/Makefile.am |    2 +-
 lib/equation_solving/newton.gel  |   20 +++++++++++++++
 lib/functions/Makefile.am        |    2 +-
 lib/functions/lambert.gel        |   51 ++++++++++++++++++++++++++++++++++++++
 lib/library-strings.c            |    3 ++
 src/geniustests.txt              |    6 ++++
 7 files changed, 89 insertions(+), 2 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 345b3a9..a42ab7b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+Tue Oct 15 18:17:58 2013  Jiri (George) Lebl <jirka 5z com>
+
+       * lib/functions/lambert.gel, lib/equation_solving/newton.gel: Add
+         LambertW, LambertWm1, NewtonsMethod
+       
+       * src/geniustests.txt: add tests for the above
+
 Tue Oct 08 17:28:57 2013  Jiri (George) Lebl <jirka 5z com>
 
        * src/graphing.[ch], src/eval.[ch], src/genius.c, src/gnome-genius.c:
diff --git a/lib/equation_solving/Makefile.am b/lib/equation_solving/Makefile.am
index 6664bf9..90df0d8 100644
--- a/lib/equation_solving/Makefile.am
+++ b/lib/equation_solving/Makefile.am
@@ -1,7 +1,7 @@
 Librarydir = $(datadir)/genius/gel/equation_solving
 SUBDIRS =
 
-GELFILES = find_root.gel diffeqs.gel formulas.gel
+GELFILES = find_root.gel diffeqs.gel formulas.gel newton.gel
 
 EXTRA_DIST = $(GELFILES)
 
diff --git a/lib/equation_solving/newton.gel b/lib/equation_solving/newton.gel
new file mode 100644
index 0000000..3975311
--- /dev/null
+++ b/lib/equation_solving/newton.gel
@@ -0,0 +1,20 @@
+# Newton's method
+
+SetHelp("NewtonsMethod","equation_solving","Run Newton's method on a function to attempt to find a root, 
returns after two successive values are within epsilon or after maxn tries (then returns null)")
+function NewtonsMethod(f,df,guess,epsilon,maxn) = (
+       guess := float(guess);
+       for n=1 to maxn do (
+               dfg := df(guess);
+               if dfg == 0 then (
+                       error ("NewtonsMethod: division by zero");
+                       bailout
+               );
+               guessn := guess - f(guess)/dfg;
+               if |guessn-guess| <= epsilon then
+                       return guessn;
+               guess := guessn
+       );
+       null
+)
+
+
diff --git a/lib/functions/Makefile.am b/lib/functions/Makefile.am
index 1beebff..09a9f2f 100644
--- a/lib/functions/Makefile.am
+++ b/lib/functions/Makefile.am
@@ -1,7 +1,7 @@
 Librarydir = $(datadir)/genius/gel/functions
 SUBDIRS =
 
-GELFILES = delta.gel elementary.gel numerical.gel orthogonal_polynomials.gel kernels.gel complex_numbers.gel
+GELFILES = delta.gel elementary.gel numerical.gel orthogonal_polynomials.gel kernels.gel complex_numbers.gel 
lambert.gel
 
 EXTRA_DIST = $(GELFILES)
 
diff --git a/lib/functions/lambert.gel b/lib/functions/lambert.gel
new file mode 100644
index 0000000..d143e90
--- /dev/null
+++ b/lib/functions/lambert.gel
@@ -0,0 +1,51 @@
+# principal branch of the Lambert W function and its minus one branch
+
+SetHelp ("LambertW", "functions", "Principal branch of the Lambert W function for real values greater than 
or equal to -1/e");
+function LambertW(x) = (
+  if not IsReal(x) or x < -1/e then (error ("LambertW: only defined for real x >= -1/e");bailout);
+  if (x < 0.367) then
+    y := x-x^2+(3/2.0)*x^3
+  else if (x < 24) then 
+    (lnxm1 := ln(x)-1; y := 1+lnxm1/2.0+(lnxm1^2)/16.0)
+  else
+    (l1:=ln(x);l2:=ln(ln(x));y:=l1-l2+l2/l1);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  for n=1 to 100 do (
+    ylast := y;
+    y := y-(y*e^y-x)/(e^y+y*e^y);
+    if y == ylast then return y;
+    ylast2 := y;
+    y := y-(y*e^y-x)/(e^y+y*e^y);
+    if (y == ylast2 or y == ylast) then return y 
+  );
+  y
+)
+
+
+SetHelp ("LambertWm1", "functions", "The minus-one branch of the Lambert W function for real values between 
-1/e and 0");
+function LambertWm1(x) = (
+  if not IsReal(x) or x < -1/e or x >= 0.0 then (error ("LambertWm1: only defined for real x with -1/e <= x 
< 0");bailout);
+  if x < -0.116 then y := -8.22*x-4.23 else (l1:=ln(-x);l2:=ln(-ln(-x));y:=l1-l2+l2/l1);
+  for n=1 to 10 do (
+    ylast := y;
+    y := y-(y*e^y-x)/(e^y+y*e^y);
+    if (y == ylast) then return y 
+  );
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  y := y-(y*e^y-x)/(e^y+y*e^y);
+  for n=1 to 100 do (
+    ylast := y;
+    y := y-(y*e^y-x)/(e^y+y*e^y);
+    if y == ylast then return y;
+    ylast2 := y;
+    y := y-(y*e^y-x)/(e^y+y*e^y);
+    if (y == ylast2 or y == ylast) then return y 
+  );
+  y
+)
+
diff --git a/lib/library-strings.c b/lib/library-strings.c
index 2aefb17..0f5d62d 100644
--- a/lib/library-strings.c
+++ b/lib/library-strings.c
@@ -216,6 +216,8 @@ char *fake = N_("Returns 1 if and only if all elements are zero");
 char *fake = N_("The error function, 2/sqrt(pi) * int_0^x e^(-t^2) dt");
 char *fake = N_("Fejer kernel of order n");
 char *fake = N_("Returns 1 if and only if all elements are equal");
+char *fake = N_("Principal branch of the Lambert W function for real values greater than or equal to -1/e");
+char *fake = N_("The minus-one branch of the Lambert W function for real values between -1/e and 0");
 char *fake = N_("Find the first value where f(x)=0");
 char *fake = N_("Moebius mapping of the disk to itself mapping a to 0");
 char *fake = N_("Moebius mapping using the cross ratio taking z2,z3,z4 to 1,0, and infinity respectively");
@@ -235,6 +237,7 @@ char *fake = N_("Find root of a function using the bisection method to within TO
 char *fake = N_("Find root of a function using the method of false position to within TOL tolerance in up to 
N iterations.  f(a) and f(b) must have opposite signs.");
 char *fake = N_("Find root of a function using the Muller's method");
 char *fake = N_("Find root of a function using the secant method to within TOL tolerance in up to N 
iterations.  f(a) and f(b) must have opposite signs.");
+char *fake = N_("Run Newton's method on a function to attempt to find a root, returns after two successive 
values are within epsilon or after maxn tries (then returns null)");
 char *fake = N_("Find roots of a polynomial (given as vector of coefficients)");
 char *fake = N_("Find roots of a quartic polynomial (given as vector of coefficients)");
 char *fake = N_("Use classical non-adaptive Runge-Kutta of fourth order method to numerically solve 
y'=f(x,y) for initial x0,y0 going to x1 with n increments, returns y at x1");
diff --git a/src/geniustests.txt b/src/geniustests.txt
index 7294d77..9d49d36 100644
--- a/src/geniustests.txt
+++ b/src/geniustests.txt
@@ -1135,6 +1135,12 @@ KroneckerProduct([11,22;33,44],[1,2;3,4])                        
[11,22,22,44;33,44,66,88;33,66,44,88
 KroneckerProduct([-1,2],[1,2;3,4])                             [-1,-2,2,4;-3,-4,6,8]
 KroneckerProduct(null,[1,2;3,4])+0                             ((null)+0)
 KroneckerProduct([0],null)+0                                   ((null)+0)
+|NewtonsMethod (`(x)=x^2-10,`(x)=2*x,1.0,0.0001,100)-sqrt(10)|<0.0001  true
+|NewtonsMethod (`(x)=x^2-2,`(x)=2*x,1.0,0.00001,100)-sqrt(2)|<0.00001  true
+|NewtonsMethod (`(x)=x*e^x-0.5,`(x)=x*e^x+e^x,1.0,0.0001,100) - LambertW(0.5)|<0.0001  true
+|NewtonsMethod (`(x)=x*e^x+0.1,`(x)=x*e^x+e^x,1.0,0.0001,100) - LambertW(-0.1)|<0.0001 true
+|NewtonsMethod (`(x)=x*e^x-10,`(x)=x*e^x+e^x,1.0,0.0001,100) - LambertW(10)|<0.0001    true
+|NewtonsMethod (`(x)=x*e^x+0.1,`(x)=x*e^x+e^x,-2.0,0.0001,100) - LambertWm1(-0.1)|<0.0001      true
 load "nullspacetest.gel"                                       true
 load "longtest.gel"                                            true
 load "testprec.gel"                                            true
[Date Prev][Date Next] [Thread Prev][Thread Next] [Thread Index] [Date Index] [Author Index]