GiNaC 1.8.7
integral.cpp
Go to the documentation of this file.
1
5/*
6 * GiNaC Copyright (C) 1999-2023 Johannes Gutenberg University Mainz, Germany
7 *
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.
12 *
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.
17 *
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
21 */
22
23#include "integral.h"
24#include "numeric.h"
25#include "symbol.h"
26#include "add.h"
27#include "mul.h"
28#include "power.h"
29#include "inifcns.h"
30#include "wildcard.h"
31#include "archive.h"
32#include "registrar.h"
33#include "utils.h"
34#include "operators.h"
35#include "relational.h"
36
37using namespace std;
38
39namespace GiNaC {
40
43 print_func<print_python>(&integral::do_print).
44 print_func<print_latex>(&integral::do_print_latex))
45
46
47
48// default constructor
50
52 : x(dynallocate<symbol>())
53{}
54
56// other constructors
58
59// public
60
61integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
62 : x(x_), a(a_), b(b_), f(f_)
63{
64 if (!is_a<symbol>(x)) {
65 throw(std::invalid_argument("first argument of integral must be of type symbol"));
66 }
67}
68
70// archiving
72
74{
75 inherited::read_archive(n, sym_lst);
76 n.find_ex("x", x, sym_lst);
77 n.find_ex("a", a, sym_lst);
78 n.find_ex("b", b, sym_lst);
79 n.find_ex("f", f, sym_lst);
80}
81
83{
84 inherited::archive(n);
85 n.add_ex("x", x);
86 n.add_ex("a", a);
87 n.add_ex("b", b);
88 n.add_ex("f", f);
89}
90
92// functions overriding virtual functions from base classes
94
95void integral::do_print(const print_context & c, unsigned level) const
96{
97 c.s << "integral(";
98 x.print(c);
99 c.s << ",";
100 a.print(c);
101 c.s << ",";
102 b.print(c);
103 c.s << ",";
104 f.print(c);
105 c.s << ")";
106}
107
108void integral::do_print_latex(const print_latex & c, unsigned level) const
109{
110 string varname = ex_to<symbol>(x).get_name();
111 if (level > precedence())
112 c.s << "\\left(";
113 c.s << "\\int_{";
114 a.print(c);
115 c.s << "}^{";
116 b.print(c);
117 c.s << "} d";
118 if (varname.size() > 1)
119 c.s << "\\," << varname << "\\:";
120 else
121 c.s << varname << "\\,";
122 f.print(c,precedence());
123 if (level > precedence())
124 c.s << "\\right)";
125}
126
127int integral::compare_same_type(const basic & other) const
128{
129 GINAC_ASSERT(is_exactly_a<integral>(other));
130 const integral &o = static_cast<const integral &>(other);
131
132 int cmpval = x.compare(o.x);
133 if (cmpval)
134 return cmpval;
135 cmpval = a.compare(o.a);
136 if (cmpval)
137 return cmpval;
138 cmpval = b.compare(o.b);
139 if (cmpval)
140 return cmpval;
141 return f.compare(o.f);
142}
143
145{
147 return *this;
148
149 if (!f.has(x) && !haswild(f))
150 return b*f-a*f;
151
152 if (a==b)
153 return _ex0;
154
155 return this->hold();
156}
157
159{
160 ex ea = a.evalf();
161 ex eb = b.evalf();
162 ex ef = f.evalf();
163
164 // 12.34 is just an arbitrary number used to check whether a number
165 // results after substituting a number for the integration variable.
166 if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb) &&
167 is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
168 return adaptivesimpson(x, ea, eb, ef);
169 }
170
173 return *this;
174 else
175 return dynallocate<integral>(x, ea, eb, ef);
176}
177
180
181ex subsvalue(const ex & var, const ex & value, const ex & fun)
182{
183 ex result = fun.subs(var==value).evalf();
184 if (is_a<numeric>(result))
185 return result;
186 throw logic_error("integrand does not evaluate to numeric");
187}
188
190{
191 error_and_integral(const ex &err, const ex &integ)
192 :error(err), integral(integ){}
195};
196
198{
199 bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
200 {
201 int c = e1.integral.compare(e2.integral);
202 if(c < 0)
203 return true;
204 if(c > 0)
205 return false;
206 return ex_is_less()(e1.error, e2.error);
207 }
208};
209
210typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
211
219ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
220{
221 // Check whether boundaries and error are numbers.
222 ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
223 ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
224 if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
225 throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
226 if(!is_exactly_a<numeric>(error))
227 throw std::runtime_error("For numerical integration the error should be a number.");
228
229 // Use lookup table to be potentially much faster.
230 static lookup_map lookup;
231 static symbol ivar("ivar");
232 ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
233 auto emi = lookup.find(error_and_integral(error, lookupex));
234 if (emi!=lookup.end())
235 return emi->second;
236
237 ex app = 0;
238 int i = 1;
246 vector<int> lvec(integral::max_integration_level+1);
247
248 avec[i] = a;
249 hvec[i] = (b-a)/2;
250 favec[i] = subsvalue(x, a, f);
251 fcvec[i] = subsvalue(x, a+hvec[i], f);
252 fbvec[i] = subsvalue(x, b, f);
253 svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
254 lvec[i] = 1;
255 errorvec[i] = error*abs(svec[i]);
256
257 while (i>0) {
258 ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
259 ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f);
260 ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6;
261 ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6;
262 ex nu1 = avec[i];
263 ex nu2 = favec[i];
264 ex nu3 = fcvec[i];
265 ex nu4 = fbvec[i];
266 ex nu5 = hvec[i];
267 // hopefully prevents a crash if the function is zero sometimes.
268 ex nu6 = max(errorvec[i], abs(s1+s2)*error);
269 ex nu7 = svec[i];
270 int nu8 = lvec[i];
271 --i;
272 if (abs(ex_to<numeric>(s1+s2-nu7)) <= nu6)
273 app+=(s1+s2);
274 else {
276 throw runtime_error("max integration level reached");
277 ++i;
278 avec[i] = nu1+nu5;
279 favec[i] = nu3;
280 fcvec[i] = fe;
281 fbvec[i] = nu4;
282 hvec[i] = nu5/2;
283 errorvec[i]=nu6/2;
284 svec[i] = s2;
285 lvec[i] = nu8+1;
286 ++i;
287 avec[i] = nu1;
288 favec[i] = nu2;
289 fcvec[i] = fd;
290 fbvec[i] = nu3;
291 hvec[i] = hvec[i-1];
292 errorvec[i]=errorvec[i-1];
293 svec[i] = s1;
294 lvec[i] = lvec[i-1];
295 }
296 }
297
298 lookup[error_and_integral(error, lookupex)]=app;
299 return app;
300}
301
302int integral::degree(const ex & s) const
303{
304 return ((b-a)*f).degree(s);
305}
306
307int integral::ldegree(const ex & s) const
308{
309 return ((b-a)*f).ldegree(s);
310}
311
313{
314 return f.eval_ncmul(v);
315}
316
317size_t integral::nops() const
318{
319 return 4;
320}
321
322ex integral::op(size_t i) const
323{
324 GINAC_ASSERT(i<4);
325
326 switch (i) {
327 case 0:
328 return x;
329 case 1:
330 return a;
331 case 2:
332 return b;
333 case 3:
334 return f;
335 default:
336 throw (std::out_of_range("integral::op() out of range"));
337 }
338}
339
341{
343 switch (i) {
344 case 0:
345 return x;
346 case 1:
347 return a;
348 case 2:
349 return b;
350 case 3:
351 return f;
352 default:
353 throw (std::out_of_range("integral::let_op() out of range"));
354 }
355}
356
358{
360 return *this;
361
362 ex newa = a.expand(options);
363 ex newb = b.expand(options);
364 ex newf = f.expand(options);
365
366 if (is_a<add>(newf)) {
367 exvector v;
368 v.reserve(newf.nops());
369 for (size_t i=0; i<newf.nops(); ++i)
370 v.push_back(integral(x, newa, newb, newf.op(i)).expand(options));
371 return ex(add(v)).expand(options);
372 }
373
374 if (is_a<mul>(newf)) {
375 ex prefactor = 1;
376 ex rest = 1;
377 for (size_t i=0; i<newf.nops(); ++i)
378 if (newf.op(i).has(x))
379 rest *= newf.op(i);
380 else
381 prefactor *= newf.op(i);
382 if (prefactor != 1)
383 return (prefactor*integral(x, newa, newb, rest)).expand(options);
384 }
385
386 if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb) &&
387 are_ex_trivially_equal(f, newf)) {
388 if (options==0)
389 this->setflag(status_flags::expanded);
390 return *this;
391 }
392
393 const integral & newint = dynallocate<integral>(x, newa, newb, newf);
394 if (options == 0)
396 return newint;
397}
398
400{
401 if (s==x)
402 throw(logic_error("differentiation with respect to dummy variable"));
403 return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s));
404}
405
406unsigned integral::return_type() const
407{
408 return f.return_type();
409}
410
412{
413 return f.return_type_tinfo();
414}
415
417{
418 ex conja = a.conjugate();
419 ex conjb = b.conjugate();
420 ex conjf = f.conjugate().subs(x.conjugate()==x);
421
422 if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb) &&
424 return *this;
425
426 return dynallocate<integral>(x, conja, conjb, conjf);
427}
428
430{
432 return this->expand().eval_integ();
433
434 if (f==x)
435 return b*b/2-a*a/2;
436 if (is_a<power>(f) && f.op(0)==x) {
437 if (f.op(1)==-1)
438 return log(b/a);
439 if (!f.op(1).has(x)) {
440 ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
441 return primit.subs(x==b)-primit.subs(x==a);
442 }
443 }
444
445 return *this;
446}
447
449} // namespace GiNaC
Interface to GiNaC's sums of expressions.
Archiving of GiNaC expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
Sum of expressions.
Definition: add.h:32
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:49
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
void ensure_if_modifiable() const
Ensure the object may be modified without hurting others, throws if this is not the case.
Definition: basic.cpp:894
friend class ex
Definition: basic.h:108
unsigned flags
of type status_flags
Definition: basic.h:302
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex eval_integ() const
Definition: ex.h:124
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:88
ex expand(unsigned options=0) const
Expand an expression.
Definition: ex.cpp:75
bool has(const ex &pattern, unsigned options=0) const
Definition: ex.h:151
ex evalf() const
Definition: ex.h:121
ex eval_ncmul(const exvector &v) const
Definition: ex.h:123
ex conjugate() const
Definition: ex.h:146
unsigned return_type() const
Definition: ex.h:230
return_type_t return_type_tinfo() const
Definition: ex.h:231
size_t nops() const
Definition: ex.h:135
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:841
int compare(const ex &other) const
Definition: ex.h:322
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
Definition: ex.cpp:56
ex op(size_t i) const
Definition: ex.h:136
Symbolic integral.
Definition: integral.h:34
ex evalf() const override
Evaluate object numerically.
Definition: integral.cpp:158
void do_print(const print_context &c, unsigned level) const
Definition: integral.cpp:95
static ex relative_integration_error
Definition: integral.h:75
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: integral.cpp:357
unsigned return_type() const override
Definition: integral.cpp:406
int degree(const ex &s) const override
Return degree of highest power in object s.
Definition: integral.cpp:302
ex eval() const override
Perform automatic non-interruptive term rewriting rules.
Definition: integral.cpp:144
size_t nops() const override
Number of operands/members.
Definition: integral.cpp:317
ex conjugate() const override
Definition: integral.cpp:416
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
Definition: integral.h:43
ex eval_integ() const override
Evaluate integrals, if result is known.
Definition: integral.cpp:429
ex eval_ncmul(const exvector &v) const override
Definition: integral.cpp:312
int ldegree(const ex &s) const override
Return degree of lowest power in object s.
Definition: integral.cpp:307
void do_print_latex(const print_latex &c, unsigned level) const
Definition: integral.cpp:108
integral(const ex &x_, const ex &a_, const ex &b_, const ex &f_)
Definition: integral.cpp:61
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
Definition: integral.cpp:340
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: integral.cpp:73
ex op(size_t i) const override
Return operand/member at position i.
Definition: integral.cpp:322
void archive(archive_node &n) const override
Save (a.k.a.
Definition: integral.cpp:82
return_type_t return_type_tinfo() const override
Definition: integral.cpp:411
static int max_integration_level
Definition: integral.h:74
ex derivative(const symbol &s) const override
Default implementation of ex::diff().
Definition: integral.cpp:399
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
Base class for print_contexts.
Definition: print.h:103
Context for latex-parsable output.
Definition: print.h:123
@ expanded
.expand(0) has already done its job (other expand() options ignore this flag)
Definition: flags.h:204
@ evaluated
.eval() has already done its job
Definition: flags.h:203
Basic CAS symbol.
Definition: symbol.h:39
static const bool value
Definition: factor.cpp:199
unsigned options
Definition: factor.cpp:2475
size_t n
Definition: factor.cpp:1432
size_t c
Definition: factor.cpp:757
ex x
Definition: factor.cpp:1610
Interface to GiNaC's initially known functions.
Interface to GiNaC's symbolic integral.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
bool haswild(const ex &x)
Check whether x has a wildcard anywhere as a subexpression.
Definition: wildcard.cpp:124
ex adaptivesimpson(const ex &x, const ex &a_in, const ex &b_in, const ex &f, const ex &error)
Numeric integration routine based upon the "Adaptive Quadrature" one in "Numerical Analysis" by Burde...
Definition: integral.cpp:219
B & dynallocate(Args &&... args)
Constructs a new (class basic or derived) B object on the heap.
Definition: basic.h:334
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2320
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
Definition: ex.h:699
print_func< print_dflt >(&diracone::do_print). print_func< print_latex >(&diracone
Definition: clifford.cpp:51
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(add, expairseq, print_func< print_context >(&add::do_print). print_func< print_latex >(&add::do_print_latex). print_func< print_csrc >(&add::do_print_csrc). print_func< print_tree >(&add::do_print_tree). print_func< print_python_repr >(&add::do_print_python_repr)) add
Definition: add.cpp:40
ex subsvalue(const ex &var, const ex &value, const ex &fun)
Definition: integral.cpp:181
const numeric log(const numeric &x)
Natural logarithm.
Definition: numeric.cpp:1450
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT_T(lst, basic, print_func< print_context >(&lst::do_print). print_func< print_tree >(&lst::do_print_tree)) template<> bool lst GINAC_BIND_UNARCHIVER(lst)
Specialization of container::info() for lst.
Definition: lst.cpp:42
const ex _ex0
Definition: utils.cpp:369
std::vector< ex > exvector
Definition: basic.h:48
map< error_and_integral, ex, error_and_integral_is_less > lookup_map
Definition: integral.cpp:210
Definition: ex.h:987
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
GiNaC's class registrar (for class basic and all classes derived from it).
Interface to relations between expressions.
bool operator()(const error_and_integral &e1, const error_and_integral &e2) const
Definition: integral.cpp:199
error_and_integral(const ex &err, const ex &integ)
Definition: integral.cpp:191
To distinguish between different kinds of non-commutative objects.
Definition: registrar.h:44
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
Interface to GiNaC's wildcard objects.

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.