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
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;
#ifndef __GINAC_INIFCNS_H__
#define __GINAC_INIFCNS_H__
+#include "numeric.h"
#include "function.h"
#include "ex.h"
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)
{
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();
{"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},
git://www.ginac.de / ginac.git / commitdiff