3 * Implementation of GiNaC's symbolic integral. */
6 * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
32 #include "registrar.h"
34 #include "operators.h"
35 #include "relational.h"
41 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic,
42 print_func<print_dflt>(&integral::do_print).
43 print_func<print_latex>(&integral::do_print_latex))
47 // default constructor
51 : inherited(TINFO_integral),
52 x((new symbol())->setflag(status_flags::dynallocated))
61 integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
62 : inherited(TINFO_integral), x(x_), a(a_), b(b_), f(f_)
64 if (!is_a<symbol>(x)) {
65 throw(std::invalid_argument("first argument of integral must be of type symbol"));
73 integral::integral(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
75 n.find_ex("x", x, sym_lst);
76 n.find_ex("a", a, sym_lst);
77 n.find_ex("b", b, sym_lst);
78 n.find_ex("f", f, sym_lst);
81 void integral::archive(archive_node & n) const
83 inherited::archive(n);
90 DEFAULT_UNARCHIVE(integral)
93 // functions overriding virtual functions from base classes
96 void integral::do_print(const print_context & c, unsigned level) const
109 void integral::do_print_latex(const print_latex & c, unsigned level) const
111 string varname = ex_to<symbol>(x).get_name();
112 if (level > precedence())
119 if (varname.size() > 1)
120 c.s << "\\," << varname << "\\:";
122 c.s << varname << "\\,";
123 f.print(c,precedence());
124 if (level > precedence())
128 int integral::compare_same_type(const basic & other) const
130 GINAC_ASSERT(is_exactly_a<integral>(other));
131 const integral &o = static_cast<const integral &>(other);
133 int cmpval = x.compare(o.x);
136 cmpval = a.compare(o.a);
139 cmpval = b.compare(o.b);
142 return f.compare(o.f);
145 ex integral::eval(int level) const
147 if ((level==1) && (flags & status_flags::evaluated))
149 if (level == -max_recursion_level)
150 throw(std::runtime_error("max recursion level reached"));
152 ex eintvar = (level==1) ? x : x.eval(level-1);
153 ex ea = (level==1) ? a : a.eval(level-1);
154 ex eb = (level==1) ? b : b.eval(level-1);
155 ex ef = (level==1) ? f : f.eval(level-1);
157 if (!ef.has(eintvar) && !haswild(ef))
163 if (are_ex_trivially_equal(eintvar,x) && are_ex_trivially_equal(ea,a)
164 && are_ex_trivially_equal(eb,b) && are_ex_trivially_equal(ef,f))
166 return (new integral(eintvar, ea, eb, ef))
167 ->setflag(status_flags::dynallocated | status_flags::evaluated);
170 ex integral::evalf(int level) const
180 } else if (level == -max_recursion_level) {
181 throw(runtime_error("max recursion level reached"));
183 ea = a.evalf(level-1);
184 eb = b.evalf(level-1);
185 ef = f.evalf(level-1);
188 // 12.34 is just an arbitrary number used to check whether a number
189 // results after subsituting a number for the integration variable.
190 if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb)
191 && is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
193 return adaptivesimpson(x, ea, eb, ef);
194 } catch (runtime_error &rte) {}
197 if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb)
198 && are_ex_trivially_equal(f, ef))
201 return (new integral(x, ea, eb, ef))
202 ->setflag(status_flags::dynallocated);
205 int integral::max_integration_level = 15;
206 ex integral::relative_integration_error = 1e-8;
208 ex subsvalue(const ex & var, const ex & value, const ex & fun)
210 ex result = fun.subs(var==value).evalf();
211 if (is_a<numeric>(result))
213 throw logic_error("integrand does not evaluate to numeric");
216 struct error_and_integral
218 error_and_integral(const ex &err, const ex &integ)
219 :error(err), integral(integ){}
224 struct error_and_integral_is_less
226 bool operator()(const error_and_integral &e1,const error_and_integral &e2)
228 int c = e1.integral.compare(e2.integral);
233 return ex_is_less()(e1.error, e2.error);
237 typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
239 /** Numeric integration routine based upon the "Adaptive Quadrature" one
240 * in "Numerical Analysis" by Burden and Faires. Parameters are integration
241 * variable, left boundary, right boundary, function to be integrated and
242 * the relative integration error. The function should evalf into a number
243 * after substituting the integration variable by a number. Another thing
244 * to note is that this implementation is no good at integrating functions
245 * with discontinuities. */
246 ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
248 // Check whether boundaries and error are numbers.
249 ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
250 ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
251 if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
252 throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
253 if(!is_exactly_a<numeric>(error))
254 throw std::runtime_error("For numerical integration the error should be a number.");
256 // Use lookup table to be potentially much faster.
257 static lookup_map lookup;
258 static symbol ivar("ivar");
259 ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
260 lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
261 if (emi!=lookup.end())
266 exvector avec(integral::max_integration_level+1);
267 exvector hvec(integral::max_integration_level+1);
268 exvector favec(integral::max_integration_level+1);
269 exvector fbvec(integral::max_integration_level+1);
270 exvector fcvec(integral::max_integration_level+1);
271 exvector svec(integral::max_integration_level+1);
272 exvector errorvec(integral::max_integration_level+1);
273 vector<int> lvec(integral::max_integration_level+1);
277 favec[i] = subsvalue(x, a, f);
278 fcvec[i] = subsvalue(x, a+hvec[i], f);
279 fbvec[i] = subsvalue(x, b, f);
280 svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
282 errorvec[i] = error*svec[i];
285 ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
286 ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f);
287 ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6;
288 ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6;
294 // hopefully prevents a crash if the function is zero sometimes.
295 ex nu6 = max(errorvec[i], (s1+s2)*error);
299 if (abs(ex_to<numeric>(s1+s2-nu7)) <= nu6)
302 if (nu8>=integral::max_integration_level)
303 throw runtime_error("max integration level reached");
319 errorvec[i]=errorvec[i-1];
325 lookup[error_and_integral(error, lookupex)]=app;
329 int integral::degree(const ex & s) const
331 return ((b-a)*f).degree(s);
334 int integral::ldegree(const ex & s) const
336 return ((b-a)*f).ldegree(s);
339 ex integral::eval_ncmul(const exvector & v) const
341 return f.eval_ncmul(v);
344 size_t integral::nops() const
349 ex integral::op(size_t i) const
363 throw (std::out_of_range("integral::op() out of range"));
367 ex & integral::let_op(size_t i)
369 ensure_if_modifiable();
380 throw (std::out_of_range("integral::let_op() out of range"));
384 ex integral::expand(unsigned options) const
386 if (options==0 && (flags & status_flags::expanded))
389 ex newa = a.expand(options);
390 ex newb = b.expand(options);
391 ex newf = f.expand(options);
393 if (is_a<add>(newf)) {
395 v.reserve(newf.nops());
396 for (size_t i=0; i<newf.nops(); ++i)
397 v.push_back(integral(x, newa, newb, newf.op(i)).expand(options));
398 return ex(add(v)).expand(options);
401 if (is_a<mul>(newf)) {
404 for (size_t i=0; i<newf.nops(); ++i)
405 if (newf.op(i).has(x))
408 prefactor *= newf.op(i);
410 return (prefactor*integral(x, newa, newb, rest)).expand(options);
413 if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb)
414 && are_ex_trivially_equal(f, newf)) {
416 this->setflag(status_flags::expanded);
420 const basic & newint = (new integral(x, newa, newb, newf))
421 ->setflag(status_flags::dynallocated);
423 newint.setflag(status_flags::expanded);
427 ex integral::derivative(const symbol & s) const
430 throw(logic_error("differentiation with respect to dummy variable"));
431 return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s));
434 unsigned integral::return_type() const
436 return f.return_type();
439 unsigned integral::return_type_tinfo() const
441 return f.return_type_tinfo();
444 ex integral::conjugate() const
446 ex conja = a.conjugate();
447 ex conjb = b.conjugate();
448 ex conjf = f.conjugate().subs(x.conjugate()==x);
450 if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb)
451 && are_ex_trivially_equal(f, conjf))
454 return (new integral(x, conja, conjb, conjf))
455 ->setflag(status_flags::dynallocated);
458 ex integral::eval_integ() const
460 if (!(flags & status_flags::expanded))
461 return this->expand().eval_integ();
465 if (is_a<power>(f) && f.op(0)==x) {
468 if (!f.op(1).has(x)) {
469 ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
470 return primit.subs(x==b)-primit.subs(x==a);