X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fintegral.cpp;h=016be2ca726f0f14a714e9c754d2a022dad26c55;hp=d2815c36f3b713e7b55a480cc12597ba38beac64;hb=2bf56ec52a7bed4ac3d02be8887b0287b5acd189;hpb=1af3517f2b664bb8f4d5b5fb066649492dc5e27e diff --git a/ginac/integral.cpp b/ginac/integral.cpp index d2815c36..016be2ca 100644 --- a/ginac/integral.cpp +++ b/ginac/integral.cpp @@ -3,7 +3,7 @@ * 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 @@ -17,7 +17,7 @@ * * 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" @@ -48,7 +48,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic, ////////// integral::integral() - : inherited(TINFO_integral), + : x((new symbol())->setflag(status_flags::dynallocated)) {} @@ -59,7 +59,7 @@ integral::integral() // 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(x)) { throw(std::invalid_argument("first argument of integral must be of type symbol")); @@ -70,8 +70,9 @@ integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_) // 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); @@ -87,8 +88,6 @@ void integral::archive(archive_node & n) const n.add_ex("f", f); } -DEFAULT_UNARCHIVE(integral) - ////////// // functions overriding virtual functions from base classes ////////// @@ -160,8 +159,8 @@ ex integral::eval(int level) const 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); @@ -186,16 +185,14 @@ 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(ea) && is_exactly_a(eb) - && is_exactly_a(ef.subs(x==12.34).evalf())) { - try { + // results after substituting a number for the integration variable. + if (is_exactly_a(ea) && is_exactly_a(eb) && + is_exactly_a(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)) @@ -203,16 +200,39 @@ ex integral::evalf(int level) const } 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(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 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 @@ -220,13 +240,21 @@ ex subsvalue(const ex & var, const ex & value, const ex & fun) * 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(a_in) ? a_in : a_in.evalf(); + ex b = is_exactly_a(b_in) ? b_in : b_in.evalf(); + if(!is_exactly_a(a) || !is_exactly_a(b)) + throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers."); + if(!is_exactly_a(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; @@ -248,7 +276,7 @@ ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const 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); @@ -261,7 +289,7 @@ ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const 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; @@ -291,7 +319,7 @@ ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const } } - lookup[lookupex]=app; + lookup[error_and_integral(error, lookupex)]=app; return app; } @@ -379,8 +407,8 @@ ex integral::expand(unsigned options) const 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; @@ -405,7 +433,7 @@ unsigned integral::return_type() const return f.return_type(); } -unsigned integral::return_type_tinfo() const +return_type_t integral::return_type_tinfo() const { return f.return_type_tinfo(); } @@ -416,8 +444,8 @@ ex integral::conjugate() const 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)) @@ -443,4 +471,5 @@ ex integral::eval_integ() const return *this; } +GINAC_BIND_UNARCHIVER(integral); } // namespace GiNaC