author Richard Kreckel Mon, 12 Sep 2005 02:33:12 +0000 (02:33 +0000) committer Richard Kreckel Mon, 12 Sep 2005 02:33:12 +0000 (02:33 +0000)
 doc/tutorial/ginac.texi patch | blob | history ginac/inifcns.cpp patch | blob | history ginac/inifcns.h patch | blob | history ginsh/ginsh.1.in patch | blob | history ginsh/ginsh_parser.yy patch | blob | history

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
@@ -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 == left interval limit, xx == right interval limit.
+       // fx == f(xx), fx == f(xx).
+       // We keep the root bracketed: xx<xx and fx*fx<0.
+       numeric xx = { x1<x2 ? x1 : x2,
+                         x1<x2 ? x2 : x1 };
+       const ex fx_ = { f.subs(x==xx).evalf(),
+                           f.subs(x==xx).evalf() };
+       if (!is_a<numeric>(fx_) || !is_a<numeric>(fx_)) {
+               throw std::runtime_error("fsolve(): function does not evaluate numerically");
+       }
+       numeric fx = { ex_to<numeric>(fx_),
+                         ex_to<numeric>(fx_) };
+       if (!fx.is_real() || !fx.is_real()) {
+               throw std::runtime_error("fsolve(): function evaluates to complex values at interval boundaries");
+       }
+       if (fx*fx>=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<xxprev) || (side==1 && xx>xxprev) || xx>xx) {
+                       // 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 and xx.
+                       // 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+xx)
+                           + secant_weight*(xx+fx*(xx-xx)/(fx-fx));
+                       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, ex_to<symbol>(e), ex_to<numeric>(e), ex_to<numeric>(e));
+}
+
static ex f_integer_content(const exprseq &e)
{
return e.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},