ginac/integral.cpp

   1 /** @file integral.cpp
   2  *
   3  *  Implementation of GiNaC's symbolic  integral. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2016 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
  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))
  44
  45
  46 //////////
  47 // default constructor
  48 //////////
  49
  50 integral::integral()
  51                 : x(dynallocate<symbol>())
  52 {}
  53
  54 //////////
  55 // other constructors
  56 //////////
  57
  58 // public
  59
  60 integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
  61                 :  x(x_), a(a_), b(b_), f(f_)
  62 {
  63         if (!is_a<symbol>(x)) {
  64                 throw(std::invalid_argument("first argument of integral must be of type symbol"));
  65         }
  66 }
  67
  68 //////////
  69 // archiving
  70 //////////
  71
  72 void integral::read_archive(const archive_node& n, lst& sym_lst)
  73 {
  74         inherited::read_archive(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);
  79 }
  80
  81 void integral::archive(archive_node & n) const
  82 {
  83         inherited::archive(n);
  84         n.add_ex("x", x);
  85         n.add_ex("a", a);
  86         n.add_ex("b", b);
  87         n.add_ex("f", f);
  88 }
  89
  90 //////////
  91 // functions overriding virtual functions from base classes
  92 //////////
  93
  94 void integral::do_print(const print_context & c, unsigned level) const
  95 {
  96         c.s << "integral(";
  97         x.print(c);
  98         c.s << ",";
  99         a.print(c);
 100         c.s << ",";
 101         b.print(c);
 102         c.s << ",";
 103         f.print(c);
 104         c.s << ")";
 105 }
 106
 107 void integral::do_print_latex(const print_latex & c, unsigned level) const
 108 {
 109         string varname = ex_to<symbol>(x).get_name();
 110         if (level > precedence())
 111                 c.s << "\\left(";
 112         c.s << "\\int_{";
 113         a.print(c);
 114         c.s << "}^{";
 115         b.print(c);
 116         c.s << "} d";
 117         if (varname.size() > 1)
 118                 c.s << "\\," << varname << "\\:";
 119         else
 120                 c.s << varname << "\\,";
 121         f.print(c,precedence());
 122         if (level > precedence())
 123                 c.s << "\\right)";
 124 }
 125
 126 int integral::compare_same_type(const basic & other) const
 127 {
 128         GINAC_ASSERT(is_exactly_a<integral>(other));
 129         const integral &o = static_cast<const integral &>(other);
 130
 131         int cmpval = x.compare(o.x);
 132         if (cmpval)
 133                 return cmpval;
 134         cmpval = a.compare(o.a);
 135         if (cmpval)
 136                 return cmpval;
 137         cmpval = b.compare(o.b);
 138         if (cmpval)
 139                 return cmpval;
 140         return f.compare(o.f);
 141 }
 142
 143 ex integral::eval() const
 144 {
 145         if (flags & status_flags::evaluated)
 146                 return *this;
 147
 148         if (!f.has(x) && !haswild(f))
 149                 return b*f-a*f;
 150
 151         if (a==b)
 152                 return _ex0;
 153
 154         return this->hold();
 155 }
 156
 157 ex integral::evalf() const
 158 {
 159         ex ea = a.evalf();
 160         ex eb = b.evalf();
 161         ex ef = f.evalf();
 162
 163         // 12.34 is just an arbitrary number used to check whether a number
 164         // results after substituting a number for the integration variable.
 165         if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb) &&
 166             is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
 167                         return adaptivesimpson(x, ea, eb, ef);
 168         }
 169
 170         if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb) &&
 171             are_ex_trivially_equal(f, ef))
 172                         return *this;
 173                 else
 174                         return dynallocate<integral>(x, ea, eb, ef);
 175 }
 176
 177 int integral::max_integration_level = 15;
 178 ex integral::relative_integration_error = 1e-8;
 179
 180 ex subsvalue(const ex & var, const ex & value, const ex & fun)
 181 {
 182         ex result = fun.subs(var==value).evalf();
 183         if (is_a<numeric>(result))
 184                 return result;
 185         throw logic_error("integrand does not evaluate to numeric");
 186 }
 187
 188 struct error_and_integral
 189 {
 190         error_and_integral(const ex &err, const ex &integ)
 191                 :error(err), integral(integ){}
 192         ex error;
 193         ex integral;
 194 };
 195
 196 struct error_and_integral_is_less
 197 {
 198         bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
 199         {
 200                 int c = e1.integral.compare(e2.integral);
 201                 if(c < 0)
 202                         return true;
 203                 if(c > 0)
 204                         return false;
 205                 return ex_is_less()(e1.error, e2.error);
 206         }
 207 };
 208
 209 typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
 210
 211 /** Numeric integration routine based upon the "Adaptive Quadrature" one
 212   * in "Numerical Analysis" by Burden and Faires. Parameters are integration
 213   * variable, left boundary, right boundary, function to be integrated and
 214   * the relative integration error. The function should evalf into a number
 215   * after substituting the integration variable by a number. Another thing
 216   * to note is that this implementation is no good at integrating functions
 217   * with discontinuities. */
 218 ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
 219 {
 220         // Check whether boundaries and error are numbers.
 221         ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
 222         ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
 223         if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
 224                 throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
 225         if(!is_exactly_a<numeric>(error))
 226                 throw std::runtime_error("For numerical integration the error should be a number.");
 227
 228         // Use lookup table to be potentially much faster.
 229         static lookup_map lookup;
 230         static symbol ivar("ivar");
 231         ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
 232         lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
 233         if (emi!=lookup.end())
 234                 return emi->second;
 235
 236         ex app = 0;
 237         int i = 1;
 238         exvector avec(integral::max_integration_level+1);
 239         exvector hvec(integral::max_integration_level+1);
 240         exvector favec(integral::max_integration_level+1);
 241         exvector fbvec(integral::max_integration_level+1);
 242         exvector fcvec(integral::max_integration_level+1);
 243         exvector svec(integral::max_integration_level+1);
 244         exvector errorvec(integral::max_integration_level+1);
 245         vector<int> lvec(integral::max_integration_level+1);
 246
 247         avec[i] = a;
 248         hvec[i] = (b-a)/2;
 249         favec[i] = subsvalue(x, a, f);
 250         fcvec[i] = subsvalue(x, a+hvec[i], f);
 251         fbvec[i] = subsvalue(x, b, f);
 252         svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
 253         lvec[i] = 1;
 254         errorvec[i] = error*abs(svec[i]);
 255
 256         while (i>0) {
 257                 ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
 258                 ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f);
 259                 ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6;
 260                 ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6;
 261                 ex nu1 = avec[i];
 262                 ex nu2 = favec[i];
 263                 ex nu3 = fcvec[i];
 264                 ex nu4 = fbvec[i];
 265                 ex nu5 = hvec[i];
 266                 // hopefully prevents a crash if the function is zero sometimes.
 267                 ex nu6 = max(errorvec[i], abs(s1+s2)*error);
 268                 ex nu7 = svec[i];
 269                 int nu8 = lvec[i];
 270                 --i;
 271                 if (abs(ex_to<numeric>(s1+s2-nu7)) <= nu6)
 272                         app+=(s1+s2);
 273                 else {
 274                         if (nu8>=integral::max_integration_level)
 275                                 throw runtime_error("max integration level reached");
 276                         ++i;
 277                         avec[i] = nu1+nu5;
 278                         favec[i] = nu3;
 279                         fcvec[i] = fe;
 280                         fbvec[i] = nu4;
 281                         hvec[i] = nu5/2;
 282                         errorvec[i]=nu6/2;
 283                         svec[i] = s2;
 284                         lvec[i] = nu8+1;
 285                         ++i;
 286                         avec[i] = nu1;
 287                         favec[i] = nu2;
 288                         fcvec[i] = fd;
 289                         fbvec[i] = nu3;
 290                         hvec[i] = hvec[i-1];
 291                         errorvec[i]=errorvec[i-1];
 292                         svec[i] = s1;
 293                         lvec[i] = lvec[i-1];
 294                 }
 295         }
 296
 297         lookup[error_and_integral(error, lookupex)]=app;
 298         return app;
 299 }
 300
 301 int integral::degree(const ex & s) const
 302 {
 303         return ((b-a)*f).degree(s);
 304 }
 305
 306 int integral::ldegree(const ex & s) const
 307 {
 308         return ((b-a)*f).ldegree(s);
 309 }
 310
 311 ex integral::eval_ncmul(const exvector & v) const
 312 {
 313         return f.eval_ncmul(v);
 314 }
 315
 316 size_t integral::nops() const
 317 {
 318         return 4;
 319 }
 320
 321 ex integral::op(size_t i) const
 322 {
 323         GINAC_ASSERT(i<4);
 324
 325         switch (i) {
 326                 case 0:
 327                         return x;
 328                 case 1:
 329                         return a;
 330                 case 2:
 331                         return b;
 332                 case 3:
 333                         return f;
 334                 default:
 335                         throw (std::out_of_range("integral::op() out of range"));
 336         }
 337 }
 338
 339 ex & integral::let_op(size_t i)
 340 {
 341         ensure_if_modifiable();
 342         switch (i) {
 343                 case 0:
 344                         return x;
 345                 case 1:
 346                         return a;
 347                 case 2:
 348                         return b;
 349                 case 3:
 350                         return f;
 351                 default:
 352                         throw (std::out_of_range("integral::let_op() out of range"));
 353         }
 354 }
 355
 356 ex integral::expand(unsigned options) const
 357 {
 358         if (options==0 && (flags & status_flags::expanded))
 359                 return *this;
 360
 361         ex newa = a.expand(options);
 362         ex newb = b.expand(options);
 363         ex newf = f.expand(options);
 364
 365         if (is_a<add>(newf)) {
 366                 exvector v;
 367                 v.reserve(newf.nops());
 368                 for (size_t i=0; i<newf.nops(); ++i)
 369                         v.push_back(integral(x, newa, newb, newf.op(i)).expand(options));
 370                 return ex(add(v)).expand(options);
 371         }
 372
 373         if (is_a<mul>(newf)) {
 374                 ex prefactor = 1;
 375                 ex rest = 1;
 376                 for (size_t i=0; i<newf.nops(); ++i)
 377                         if (newf.op(i).has(x))
 378                                 rest *= newf.op(i);
 379                         else
 380                                 prefactor *= newf.op(i);
 381                 if (prefactor != 1)
 382                         return (prefactor*integral(x, newa, newb, rest)).expand(options);
 383         }
 384
 385         if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb) &&
 386             are_ex_trivially_equal(f, newf)) {
 387                 if (options==0)
 388                         this->setflag(status_flags::expanded);
 389                 return *this;
 390         }
 391
 392         const integral & newint = dynallocate<integral>(x, newa, newb, newf);
 393         if (options == 0)
 394                 newint.setflag(status_flags::expanded);
 395         return newint;
 396 }
 397
 398 ex integral::derivative(const symbol & s) const
 399 {
 400         if (s==x)
 401                 throw(logic_error("differentiation with respect to dummy variable"));
 402         return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s));
 403 }
 404
 405 unsigned integral::return_type() const
 406 {
 407         return f.return_type();
 408 }
 409
 410 return_type_t integral::return_type_tinfo() const
 411 {
 412         return f.return_type_tinfo();
 413 }
 414
 415 ex integral::conjugate() const
 416 {
 417         ex conja = a.conjugate();
 418         ex conjb = b.conjugate();
 419         ex conjf = f.conjugate().subs(x.conjugate()==x);
 420
 421         if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb) &&
 422             are_ex_trivially_equal(f, conjf))
 423                 return *this;
 424
 425         return dynallocate<integral>(x, conja, conjb, conjf);
 426 }
 427
 428 ex integral::eval_integ() const
 429 {
 430         if (!(flags & status_flags::expanded))
 431                 return this->expand().eval_integ();
 432
 433         if (f==x)
 434                 return b*b/2-a*a/2;
 435         if (is_a<power>(f) && f.op(0)==x) {
 436                 if (f.op(1)==-1)
 437                         return log(b/a);
 438                 if (!f.op(1).has(x)) {
 439                         ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
 440                         return primit.subs(x==b)-primit.subs(x==a);
 441                 }
 442         }
 443
 444         return *this;
 445 }
 446
 447 GINAC_BIND_UNARCHIVER(integral);
 448 } // namespace GiNaC