* Implementation of GiNaC's symbolic integral. */
/*
- * GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include "integral.h"
//////////
integral::integral()
- : inherited(TINFO_integral),
- x((new symbol())->setflag(status_flags::dynallocated))
+ : x(dynallocate<symbol>())
{}
//////////
// public
integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
- : inherited(TINFO_integral), x(x_), a(a_), b(b_), f(f_)
+ : x(x_), a(a_), b(b_), f(f_)
{
if (!is_a<symbol>(x)) {
throw(std::invalid_argument("first argument of integral must be of type symbol"));
// archiving
//////////
-integral::integral(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
+void integral::read_archive(const archive_node& n, lst& sym_lst)
{
+ inherited::read_archive(n, sym_lst);
n.find_ex("x", x, sym_lst);
n.find_ex("a", a, sym_lst);
n.find_ex("b", b, sym_lst);
n.add_ex("f", f);
}
-DEFAULT_UNARCHIVE(integral)
-
//////////
// functions overriding virtual functions from base classes
//////////
if (ea==eb)
return _ex0;
- if (are_ex_trivially_equal(eintvar,x) && are_ex_trivially_equal(ea,a)
- && are_ex_trivially_equal(eb,b) && are_ex_trivially_equal(ef,f))
+ if (are_ex_trivially_equal(eintvar,x) && are_ex_trivially_equal(ea,a) &&
+ are_ex_trivially_equal(eb,b) && are_ex_trivially_equal(ef,f))
return this->hold();
- return (new integral(eintvar, ea, eb, ef))
- ->setflag(status_flags::dynallocated | status_flags::evaluated);
+ return dynallocate<integral>(eintvar, ea, eb, ef).setflag(status_flags::evaluated);
}
ex integral::evalf(int level) const
}
// 12.34 is just an arbitrary number used to check whether a number
- // results after subsituting a number for the integration variable.
- if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb)
- && is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
- try {
+ // results after substituting a number for the integration variable.
+ if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb) &&
+ is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
return adaptivesimpson(x, ea, eb, ef);
- } catch (runtime_error &rte) {}
}
- if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb)
- && are_ex_trivially_equal(f, ef))
+ if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb) &&
+ are_ex_trivially_equal(f, ef))
return *this;
else
- return (new integral(x, ea, eb, ef))
- ->setflag(status_flags::dynallocated);
+ return dynallocate<integral>(x, ea, eb, ef);
}
int integral::max_integration_level = 15;
-ex integral::relative_integration_error = power(10,-8).evalf();
+ex integral::relative_integration_error = 1e-8;
ex subsvalue(const ex & var, const ex & value, const ex & fun)
{
ex result = fun.subs(var==value).evalf();
if (is_a<numeric>(result))
return result;
- throw logic_error("integrant does not evaluate to numeric");
+ throw logic_error("integrand does not evaluate to numeric");
}
+struct error_and_integral
+{
+ error_and_integral(const ex &err, const ex &integ)
+ :error(err), integral(integ){}
+ ex error;
+ ex integral;
+};
+
+struct error_and_integral_is_less
+{
+ bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
+ {
+ int c = e1.integral.compare(e2.integral);
+ if(c < 0)
+ return true;
+ if(c > 0)
+ return false;
+ return ex_is_less()(e1.error, e2.error);
+ }
+};
+
+typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
+
/** Numeric integration routine based upon the "Adaptive Quadrature" one
* in "Numerical Analysis" by Burden and Faires. Parameters are integration
* variable, left boundary, right boundary, function to be integrated and
* after substituting the integration variable by a number. Another thing
* to note is that this implementation is no good at integrating functions
* with discontinuities. */
-ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const ex & error)
+ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
{
- // use lookup table to be potentially much faster.
- static exmap lookup;
+ // Check whether boundaries and error are numbers.
+ ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
+ ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
+ if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
+ throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
+ if(!is_exactly_a<numeric>(error))
+ throw std::runtime_error("For numerical integration the error should be a number.");
+
+ // Use lookup table to be potentially much faster.
+ static lookup_map lookup;
static symbol ivar("ivar");
ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
- exmap::iterator emi = lookup.find(lookupex);
+ lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
if (emi!=lookup.end())
return emi->second;
fbvec[i] = subsvalue(x, b, f);
svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
lvec[i] = 1;
- errorvec[i] = integral::relative_integration_error*svec[i];
+ errorvec[i] = error*abs(svec[i]);
while (i>0) {
ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
ex nu4 = fbvec[i];
ex nu5 = hvec[i];
// hopefully prevents a crash if the function is zero sometimes.
- ex nu6 = max(errorvec[i], (s1+s2)*integral::relative_integration_error);
+ ex nu6 = max(errorvec[i], abs(s1+s2)*error);
ex nu7 = svec[i];
int nu8 = lvec[i];
--i;
}
}
- lookup[lookupex]=app;
+ lookup[error_and_integral(error, lookupex)]=app;
return app;
}
return (prefactor*integral(x, newa, newb, rest)).expand(options);
}
- if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb)
- && are_ex_trivially_equal(f, newf)) {
+ if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb) &&
+ are_ex_trivially_equal(f, newf)) {
if (options==0)
this->setflag(status_flags::expanded);
return *this;
}
- const basic & newint = (new integral(x, newa, newb, newf))
- ->setflag(status_flags::dynallocated);
+ const integral & newint = dynallocate<integral>(x, newa, newb, newf);
if (options == 0)
newint.setflag(status_flags::expanded);
return newint;
return f.return_type();
}
-unsigned integral::return_type_tinfo() const
+return_type_t integral::return_type_tinfo() const
{
return f.return_type_tinfo();
}
ex conjb = b.conjugate();
ex conjf = f.conjugate().subs(x.conjugate()==x);
- if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb)
- && are_ex_trivially_equal(f, conjf))
+ if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb) &&
+ are_ex_trivially_equal(f, conjf))
return *this;
- return (new integral(x, conja, conjb, conjf))
- ->setflag(status_flags::dynallocated);
+ return dynallocate<integral>(x, conja, conjb, conjf);
}
ex integral::eval_integ() const
return *this;
}
+GINAC_BIND_UNARCHIVER(integral);
} // namespace GiNaC
git://www.ginac.de / ginac.git / blobdiff