From: Richard Kreckel Date: Mon, 12 Sep 2005 02:33:12 +0000 (+0000) Subject: * Added fsolve() numerical univariate real-valued function solver. X-Git-Tag: release_1-4-0~155 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=4d59c02d51fbf50ff24d616b00296aa4cbb1ea5e;ds=sidebyside * Added fsolve() numerical univariate real-valued function solver. --- diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 02d9e508..0d337cdb 100644 --- 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. +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 diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index ebb38465..a9ad5ce5 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -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](fx_[0]) || !is_a(fx_[1])) { + throw std::runtime_error("fsolve(): function does not evaluate numerically"); + } + numeric fx[2] = { ex_to(fx_[0]), + ex_to(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(ff.subs(x==xx[side]).evalf()); + fx[side] = ex_to(f.subs(x==xx[side]).evalf()); + if ((side==0 && xx[0]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(ff.subs(x==xx[side]).evalf()); + fx[side] = ex_to(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(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; diff --git a/ginac/inifcns.h b/ginac/inifcns.h index 9ce66942..05043819 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -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) { diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in index bfdc1639..d9a67e4b 100644 --- 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 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 diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index 79ae94ce..005ba9e9 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -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(e[1]), ex_to(e[2]), ex_to(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},