-ex lsolve(ex const &eqns, ex const &symbols)
-{
- // solve a system of linear equations
- if (eqns.info(info_flags::relation_equal)) {
- if (!symbols.info(info_flags::symbol)) {
- throw(std::invalid_argument("lsolve: 2nd argument must be a symbol"));
- }
- ex sol=lsolve(lst(eqns),lst(symbols));
-
- GINAC_ASSERT(sol.nops()==1);
- GINAC_ASSERT(is_ex_exactly_of_type(sol.op(0),relational));
-
- return sol.op(0).op(1); // return rhs of first solution
- }
-
- // syntax checks
- if (!eqns.info(info_flags::list)) {
- throw(std::invalid_argument("lsolve: 1st argument must be a list"));
- }
- for (unsigned i=0; i<eqns.nops(); i++) {
- if (!eqns.op(i).info(info_flags::relation_equal)) {
- throw(std::invalid_argument("lsolve: 1st argument must be a list of equations"));
- }
- }
- if (!symbols.info(info_flags::list)) {
- throw(std::invalid_argument("lsolve: 2nd argument must be a list"));
- }
- for (unsigned i=0; i<symbols.nops(); i++) {
- if (!symbols.op(i).info(info_flags::symbol)) {
- throw(std::invalid_argument("lsolve: 2nd argument must be a list of symbols"));
- }
- }
-
- // build matrix from equation system
- matrix sys(eqns.nops(),symbols.nops());
- matrix rhs(eqns.nops(),1);
- matrix vars(symbols.nops(),1);
-
- for (unsigned r=0; r<eqns.nops(); r++) {
- ex eq=eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
- ex linpart=eq;
- for (unsigned c=0; c<symbols.nops(); c++) {
- ex co=eq.coeff(ex_to_symbol(symbols.op(c)),1);
- linpart -= co*symbols.op(c);
- sys.set(r,c,co);
- }
- linpart=linpart.expand();
- rhs.set(r,0,-linpart);
- }
-
- // test if system is linear and fill vars matrix
- for (unsigned i=0; i<symbols.nops(); i++) {
- vars.set(i,0,symbols.op(i));
- if (sys.has(symbols.op(i))) {
- throw(std::logic_error("lsolve: system is not linear"));
- }
- if (rhs.has(symbols.op(i))) {
- throw(std::logic_error("lsolve: system is not linear"));
- }
- }
-
- //matrix solution=sys.solve(rhs);
- matrix solution;
- try {
- solution=sys.fraction_free_elim(vars,rhs);
- } catch (runtime_error const & e) {
- // probably singular matrix (or other error)
- // return empty solution list
- // cerr << e.what() << endl;
- return lst();
- }
-
- // return a list of equations
- if (solution.cols()!=1) {
- throw(std::runtime_error("lsolve: strange number of columns returned from matrix::solve"));
- }
- if (solution.rows()!=symbols.nops()) {
- cout << "symbols.nops()=" << symbols.nops() << endl;
- cout << "solution.rows()=" << solution.rows() << endl;
- throw(std::runtime_error("lsolve: strange number of rows returned from matrix::solve"));
- }
-
- // return list of the form lst(var1==sol1,var2==sol2,...)
- lst sollist;
- for (unsigned i=0; i<symbols.nops(); i++) {
- sollist.append(symbols.op(i)==solution(i,0));
- }
-
- return sollist;
-}
-
-/** non-commutative power. */
-ex ncpower(ex const &basis, unsigned exponent)
-{
- if (exponent==0) {
- return _ex1();
- }
-
- exvector v;
- v.reserve(exponent);
- for (unsigned i=0; i<exponent; ++i) {
- v.push_back(basis);
- }
-
- return ncmul(v,1);
-}
-
-/** Force inclusion of functions from initcns_gamma and inifcns_zeta
- * for static lib (so ginsh will see them). */
-unsigned force_include_gamma = function_index_gamma;
-unsigned force_include_zeta1 = function_index_zeta1;
-
-#ifndef NO_GINAC_NAMESPACE
+ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
+{
+ // solve a system of linear equations
+ if (eqns.info(info_flags::relation_equal)) {
+ if (!symbols.info(info_flags::symbol))
+ throw(std::invalid_argument("lsolve(): 2nd argument must be a symbol"));
+ const ex sol = lsolve(lst(eqns),lst(symbols));
+
+ GINAC_ASSERT(sol.nops()==1);
+ GINAC_ASSERT(is_exactly_a<relational>(sol.op(0)));
+
+ return sol.op(0).op(1); // return rhs of first solution
+ }
+
+ // syntax checks
+ if (!eqns.info(info_flags::list)) {
+ throw(std::invalid_argument("lsolve(): 1st argument must be a list"));
+ }
+ for (size_t i=0; i<eqns.nops(); i++) {
+ if (!eqns.op(i).info(info_flags::relation_equal)) {
+ throw(std::invalid_argument("lsolve(): 1st argument must be a list of equations"));
+ }
+ }
+ if (!symbols.info(info_flags::list)) {
+ throw(std::invalid_argument("lsolve(): 2nd argument must be a list"));
+ }
+ for (size_t i=0; i<symbols.nops(); i++) {
+ if (!symbols.op(i).info(info_flags::symbol)) {
+ throw(std::invalid_argument("lsolve(): 2nd argument must be a list of symbols"));
+ }
+ }
+
+ // build matrix from equation system
+ matrix sys(eqns.nops(),symbols.nops());
+ matrix rhs(eqns.nops(),1);
+ matrix vars(symbols.nops(),1);
+
+ for (size_t r=0; r<eqns.nops(); r++) {
+ const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
+ ex linpart = eq;
+ for (size_t c=0; c<symbols.nops(); c++) {
+ const ex co = eq.coeff(ex_to<symbol>(symbols.op(c)),1);
+ linpart -= co*symbols.op(c);
+ sys(r,c) = co;
+ }
+ linpart = linpart.expand();
+ rhs(r,0) = -linpart;
+ }
+
+ // test if system is linear and fill vars matrix
+ for (size_t i=0; i<symbols.nops(); i++) {
+ vars(i,0) = symbols.op(i);
+ if (sys.has(symbols.op(i)))
+ throw(std::logic_error("lsolve: system is not linear"));
+ if (rhs.has(symbols.op(i)))
+ throw(std::logic_error("lsolve: system is not linear"));
+ }
+
+ matrix solution;
+ try {
+ solution = sys.solve(vars,rhs,options);
+ } catch (const std::runtime_error & e) {
+ // Probably singular matrix or otherwise overdetermined system:
+ // It is consistent to return an empty list
+ return lst();
+ }
+ GINAC_ASSERT(solution.cols()==1);
+ GINAC_ASSERT(solution.rows()==symbols.nops());
+
+ // return list of equations of the form lst(var1==sol1,var2==sol2,...)
+ lst sollist;
+ for (size_t i=0; i<symbols.nops(); i++)
+ sollist.append(symbols.op(i)==solution(i,0));
+
+ 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;
+unsigned force_include_zeta1 = zeta1_SERIAL::serial;
+