* Added fsolve() numerical univariate real-valued function solver.

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 02d9e508a34b5b94d33f39162e4eaa1f6d45bf97..0d337cdb74275d10cb7961edb3d2c0a87f9e830a 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -417,6 +417,27 @@ x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
 Here we have made use of the @command{ginsh}-command @code{%} to pop the
 previously evaluated element from @command{ginsh}'s internal stack.
 
 Here we have made use of the @command{ginsh}-command @code{%} to pop the
 previously evaluated element from @command{ginsh}'s internal stack.
 

+Often, functions don't have roots in closed form.  Nevertheless, it's
+quite easy to compute a solution numerically, to arbitrary precision:
+
+@cindex fsolve
+@example
+> Digits=50:
+> fsolve(cos(x)-x,x,0,2);
+0.7390851332151606416553120876738734040134117589007574649658
+> f=exp(sin(x))-x:
+> X=fsolve(f,x,-10,10);
+2.2191071489137460325957851882042901681753665565320678854155
+> subs(f,x==X);
+-6.372367644529809108115521591070847222364418220770475144296E-58
+@end example
+
+Notice how the final result above differs slightly from zero by about
+@math{6*10^(-58)}.  This is because with 50 decimal digits precision the
+root cannot be represented more accurately than @code{X}.  Such
+inaccuracies are to be expected when computing with finite floating
+point values.
+

 If you ever wanted to convert units in C or C++ and found this is
 cumbersome, here is the solution.  Symbolic types can always be used as
 tags for different types of objects.  Converting from wrong units to the
 If you ever wanted to convert units in C or C++ and found this is
 cumbersome, here is the solution.  Symbolic types can always be used as
 tags for different types of objects.  Converting from wrong units to the

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index ebb384656ef9c7b7ab673e668f69f6fea072b39e..a9ad5ce5223ec84c2e0b0d4e55813d01a8d3a07a 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -663,6 +663,100 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
        return sollist;
 }
 
        return sollist;
 }
 

+//////////
+// Find real root of f(x) numerically
+//////////
+
+const numeric
+fsolve(const ex& f, const symbol& x, const numeric& x1, const numeric& x2)
+{
+       if (!x1.is_real() || !x2.is_real()) {
+               throw std::runtime_error("fsolve(): interval not bounded by real numbers");
+       }
+       if (x1==x2) {
+               throw std::runtime_error("fsolve(): vanishing interval");
+       }
+       // xx[0] == left interval limit, xx[1] == right interval limit.
+       // fx[0] == f(xx[0]), fx[1] == f(xx[1]).
+       // We keep the root bracketed: xx[0]<xx[1] and fx[0]*fx[1]<0.
+       numeric xx[2] = { x1<x2 ? x1 : x2,
+                         x1<x2 ? x2 : x1 };
+       const ex fx_[2] = { f.subs(x==xx[0]).evalf(),
+                           f.subs(x==xx[1]).evalf() };
+       if (!is_a<numeric>(fx_[0]) || !is_a<numeric>(fx_[1])) {
+               throw std::runtime_error("fsolve(): function does not evaluate numerically");
+       }
+       numeric fx[2] = { ex_to<numeric>(fx_[0]),
+                         ex_to<numeric>(fx_[1]) };
+       if (!fx[0].is_real() || !fx[1].is_real()) {
+               throw std::runtime_error("fsolve(): function evaluates to complex values at interval boundaries");
+       }
+       if (fx[0]*fx[1]>=0) {
+               throw std::runtime_error("fsolve(): function does not change sign at interval boundaries");
+       }
+
+       // The Newton-Raphson method has quadratic convergence!  Simply put, it
+       // replaces x with x-f(x)/f'(x) at each step.  -f/f' is the delta:
+       const ex ff = normal(-f/f.diff(x));
+       int side = 0;  // Start at left interval limit.
+       numeric xxprev;
+       numeric fxprev;
+       do {
+               xxprev = xx[side];
+               fxprev = fx[side];
+               xx[side] += ex_to<numeric>(ff.subs(x==xx[side]).evalf());
+               fx[side] = ex_to<numeric>(f.subs(x==xx[side]).evalf());
+               if ((side==0 && xx[0]<xxprev) || (side==1 && xx[1]>xxprev) || xx[0]>xx[1]) {
+                       // Oops, Newton-Raphson method shot out of the interval.
+                       // Restore, and try again with the other side instead!
+                       xx[side] = xxprev;
+                       fx[side] = fxprev;
+                       side = !side;
+                       xxprev = xx[side];
+                       fxprev = fx[side];
+                       xx[side] += ex_to<numeric>(ff.subs(x==xx[side]).evalf());
+                       fx[side] = ex_to<numeric>(f.subs(x==xx[side]).evalf());
+               }
+               if ((fx[side]<0 && fx[!side]<0) || (fx[side]>0 && fx[!side]>0)) {
+                       // Oops, the root isn't bracketed any more.
+                       // Restore, and perform a bisection!
+                       xx[side] = xxprev;
+                       fx[side] = fxprev;
+
+                       // Ah, the bisection! Bisections converge linearly. Unfortunately,
+                       // they occur pretty often when Newton-Raphson arrives at an x too
+                       // close to the result on one side of the interval and
+                       // f(x-f(x)/f'(x)) turns out to have the same sign as f(x) due to
+                       // precision errors! Recall that this function does not have a
+                       // precision goal as one of its arguments but instead relies on
+                       // x converging to a fixed point. We speed up the (safe but slow)
+                       // bisection method by mixing in a dash of the (unsafer but faster)
+                       // secant method: Instead of splitting the interval at the
+                       // arithmetic mean (bisection), we split it nearer to the root as
+                       // determined by the secant between the values xx[0] and xx[1].
+                       // Don't set the secant_weight to one because that could disturb
+                       // the convergence in some corner cases!
+                       static const double secant_weight = 0.96875;  // == 31/32 < 1
+                       numeric xxmid = (1-secant_weight)*0.5*(xx[0]+xx[1])
+                           + secant_weight*(xx[0]+fx[0]*(xx[0]-xx[1])/(fx[1]-fx[0]));
+                       numeric fxmid = ex_to<numeric>(f.subs(x==xxmid).evalf());
+                       if (fxmid.is_zero()) {
+                               // Luck strikes...
+                               return xxmid;
+                       }
+                       if ((fxmid<0 && fx[side]>0) || (fxmid>0 && fx[side]<0)) {
+                               side = !side;
+                       }
+                       xxprev = xx[side];
+                       fxprev = fx[side];
+                       xx[side] = xxmid;
+                       fx[side] = fxmid;
+               }
+       } while (xxprev!=xx[side]);
+       return xxprev;
+}
+
+

 /* Force inclusion of functions from inifcns_gamma and inifcns_zeta
  * for static lib (so ginsh will see them). */
 unsigned force_include_tgamma = tgamma_SERIAL::serial;
 /* Force inclusion of functions from inifcns_gamma and inifcns_zeta
  * for static lib (so ginsh will see them). */
 unsigned force_include_tgamma = tgamma_SERIAL::serial;

diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 9ce66942b9556cd524395d40a28667693cc32de1..05043819e9deeceb6a42d54cc03c5140238399e9 100644 (file)
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -23,6 +23,7 @@
 #ifndef __GINAC_INIFCNS_H__
 #define __GINAC_INIFCNS_H__
 
 #ifndef __GINAC_INIFCNS_H__
 #define __GINAC_INIFCNS_H__
 

+#include "numeric.h"

 #include "function.h"
 #include "ex.h"
 
 #include "function.h"
 #include "ex.h"
 


@@ -178,6 +179,17 @@ DECLARE_FUNCTION_1P(Order)
 
 ex lsolve(const ex &eqns, const ex &symbols, unsigned options = solve_algo::automatic);
 
 
 ex lsolve(const ex &eqns, const ex &symbols, unsigned options = solve_algo::automatic);
 

+/** Find a real root of real-valued function f(x) numerically within a given
+ *  interval. The function must change sign across interval. Uses Newton-
+ *  Raphson method combined with bisection in order to guarantee convergence.
+ *
+ *  @param f  Function f(x)
+ *  @param x  Symbol f(x)
+ *  @param x1  lower interval limit
+ *  @param x2  upper interval limit
+ *  @exception runtime_error (if interval is invalid). */
+const numeric fsolve(const ex& f, const symbol& x, const numeric& x1, const numeric& x2);
+

 /** Check whether a function is the Order (O(n)) function. */
 inline bool is_order_function(const ex & e)
 {
 /** Check whether a function is the Order (O(n)) function. */
 inline bool is_order_function(const ex & e)
 {

diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in
index bfdc1639268048425ff0d53c1a9ee60255bc54e5..d9a67e4b7382c48c5b5ef88b3a6568713fd90a85 100644 (file)
--- a/ginsh/ginsh.1.in
+++ b/ginsh/ginsh.1.in
@@ -293,6 +293,9 @@ detail here. Please refer to the GiNaC documentation.
 .BI find( expression ", " pattern )
 \- returns a list of all occurrences of a pattern in an expression
 .br
 .BI find( expression ", " pattern )
 \- returns a list of all occurrences of a pattern in an expression
 .br

+.BI fsolve( expression ", " symbol ", " number ", " number )
+\- numerically find root of a real-valued function within an interval
+.br

 .BI gcd( expression ", " expression )
 \- greatest common divisor
 .br
 .BI gcd( expression ", " expression )
 \- greatest common divisor
 .br

diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy
index 79ae94cee86b8d576eac497ec44e1b26d2705272..005ba9e9082b9dc460707b98324be37181baeac3 100644 (file)
--- a/ginsh/ginsh_parser.yy
+++ b/ginsh/ginsh_parser.yy
@@ -411,6 +411,14 @@ static ex f_find(const exprseq &e)
        return found;
 }
 
        return found;
 }
 

+static ex f_fsolve(const exprseq &e)
+{
+       CHECK_ARG(1, symbol, fsolve);
+       CHECK_ARG(2, numeric, fsolve);
+       CHECK_ARG(3, numeric, fsolve);
+       return fsolve(e[0], ex_to<symbol>(e[1]), ex_to<numeric>(e[2]), ex_to<numeric>(e[3]));
+}
+

 static ex f_integer_content(const exprseq &e)
 {
        return e[0].expand().integer_content();
 static ex f_integer_content(const exprseq &e)
 {
        return e[0].expand().integer_content();


@@ -588,6 +596,7 @@ static const fcn_init builtin_fcns[] = {
        {"eval_integ", f_eval_integ, 1},
        {"expand", f_expand, 1},
        {"find", f_find, 2},
        {"eval_integ", f_eval_integ, 1},
        {"expand", f_expand, 1},
        {"find", f_find, 2},

+       {"fsolve", f_fsolve, 4},

        {"gcd", f_gcd, 2},
        {"has", f_has, 2},
        {"integer_content", f_integer_content, 1},
        {"gcd", f_gcd, 2},
        {"has", f_has, 2},
        {"integer_content", f_integer_content, 1},