ginac/integral.cpp

   1 /** @file integral.cpp
   2  *
   3  *  Implementation of GiNaC's symbolic  integral. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2011 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                 :
  52                 x((new symbol())->setflag(status_flags::dynallocated))
  53 {}
  54
  55 //////////
  56 // other constructors
  57 //////////
  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
  69 //////////
  70 // archiving
  71 //////////
  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
  82 void integral::archive(archive_node & n) const
  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
  91 //////////
  92 // functions overriding virtual functions from base classes
  93 //////////
  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
 144 ex integral::eval(int level) const
 145 {
 146         if ((level==1) && (flags & status_flags::evaluated))
 147                 return *this;
 148         if (level == -max_recursion_level)
 149                 throw(std::runtime_error("max recursion level reached"));
 150
 151         ex eintvar = (level==1) ? x : x.eval(level-1);
 152         ex ea      = (level==1) ? a : a.eval(level-1);
 153         ex eb      = (level==1) ? b : b.eval(level-1);
 154         ex ef      = (level==1) ? f : f.eval(level-1);
 155
 156         if (!ef.has(eintvar) && !haswild(ef))
 157                 return eb*ef-ea*ef;
 158
 159         if (ea==eb)
 160                 return _ex0;
 161
 162         if (are_ex_trivially_equal(eintvar,x) && are_ex_trivially_equal(ea,a)
 163                         && are_ex_trivially_equal(eb,b) && are_ex_trivially_equal(ef,f))
 164                 return this->hold();
 165         return (new integral(eintvar, ea, eb, ef))
 166                 ->setflag(status_flags::dynallocated | status_flags::evaluated);
 167 }
 168
 169 ex integral::evalf(int level) const
 170 {
 171         ex ea;
 172         ex eb;
 173         ex ef;
 174
 175         if (level==1) {
 176                 ea = a;
 177                 eb = b;
 178                 ef = f;
 179         } else if (level == -max_recursion_level) {
 180                 throw(runtime_error("max recursion level reached"));
 181         } else {
 182                 ea = a.evalf(level-1);
 183                 eb = b.evalf(level-1);
 184                 ef = f.evalf(level-1);
 185         }
 186
 187         // 12.34 is just an arbitrary number used to check whether a number
 188         // results after subsituting a number for the integration variable.
 189         if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb)
 190                         && is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
 191                         return adaptivesimpson(x, ea, eb, ef);
 192         }
 193
 194         if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb)
 195                                 && are_ex_trivially_equal(f, ef))
 196                         return *this;
 197                 else
 198                         return (new integral(x, ea, eb, ef))
 199                                 ->setflag(status_flags::dynallocated);
 200 }
 201
 202 int integral::max_integration_level = 15;
 203 ex integral::relative_integration_error = 1e-8;
 204
 205 ex subsvalue(const ex & var, const ex & value, const ex & fun)
 206 {
 207         ex result = fun.subs(var==value).evalf();
 208         if (is_a<numeric>(result))
 209                 return result;
 210         throw logic_error("integrand does not evaluate to numeric");
 211 }
 212
 213 struct error_and_integral
 214 {
 215         error_and_integral(const ex &err, const ex &integ)
 216                 :error(err), integral(integ){}
 217         ex error;
 218         ex integral;
 219 };
 220
 221 struct error_and_integral_is_less
 222 {
 223         bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
 224         {
 225                 int c = e1.integral.compare(e2.integral);
 226                 if(c < 0)
 227                         return true;
 228                 if(c > 0)
 229                         return false;
 230                 return ex_is_less()(e1.error, e2.error);
 231         }
 232 };
 233
 234 typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
 235
 236 /** Numeric integration routine based upon the "Adaptive Quadrature" one
 237   * in "Numerical Analysis" by Burden and Faires. Parameters are integration
 238   * variable, left boundary, right boundary, function to be integrated and
 239   * the relative integration error. The function should evalf into a number
 240   * after substituting the integration variable by a number. Another thing
 241   * to note is that this implementation is no good at integrating functions
 242   * with discontinuities. */
 243 ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
 244 {
 245         // Check whether boundaries and error are numbers.
 246         ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
 247         ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
 248         if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
 249                 throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
 250         if(!is_exactly_a<numeric>(error))
 251                 throw std::runtime_error("For numerical integration the error should be a number.");
 252
 253         // Use lookup table to be potentially much faster.
 254         static lookup_map lookup;
 255         static symbol ivar("ivar");
 256         ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
 257         lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
 258         if (emi!=lookup.end())
 259                 return emi->second;
 260
 261         ex app = 0;
 262         int i = 1;
 263         exvector avec(integral::max_integration_level+1);
 264         exvector hvec(integral::max_integration_level+1);
 265         exvector favec(integral::max_integration_level+1);
 266         exvector fbvec(integral::max_integration_level+1);
 267         exvector fcvec(integral::max_integration_level+1);
 268         exvector svec(integral::max_integration_level+1);
 269         exvector errorvec(integral::max_integration_level+1);
 270         vector<int> lvec(integral::max_integration_level+1);
 271
 272         avec[i] = a;
 273         hvec[i] = (b-a)/2;
 274         favec[i] = subsvalue(x, a, f);
 275         fcvec[i] = subsvalue(x, a+hvec[i], f);
 276         fbvec[i] = subsvalue(x, b, f);
 277         svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
 278         lvec[i] = 1;
 279         errorvec[i] = error*abs(svec[i]);
 280
 281         while (i>0) {
 282                 ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
 283                 ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f);
 284                 ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6;
 285                 ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6;
 286                 ex nu1 = avec[i];
 287                 ex nu2 = favec[i];
 288                 ex nu3 = fcvec[i];
 289                 ex nu4 = fbvec[i];
 290                 ex nu5 = hvec[i];
 291                 // hopefully prevents a crash if the function is zero sometimes.
 292                 ex nu6 = max(errorvec[i], abs(s1+s2)*error);
 293                 ex nu7 = svec[i];
 294                 int nu8 = lvec[i];
 295                 --i;
 296                 if (abs(ex_to<numeric>(s1+s2-nu7)) <= nu6)
 297                         app+=(s1+s2);
 298                 else {
 299                         if (nu8>=integral::max_integration_level)
 300                                 throw runtime_error("max integration level reached");
 301                         ++i;
 302                         avec[i] = nu1+nu5;
 303                         favec[i] = nu3;
 304                         fcvec[i] = fe;
 305                         fbvec[i] = nu4;
 306                         hvec[i] = nu5/2;
 307                         errorvec[i]=nu6/2;
 308                         svec[i] = s2;
 309                         lvec[i] = nu8+1;
 310                         ++i;
 311                         avec[i] = nu1;
 312                         favec[i] = nu2;
 313                         fcvec[i] = fd;
 314                         fbvec[i] = nu3;
 315                         hvec[i] = hvec[i-1];
 316                         errorvec[i]=errorvec[i-1];
 317                         svec[i] = s1;
 318                         lvec[i] = lvec[i-1];
 319                 }
 320         }
 321
 322         lookup[error_and_integral(error, lookupex)]=app;
 323         return app;
 324 }
 325
 326 int integral::degree(const ex & s) const
 327 {
 328         return ((b-a)*f).degree(s);
 329 }
 330
 331 int integral::ldegree(const ex & s) const
 332 {
 333         return ((b-a)*f).ldegree(s);
 334 }
 335
 336 ex integral::eval_ncmul(const exvector & v) const
 337 {
 338         return f.eval_ncmul(v);
 339 }
 340
 341 size_t integral::nops() const
 342 {
 343         return 4;
 344 }
 345
 346 ex integral::op(size_t i) const
 347 {
 348         GINAC_ASSERT(i<4);
 349
 350         switch (i) {
 351                 case 0:
 352                         return x;
 353                 case 1:
 354                         return a;
 355                 case 2:
 356                         return b;
 357                 case 3:
 358                         return f;
 359                 default:
 360                         throw (std::out_of_range("integral::op() out of range"));
 361         }
 362 }
 363
 364 ex & integral::let_op(size_t i)
 365 {
 366         ensure_if_modifiable();
 367         switch (i) {
 368                 case 0:
 369                         return x;
 370                 case 1:
 371                         return a;
 372                 case 2:
 373                         return b;
 374                 case 3:
 375                         return f;
 376                 default:
 377                         throw (std::out_of_range("integral::let_op() out of range"));
 378         }
 379 }
 380
 381 ex integral::expand(unsigned options) const
 382 {
 383         if (options==0 && (flags & status_flags::expanded))
 384                 return *this;
 385
 386         ex newa = a.expand(options);
 387         ex newb = b.expand(options);
 388         ex newf = f.expand(options);
 389
 390         if (is_a<add>(newf)) {
 391                 exvector v;
 392                 v.reserve(newf.nops());
 393                 for (size_t i=0; i<newf.nops(); ++i)
 394                         v.push_back(integral(x, newa, newb, newf.op(i)).expand(options));
 395                 return ex(add(v)).expand(options);
 396         }
 397
 398         if (is_a<mul>(newf)) {
 399                 ex prefactor = 1;
 400                 ex rest = 1;
 401                 for (size_t i=0; i<newf.nops(); ++i)
 402                         if (newf.op(i).has(x))
 403                                 rest *= newf.op(i);
 404                         else
 405                                 prefactor *= newf.op(i);
 406                 if (prefactor != 1)
 407                         return (prefactor*integral(x, newa, newb, rest)).expand(options);
 408         }
 409
 410         if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb)
 411                         && are_ex_trivially_equal(f, newf)) {
 412                 if (options==0)
 413                         this->setflag(status_flags::expanded);
 414                 return *this;
 415         }
 416
 417         const basic & newint = (new integral(x, newa, newb, newf))
 418                 ->setflag(status_flags::dynallocated);
 419         if (options == 0)
 420                 newint.setflag(status_flags::expanded);
 421         return newint;
 422 }
 423
 424 ex integral::derivative(const symbol & s) const
 425 {
 426         if (s==x)
 427                 throw(logic_error("differentiation with respect to dummy variable"));
 428         return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s));
 429 }
 430
 431 unsigned integral::return_type() const
 432 {
 433         return f.return_type();
 434 }
 435
 436 return_type_t integral::return_type_tinfo() const
 437 {
 438         return f.return_type_tinfo();
 439 }
 440
 441 ex integral::conjugate() const
 442 {
 443         ex conja = a.conjugate();
 444         ex conjb = b.conjugate();
 445         ex conjf = f.conjugate().subs(x.conjugate()==x);
 446
 447         if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb)
 448                         && are_ex_trivially_equal(f, conjf))
 449                 return *this;
 450
 451         return (new integral(x, conja, conjb, conjf))
 452                 ->setflag(status_flags::dynallocated);
 453 }
 454
 455 ex integral::eval_integ() const
 456 {
 457         if (!(flags & status_flags::expanded))
 458                 return this->expand().eval_integ();
 459
 460         if (f==x)
 461                 return b*b/2-a*a/2;
 462         if (is_a<power>(f) && f.op(0)==x) {
 463                 if (f.op(1)==-1)
 464                         return log(b/a);
 465                 if (!f.op(1).has(x)) {
 466                         ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
 467                         return primit.subs(x==b)-primit.subs(x==a);
 468                 }
 469         }
 470
 471         return *this;
 472 }
 473
 474 GINAC_BIND_UNARCHIVER(integral);
 475 } // namespace GiNaC