* 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 02d9e50..0d337cd 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 ebb3846..a9ad5ce 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 9ce6694..0504381 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 bfdc163..d9a67e4 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 79ae94c..005ba9e 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},