From: Jens Vollinga Date: Tue, 12 Oct 2004 13:32:11 +0000 (+0000) Subject: Added integral class. X-Git-Tag: release_1-3-0~15 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=c0ba2b239690e1b97aaecb860930b170a2bc1df2 Added integral class. --- diff --git a/ginac/Makefile.am b/ginac/Makefile.am index c7285c9d..2db8a602 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -5,10 +5,10 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \ constant.cpp ex.cpp expair.cpp expairseq.cpp exprseq.cpp fail.cpp \ fderivative.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp \ inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \ - lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp operators.cpp \ - power.cpp registrar.cpp relational.cpp remember.cpp pseries.cpp print.cpp \ - structure.cpp symbol.cpp symmetry.cpp tensor.cpp utils.cpp wildcard.cpp \ - input_parser.yy input_lexer.ll \ + integral.cpp lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp \ + operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \ + pseries.cpp print.cpp structure.cpp symbol.cpp symmetry.cpp tensor.cpp \ + utils.cpp wildcard.cpp input_parser.yy input_lexer.ll \ input_lexer.h remember.h tostring.h utils.h libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \ -release $(LT_RELEASE) @@ -16,9 +16,9 @@ ginacincludedir = $(includedir)/ginac ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h \ clifford.h color.h constant.h container.h ex.h expair.h expairseq.h exprseq.h \ fail.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h inifcns.h \ - lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h power.h print.h \ - pseries.h ptr.h registrar.h relational.h structure.h symbol.h symmetry.h \ - tensor.h tinfos.h version.h wildcard.h + integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \ + power.h print.h pseries.h ptr.h registrar.h relational.h structure.h \ + symbol.h symmetry.h tensor.h tinfos.h version.h wildcard.h AM_LFLAGS = -Pginac_yy -olex.yy.c AM_YFLAGS = -p ginac_yy -d EXTRA_DIST = function.pl input_parser.h version.h.in diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 98ba54d2..c37ed3e3 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -482,6 +482,20 @@ ex basic::evalm() const return map(map_evalm); } +/** Function object to be applied by basic::eval_integ(). */ +struct eval_integ_map_function : public map_function { + ex operator()(const ex & e) { return eval_integ(e); } +} map_eval_integ; + +/** Evaluate integrals, if result is known. */ +ex basic::eval_integ() const +{ + if (nops() == 0) + return *this; + else + return map(map_eval_integ); +} + /** Perform automatic symbolic evaluations on indexed expression that * contains this object as the base expression. */ ex basic::eval_indexed(const basic & i) const diff --git a/ginac/basic.h b/ginac/basic.h index a5b71a59..9f88d7b8 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -136,6 +136,7 @@ public: // only const functions please (may break reference counting) virtual ex eval(int level = 0) const; virtual ex evalf(int level = 0) const; virtual ex evalm() const; + virtual ex eval_integ() const; protected: virtual ex eval_ncmul(const exvector & v) const; public: diff --git a/ginac/ex.h b/ginac/ex.h index b274b49e..1c1e799d 100644 --- a/ginac/ex.h +++ b/ginac/ex.h @@ -117,6 +117,7 @@ public: ex evalf(int level = 0) const { return bp->evalf(level); } ex evalm() const { return bp->evalm(); } ex eval_ncmul(const exvector & v) const { return bp->eval_ncmul(v); } + ex eval_integ() const { return bp->eval_integ(); } // printing void print(const print_context & c, unsigned level = 0) const; @@ -807,6 +808,9 @@ inline ex evalf(const ex & thisex, int level = 0) inline ex evalm(const ex & thisex) { return thisex.evalm(); } +inline ex eval_integ(const ex & thisex) +{ return thisex.eval_integ(); } + inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1) { return thisex.diff(s, nth); } diff --git a/ginac/ginac.h b/ginac/ginac.h index 2ba0ad7f..54af6c41 100644 --- a/ginac/ginac.h +++ b/ginac/ginac.h @@ -34,6 +34,7 @@ #include "constant.h" #include "fail.h" +#include "integral.h" #include "lst.h" #include "matrix.h" #include "numeric.h" diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index d693373c..29dd86ae 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -36,6 +36,7 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "integral.h" namespace GiNaC { @@ -522,6 +523,13 @@ exvector power::get_free_indices() const return basis_indices; } +exvector integral::get_free_indices() const +{ + if (a.get_free_indices().size() || b.get_free_indices().size()) + throw (std::runtime_error("integral::get_free_indices: boundary values should not have free indices")); + return f.get_free_indices(); +} + /** Rename dummy indices in an expression. * * @param e Expression to work on diff --git a/ginac/integral.cpp b/ginac/integral.cpp new file mode 100644 index 00000000..467cf290 --- /dev/null +++ b/ginac/integral.cpp @@ -0,0 +1,441 @@ +/** @file integral.cpp + * + * Implementation of GiNaC's symbolic integral. */ + +/* + * GiNaC Copyright (C) 1999-2004 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 + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * 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 + */ + +#include "integral.h" +#include "numeric.h" +#include "symbol.h" +#include "add.h" +#include "mul.h" +#include "power.h" +#include "inifcns.h" +#include "wildcard.h" +#include "archive.h" +#include "registrar.h" +#include "utils.h" +#include "operators.h" +#include "relational.h" + +using namespace std; + +namespace GiNaC { + +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic, + print_func(&integral::do_print). + print_func(&integral::do_print_latex)) + + +////////// +// default constructor +////////// + +integral::integral() + : inherited(TINFO_integral), + x((new symbol())->setflag(status_flags::dynallocated)) +{} + +////////// +// other constructors +////////// + +// 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_) +{ + if (!is_a(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) +{ + n.find_ex("x", x, sym_lst); + n.find_ex("a", a, sym_lst); + n.find_ex("b", b, sym_lst); + n.find_ex("f", f, sym_lst); +} + +void integral::archive(archive_node & n) const +{ + inherited::archive(n); + n.add_ex("x", x); + n.add_ex("a", a); + n.add_ex("b", b); + n.add_ex("f", f); +} + +DEFAULT_UNARCHIVE(integral) + +////////// +// functions overriding virtual functions from base classes +////////// + +void integral::do_print(const print_context & c, unsigned level) const +{ + c.s << "integral("; + x.print(c); + c.s << ","; + a.print(c); + c.s << ","; + b.print(c); + c.s << ","; + f.print(c); + c.s << ")"; +} + +void integral::do_print_latex(const print_latex & c, unsigned level) const +{ + string varname = ex_to(x).get_name(); + if (level > precedence()) + c.s << "\\left("; + c.s << "\\int_{"; + a.print(c); + c.s << "}^{"; + b.print(c); + c.s << "} d"; + if (varname.size() > 1) + c.s << "\\," << varname << "\\:"; + else + c.s << varname << "\\,"; + f.print(c,precedence()); + if (level > precedence()) + c.s << "\\right)"; +} + +int integral::compare_same_type(const basic & other) const +{ + GINAC_ASSERT(is_exactly_a(other)); + const integral &o = static_cast(other); + + int cmpval = x.compare(o.x); + if (cmpval) + return cmpval; + cmpval = a.compare(o.a); + if (cmpval) + return cmpval; + cmpval = b.compare(o.b); + if (cmpval) + return cmpval; + return f.compare(o.f); +} + +ex integral::eval(int level) const +{ + if ((level==1) && (flags & status_flags::evaluated)) + return *this; + if (level == -max_recursion_level) + throw(std::runtime_error("max recursion level reached")); + + ex eintvar = (level==1) ? x : x.eval(level-1); + ex ea = (level==1) ? a : a.eval(level-1); + ex eb = (level==1) ? b : b.eval(level-1); + ex ef = (level==1) ? f : f.eval(level-1); + + if (!ef.has(eintvar) && !haswild(ef)) + return eb*ef-ea*ef; + + 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)) + return this->hold(); + return (new integral(eintvar, ea, eb, ef)) + ->setflag(status_flags::dynallocated | status_flags::evaluated); +} + +ex integral::evalf(int level) const +{ + ex ea; + ex eb; + ex ef; + + if (level==1) { + ea = a; + eb = b; + ef = f; + } else if (level == -max_recursion_level) { + throw(runtime_error("max recursion level reached")); + } else { + ea = a.evalf(level-1); + eb = b.evalf(level-1); + ef = f.evalf(level-1); + } + + // 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 { + 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)) + return *this; + else + return (new integral(x, ea, eb, ef)) + ->setflag(status_flags::dynallocated); +} + +int integral::max_integration_level = 15; +ex integral::relative_integration_error = power(10,-8).evalf(); + +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"); +} + +/** 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 + * the relative integration error. The function should evalf into a number + * 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) +{ + // use lookup table to be potentially much faster. + static exmap lookup; + static symbol ivar("ivar"); + ex lookupex = integral(ivar,a,b,f.subs(x==ivar)); + exmap::iterator emi = lookup.find(lookupex); + if (emi!=lookup.end()) + return emi->second; + + ex app = 0; + int i = 1; + exvector avec(integral::max_integration_level+1); + exvector hvec(integral::max_integration_level+1); + exvector favec(integral::max_integration_level+1); + exvector fbvec(integral::max_integration_level+1); + exvector fcvec(integral::max_integration_level+1); + exvector svec(integral::max_integration_level+1); + exvector errorvec(integral::max_integration_level+1); + vector lvec(integral::max_integration_level+1); + + avec[i] = a; + hvec[i] = (b-a)/2; + favec[i] = subsvalue(x, a, f); + fcvec[i] = subsvalue(x, a+hvec[i], f); + 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]; + + while (i>0) { + ex fd = subsvalue(x, avec[i]+hvec[i]/2, f); + ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f); + ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6; + ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6; + ex nu1 = avec[i]; + ex nu2 = favec[i]; + ex nu3 = fcvec[i]; + 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 nu7 = svec[i]; + int nu8 = lvec[i]; + --i; + if (abs(ex_to(s1+s2-nu7)) <= nu6) + app+=(s1+s2); + else { + if (nu8>=integral::max_integration_level) + throw runtime_error("max integration level reached"); + ++i; + avec[i] = nu1+nu5; + favec[i] = nu3; + fcvec[i] = fe; + fbvec[i] = nu4; + hvec[i] = nu5/2; + errorvec[i]=nu6/2; + svec[i] = s2; + lvec[i] = nu8+1; + ++i; + avec[i] = nu1; + favec[i] = nu2; + fcvec[i] = fd; + fbvec[i] = nu3; + hvec[i] = hvec[i-1]; + errorvec[i]=errorvec[i-1]; + svec[i] = s1; + lvec[i] = lvec[i-1]; + } + } + + lookup[lookupex]=app; + return app; +} + +int integral::degree(const ex & s) const +{ + return ((b-a)*f).degree(s); +} + +int integral::ldegree(const ex & s) const +{ + return ((b-a)*f).ldegree(s); +} + +ex integral::eval_ncmul(const exvector & v) const +{ + return f.eval_ncmul(v); +} + +size_t integral::nops() const +{ + return 4; +} + +ex integral::op(size_t i) const +{ + GINAC_ASSERT(i<4); + + switch(i) { + case(0): + return x; + case(1): + return a; + case(2): + return b; + case(3): + return f; + } +} + +ex & integral::let_op(size_t i) +{ + ensure_if_modifiable(); + switch(i) { + case(0): + return x; + case(1): + return a; + case(2): + return b; + case(3): + return f; + } +} + +ex integral::expand(unsigned options) const +{ + if (options==0 && (flags & status_flags::expanded)) + return *this; + + ex newa = a.expand(options); + ex newb = b.expand(options); + ex newf = f.expand(options); + + if (is_a(newf)) { + exvector v; + v.reserve(newf.nops()); + for (size_t i=0; i(newf)) { + ex prefactor = 1; + ex rest = 1; + for (size_t i=0; isetflag(status_flags::expanded); + return *this; + } + + const basic & newint = (new integral(x, newa, newb, newf)) + ->setflag(status_flags::dynallocated); + if (options == 0) + newint.setflag(status_flags::expanded); + return newint; +} + +ex integral::derivative(const symbol & s) const +{ if (s==x) + throw(logic_error("differentiation with respect to dummy variable")); + return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s)); +} + +unsigned integral::return_type() const +{ + return f.return_type(); +} + +unsigned integral::return_type_tinfo() const +{ + return f.return_type_tinfo(); +} + +ex integral::conjugate() const +{ + ex conja = a.conjugate(); + 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)) + return *this; + + return (new integral(x, conja, conjb, conjf)) + ->setflag(status_flags::dynallocated); +} + +ex integral::eval_integ() const +{ + if (!(flags & status_flags::expanded)) + return this->expand().eval_integ(); + + if (f==x) + return b*b/2-a*a/2; + if (is_a(f) && f.op(0)==x) { + if (f.op(1)==-1) + return log(b/a); + if (!f.op(1).has(x)) { + ex primit = power(x,f.op(1)+1)/(f.op(1)+1); + return primit.subs(x==b)-primit.subs(x==a); + } + } + + return *this; +} + +} // namespace GiNaC diff --git a/ginac/integral.h b/ginac/integral.h new file mode 100644 index 00000000..ef2e6626 --- /dev/null +++ b/ginac/integral.h @@ -0,0 +1,96 @@ +/** @file integral.h + * + * Interface to GiNaC's symbolic integral. */ + +/* + * GiNaC Copyright (C) 1999-2004 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 + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * 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 + */ + +#ifndef __GINAC_INTEGRAL_H__ +#define __GINAC_INTEGRAL_H__ + +#include "basic.h" +#include "ex.h" + +namespace GiNaC { + +/** Product of expressions. */ +class integral : public basic +{ + GINAC_DECLARE_REGISTERED_CLASS(integral, basic) + + // other constructors +public: + integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_); + + // functions overriding virtual functions from base classes +public: + unsigned precedence() const {return 45;} + ex eval(int level=0) const; + ex evalf(int level=0) const; + int degree(const ex & s) const; + int ldegree(const ex & s) const; + ex eval_ncmul(const exvector & v) const; + size_t nops() const; + ex op(size_t i) const; + ex & let_op(size_t i); + ex expand(unsigned options = 0) const; + exvector get_free_indices() const; + unsigned return_type() const; + unsigned return_type_tinfo() const; + ex conjugate() const; + ex eval_integ() const; +protected: + ex derivative(const symbol & s) const; + + // new virtual functions which can be overridden by derived classes + // none + + // non-virtual functions in this class +protected: + void do_print(const print_context & c, unsigned level) const; + void do_print_latex(const print_latex & c, unsigned level) const; +public: + static int max_integration_level; + static ex relative_integration_error; +private: + ex x; + ex a; + ex b; + ex f; +}; + +// utility functions + +/** Specialization of is_exactly_a(obj) for integral objects. */ +template<> inline bool is_exactly_a(const basic & obj) +{ + return obj.tinfo()==TINFO_integral; +} + + +GiNaC::ex adaptivesimpson( + const GiNaC::ex &x, + const GiNaC::ex &a, + const GiNaC::ex &b, + const GiNaC::ex &f, + const GiNaC::ex &error = integral::relative_integration_error +); + +} // namespace GiNaC + +#endif // ndef __GINAC_INTEGRAL_H__ diff --git a/ginac/tinfos.h b/ginac/tinfos.h index 37d681bd..f6b0551a 100644 --- a/ginac/tinfos.h +++ b/ginac/tinfos.h @@ -38,6 +38,7 @@ const unsigned TINFO_exprseq = 0x00030001U; const unsigned TINFO_function = 0x00031001U; const unsigned TINFO_fderivative = 0x00032001U; const unsigned TINFO_ncmul = 0x00031002U; +const unsigned TINFO_integral = 0x00033001U; const unsigned TINFO_lst = 0x00040001U; diff --git a/ginac/wildcard.cpp b/ginac/wildcard.cpp index 7ebc023e..1bcf2ac2 100644 --- a/ginac/wildcard.cpp +++ b/ginac/wildcard.cpp @@ -119,4 +119,14 @@ bool wildcard::match(const ex & pattern, lst & repl_lst) const return is_equal(ex_to(pattern)); } +bool haswild(const ex & x) +{ + if (is_a(x)) + return true; + for (int i=0; i