GiNaC  1.8.0
integral.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2020 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 
37 using namespace std;
38 
39 namespace GiNaC {
40 
42  print_func<print_dflt>(&integral::do_print).
43  print_func<print_python>(&integral::do_print).
44  print_func<print_latex>(&integral::do_print_latex))
45 
46 
47 // default constructor
50 
51 integral::integral()
52  : x(dynallocate<symbol>())
53 {}
54 
56 // other constructors
58 
59 // public
60 
61 integral::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 
73 void integral::read_archive(const archive_node& n, lst& sym_lst)
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 
95 void 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 
108 void 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 
127 int 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 
181 ex 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 
210 typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
211 
219 ex 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 
302 int integral::degree(const ex & s) const
303 {
304  return ((b-a)*f).degree(s);
305 }
306 
307 int 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 
317 size_t integral::nops() const
318 {
319  return 4;
320 }
321 
322 ex 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 
340 ex & integral::let_op(size_t i)
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 
357 ex integral::expand(unsigned options) const
358 {
359  if (options==0 && (flags & status_flags::expanded))
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 
406 unsigned 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) &&
423  are_ex_trivially_equal(f, conjf))
424  return *this;
425 
426  return dynallocate<integral>(x, conja, conjb, conjf);
427 }
428 
430 {
431  if (!(flags & status_flags::expanded))
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&#39;s symbolic exponentiation (basis^exponent).
int ldegree(const ex &s) const override
Return degree of lowest power in object s.
Definition: integral.cpp:307
Interface to GiNaC&#39;s symbolic objects.
static int max_integration_level
Definition: integral.h:74
print_func< print_dflt >(&diracone::do_print). print_func< print_latex >(&diracone
Definition: clifford.cpp:51
int compare(const ex &other) const
Definition: ex.h:322
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
void do_print_latex(const print_latex &c, unsigned level) const
Definition: integral.cpp:108
ex conjugate() const override
Definition: integral.cpp:416
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
ex conjugate() const
Definition: ex.h:146
friend class ex
Definition: basic.h:108
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
ex expand(unsigned options=0) const
Definition: ex.cpp:73
integral(const ex &x_, const ex &a_, const ex &b_, const ex &f_)
Definition: integral.cpp:61
Definition: add.cpp:38
bool has(const ex &pattern, unsigned options=0) const
Definition: ex.h:151
size_t nops() const override
Number of operands/members.
Definition: integral.cpp:317
ex eval_ncmul(const exvector &v) const
Definition: ex.h:123
ex evalf() const
Definition: ex.h:121
Definition: ex.h:972
Archiving of GiNaC expressions.
Interface to GiNaC&#39;s sums of expressions.
This class is the ABC (abstract base class) of GiNaC&#39;s class hierarchy.
Definition: basic.h:104
size_t nops() const
Definition: ex.h:135
const ex _ex0
Definition: utils.cpp:177
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:684
int degree(const ex &s) const override
Return degree of highest power in object s.
Definition: integral.cpp:302
ex x
Definition: factor.cpp:1641
Interface to GiNaC&#39;s products of expressions.
Interface to GiNaC&#39;s wildcard objects.
ex eval_ncmul(const exvector &v) const override
Definition: integral.cpp:312
ex derivative(const symbol &s) const override
Default implementation of ex::diff().
Definition: integral.cpp:399
ex op(size_t i) const
Definition: ex.h:136
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
ex subsvalue(const ex &var, const ex &value, const ex &fun)
Definition: integral.cpp:181
ex eval_integ() const override
Evaluate integrals, if result is known.
Definition: integral.cpp:429
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
Definition: integral.h:43
static const bool value
Definition: factor.cpp:231
bool haswild(const ex &x)
Check whether x has a wildcard anywhere as a subexpression.
Definition: wildcard.cpp:124
Context for latex-parsable output.
Definition: print.h:122
void do_print(const print_context &c, unsigned level) const
Definition: integral.cpp:95
ex eval_integ() const
Definition: ex.h:124
static ex relative_integration_error
Definition: integral.h:75
GiNaC&#39;s class registrar (for class basic and all classes derived from it).
unsigned options
Definition: factor.cpp:2480
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
B & dynallocate(Args &&... args)
Constructs a new (class basic or derived) B object on the heap.
Definition: basic.h:334
Sum of expressions.
Definition: add.h:31
const numeric log(const numeric &x)
Natural logarithm.
Definition: numeric.cpp:1450
To distinguish between different kinds of non-commutative objects.
Definition: registrar.h:43
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:86
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
std::vector< ex > exvector
Definition: basic.h:46
Interface to GiNaC&#39;s overloaded operators.
#define GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(classname, supername, options)
Macro for inclusion in the implementation of each registered class.
Definition: registrar.h:185
return_type_t return_type_tinfo() const override
Definition: integral.cpp:411
return_type_t return_type_tinfo() const
Definition: ex.h:231
size_t n
Definition: factor.cpp:1463
Interface to GiNaC&#39;s symbolic integral.
Base class for print_contexts.
Definition: print.h:102
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
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: integral.cpp:73
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
error_and_integral(const ex &err, const ex &integ)
Definition: integral.cpp:191
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
Definition: ex.cpp:56
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:48
ex op(size_t i) const override
Return operand/member at position i.
Definition: integral.cpp:322
Lightweight wrapper for GiNaC&#39;s symbolic objects.
Definition: ex.h:72
ex evalf() const override
Evaluate object numerically.
Definition: integral.cpp:158
bool operator()(const error_and_integral &e1, const error_and_integral &e2) const
Definition: integral.cpp:199
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:38
unsigned return_type() const override
Definition: integral.cpp:406
Interface to GiNaC&#39;s initially known functions.
map< error_and_integral, ex, error_and_integral_is_less > lookup_map
Definition: integral.cpp:210
Basic CAS symbol.
Definition: symbol.h:38
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
unsigned return_type() const
Definition: ex.h:230
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
.expand(0) has already done its job (other expand() options ignore this flag)
Definition: flags.h:204
Interface to relations between expressions.
Makes the interface to the underlying bignum package available.
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
Definition: integral.cpp:340
size_t c
Definition: factor.cpp:770
void archive(archive_node &n) const override
Save (a.k.a.
Definition: integral.cpp:82
.eval() has already done its job
Definition: flags.h:203
unsigned flags
of type status_flags
Definition: basic.h:302
ex eval() const override
Perform automatic non-interruptive term rewriting rules.
Definition: integral.cpp:144
GINAC_BIND_UNARCHIVER(add)
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: integral.cpp:357
Symbolic integral.
Definition: integral.h:33

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