* Added fsolve() numerical univariate real-valued function solver.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Mon, 12 Sep 2005 02:33:12 +0000 (02:33 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Mon, 12 Sep 2005 02:33:12 +0000 (02:33 +0000)
doc/tutorial/ginac.texi
ginac/inifcns.cpp
ginac/inifcns.h
ginsh/ginsh.1.in
ginsh/ginsh_parser.yy

index 02d9e508a34b5b94d33f39162e4eaa1f6d45bf97..0d337cdb74275d10cb7961edb3d2c0a87f9e830a 100644 (file)
@@ -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.
 
+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
index ebb384656ef9c7b7ab673e668f69f6fea072b39e..a9ad5ce5223ec84c2e0b0d4e55813d01a8d3a07a 100644 (file)
@@ -663,6 +663,100 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
        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;
index 9ce66942b9556cd524395d40a28667693cc32de1..05043819e9deeceb6a42d54cc03c5140238399e9 100644 (file)
@@ -23,6 +23,7 @@
 #ifndef __GINAC_INIFCNS_H__
 #define __GINAC_INIFCNS_H__
 
+#include "numeric.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);
 
+/** 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)
 {
index bfdc1639268048425ff0d53c1a9ee60255bc54e5..d9a67e4b7382c48c5b5ef88b3a6568713fd90a85 100644 (file)
@@ -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 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
index 79ae94cee86b8d576eac497ec44e1b26d2705272..005ba9e9082b9dc460707b98324be37181baeac3 100644 (file)
@@ -411,6 +411,14 @@ static ex f_find(const exprseq &e)
        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();
@@ -588,6 +596,7 @@ static const fcn_init builtin_fcns[] = {
        {"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},