sympy/integrals/integrals.py

   1
   2 from sympy.core import Basic, S, C, Symbol, Wild, Pow, sympify
   3 from sympy.core.methods import RelMeths, ArithMeths
   4
   5 from sympy.integrals.trigonometry import trigintegrate
   6 from sympy.integrals.risch import heurisch
   7 from sympy.utilities import threaded
   8 from sympy.simplify import apart
   9 from sympy.series import limit
  10 from sympy.polys import Poly
  11 from sympy.solvers import solve
  12
  13 class Integral(Basic, RelMeths, ArithMeths):
  14     """Represents unevaluated integral."""
  15
  16     def __new__(cls, function, *symbols, **assumptions):
  17         function = sympify(function)
  18
  19         if function.is_Number:
  20             if function is S.NaN:
  21                 return S.NaN
  22             elif function is S.Infinity:
  23                 return S.Infinity
  24             elif function is S.NegativeInfinity:
  25                 return S.NegativeInfinity
  26
  27         if symbols:
  28             limits = []
  29
  30             for V in symbols:
  31                 if isinstance(V, Symbol):
  32                     limits.append((V,None))
  33                     continue
  34                 elif isinstance(V, (tuple, list)):
  35                     if len(V) == 3:
  36                         limits.append( (V[0],tuple(V[1:])) )
  37                         continue
  38                     elif len(V) == 1:
  39                         if isinstance(V[0], Symbol):
  40                             limits.append((V[0],None))
  41                             continue
  42
  43                 raise ValueError("Invalid integration variable or limits")
  44         else:
  45             # no symbols provided -- let's compute full antiderivative
  46             limits = [(symb,None) for symb in function.atoms(Symbol)]
  47
  48             if not limits:
  49                 return function
  50
  51         obj = Basic.__new__(cls, **assumptions)
  52         obj._args = (function, tuple(limits))
  53
  54         return obj
  55
  56     @property
  57     def function(self):
  58         return self._args[0]
  59
  60     @property
  61     def limits(self):
  62         return self._args[1]
  63
  64     @property
  65     def variables(self):
  66         variables = []
  67
  68         for x,ab in self.limits:
  69             variables.append(x)
  70
  71         return variables
  72
  73     def transform(self, x, mapping, inverse=False):
  74         """
  75         Replace the integration variable x in the integrand with the
  76         expression given by `mapping`, e.g. 2*x or 1/x. The integrand and
  77         endpoints are rescaled to preserve the value of the original
  78         integral.
  79
  80         In effect, this performs a variable substitution (although the
  81         symbol remains unchanged; follow up with subs to obtain a
  82         new symbol.)
  83
  84         With inverse=True, the inverse transformation is performed.
  85
  86         The mapping must be uniquely invertible (e.g. a linear or linear
  87         fractional transformation).
  88         """
  89         if x not in self.variables:
  90             return self
  91         limits = self.limits
  92         function = self.function
  93         y = Symbol('y', dummy=True)
  94         inverse_mapping = solve(mapping.subs(x,y)-x, y)
  95         if len(inverse_mapping) != 1 or not inverse_mapping[0].has(x):
  96             raise ValueError("The mapping must be uniquely invertible")
  97         inverse_mapping = inverse_mapping[0]
  98         if inverse:
  99             mapping, inverse_mapping = inverse_mapping, mapping
 100         function = function.subs(x, mapping) * mapping.diff(x)
 101         newlimits = []
 102         for sym, limit in limits:
 103             if sym == x and limit and len(limit) == 2:
 104                 a, b = limit
 105                 a = inverse_mapping.subs(x, a)
 106                 b = inverse_mapping.subs(x, b)
 107                 if a == b:
 108                     raise ValueError("The mapping must transform the "
 109                         "endpoints into separate points")
 110                 if a > b:
 111                     a, b = b, a
 112                     function = -function
 113                 newlimits.append((sym, a, b))
 114             else:
 115                 newlimits.append((sym, limit))
 116         return Integral(function, *newlimits)
 117
 118     def doit(self, **hints):
 119         if not hints.get('integrals', True):
 120             return self
 121
 122         function = self.function
 123
 124         for x,ab in self.limits:
 125             antideriv = self._eval_integral(function, x)
 126
 127             if antideriv is None:
 128                 return self
 129             else:
 130                 if ab is None:
 131                     function = antideriv
 132                 else:
 133                     a,b = ab
 134                     A = antideriv.subs(x, a)
 135
 136                     if A is S.NaN:
 137                         A = limit(antideriv, x, a)
 138                     if A is S.NaN:
 139                         return self
 140
 141                     B = antideriv.subs(x, b)
 142
 143                     if B is S.NaN:
 144                         B = limit(antideriv, x, b)
 145                     if B is S.NaN:
 146                         return self
 147
 148                     function = B - A
 149
 150         return function
 151
 152     def _eval_integral(self, f, x):
 153         """Calculate the antiderivative to the function f(x).
 154
 155         This is a powerful function that should in theory be able to integrate
 156         everything that can be integrated. If you find something, that it
 157         doesn't, it is easy to implement it.
 158
 159         (1) Simple heuristics (based on pattern matching and integral table):
 160
 161          - most frequently used functions (eg. polynomials)
 162          - functions non-integrable by any of the following algorithms (eg.
 163            exp(-x**2))
 164
 165         (2) Integration of rational functions:
 166
 167          (a) using apart() - apart() is full partial fraction decomposition
 168          procedure based on Bronstein-Salvy algorithm. It gives formal
 169          decomposition with no polynomial factorization at all (so it's fast
 170          and gives the most general results). However it needs much better
 171          implementation of RootsOf class (if fact any implementation).
 172          (b) using Trager's algorithm - possibly faster than (a) but needs
 173          implementation :)
 174
 175         (3) Whichever implementation of pmInt (Mateusz, Kirill's or a
 176         combination of both).
 177
 178           - this way we can handle efficiently huge class of elementary and
 179             special functions
 180
 181         (4) Recursive Risch algorithm as described in Bronstein's integration
 182         tutorial.
 183
 184           - this way we can handle those integrable functions for which (3)
 185             fails
 186
 187         (5) Powerful heuristics based mostly on user defined rules.
 188
 189          - handle complicated, rarely used cases
 190         """
 191
 192         # if it is a poly(x) then let the polynomial integrate itself (fast)
 193         #
 194         # It is important to make this check first, otherwise the other code
 195         # will return a sympy expression instead of a Polynomial.
 196         #
 197         # see Polynomial for details.
 198         if isinstance(f, Poly):
 199             return f.integrate(x)
 200
 201         # let's cut it short if `f` does not depend on `x`
 202         if not f.has(x):
 203             return f*x
 204
 205         # try to convert to poly(x) and then integrate if successful (fast)
 206         poly = f.as_poly(x)
 207
 208         if poly is not None:
 209             return poly.integrate(x).as_basic()
 210
 211         # since Integral(f=g1+g2+...) == Integral(g1) + Integral(g2) + ...
 212         # we are going to handle Add terms separately,
 213         # if `f` is not Add -- we only have one term
 214         if not f.is_Add:
 215             f = [f]
 216
 217         parts = []
 218
 219         if isinstance(f, Basic):
 220             f = f.args
 221         for g in f:
 222             coeff, g = g.as_independent(x)
 223
 224             # g(x) = const
 225             if g is S.One:
 226                 parts.append(coeff * x)
 227                 continue
 228
 229             #               c
 230             # g(x) = (a*x+b)
 231             if g.is_Pow and not g.exp.has(x):
 232                 a = Wild('a', exclude=[x])
 233                 b = Wild('b', exclude=[x])
 234
 235                 M = g.base.match(a*x + b)
 236
 237                 if M is not None:
 238                     if g.exp == -1:
 239                         h = C.log(g.base)
 240                     else:
 241                         h = g.base**(g.exp+1) / (g.exp+1)
 242
 243                     parts.append(coeff * h / M[a])
 244                     continue
 245
 246             #        poly(x)
 247             # g(x) = -------
 248             #        poly(x)
 249             if g.is_fraction(x):
 250                 h = self._eval_integral(apart(g, x), x)
 251                 parts.append(coeff * h)
 252                 continue
 253
 254             # g(x) = Mul(trig)
 255             h = trigintegrate(g, x)
 256             if h is not None:
 257                 parts.append(coeff * h)
 258                 continue
 259
 260             # fall back to the more general algorithm
 261             h = heurisch(g, x, hints=[])
 262
 263             if h is not None:
 264                 parts.append(coeff * h)
 265             else:
 266                 return None
 267
 268         return C.Add(*parts)
 269
 270 @threaded(use_add=False)
 271 def integrate(*args, **kwargs):
 272     """integrate(f, var, ...)
 273
 274        Compute definite or indefinite integral of one or more variables
 275        using Risch-Norman algorithm and table lookup. This procedure is
 276        able to handle elementary algebraic and transcendental functions
 277        and also a huge class of special functions, including Airy,
 278        Bessel, Whittaker and Lambert.
 279
 280        var can be:
 281
 282        - a symbol                   -- indefinite integration
 283        - a tuple (symbol, a, b)     -- definite integration
 284
 285        Several variables can be specified, in which case the result is multiple
 286        integration.
 287
 288        Also, if no var is specified at all, then full-antiderivative of f is
 289        returned. This is equivalent of integrating f over all it's variables.
 290
 291        Examples
 292        --------
 293
 294        >>> from sympy import *
 295        >>> x, y = symbols('xy')
 296
 297        >>> integrate(x*y, x)
 298        (1/2)*y*x**2
 299
 300        >>> integrate(log(x), x)
 301        -x + x*log(x)
 302
 303        >>> integrate(x)
 304        (1/2)*x**2
 305
 306        >>> integrate(x*y)
 307        (1/4)*x**2*y**2
 308
 309        See also the doctest of Integral._eval_integral(), which explains
 310        thoroughly the strategy that SymPy uses for integration.
 311
 312     """
 313     integral = Integral(*args, **kwargs)
 314
 315     if isinstance(integral, Integral):
 316         return integral.doit()
 317     else:
 318         return integral