Hi Burcin,
On Wed, Jul 8, 2009 at 11:02 AM, Burcin Erocal<[email protected]> wrote:
> If you share your code creating an SFunction called 'integral', I can
> produce patches for Sage and pynac to special case that function when
> computing derivatives. (Much like what is done for Order in the pynac
> code right now.)
Thanks for looking into this. I am attaching prototype integration
class and integration code for generalized function that you can play with.
Sub-classing integration from SFunction class, makes few thing really
easy. For example numerical approximation .n() works like a charm
------
sage: integrate( exp(x)/x, x, 1, 10)
integrate(e^x/x, x, 1, 10)
sage: integrate( exp(x)/x, x, 1, 10).n()
2490.33385842557
------
Cheers,
Golam
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
"""
Symbolic Integration
AUTHORS:
* Golam Mortuza Hossain (2009-07-08)
"""
##############################################################################
#
# Copyright (C) 2009 Golam Mortuza Hossain <[email protected]>
#
# Distributed under the terms of the GNU General Public License (GPL v2+)
# http://www.gnu.org/licenses/
#
##############################################################################
from sage.symbolic.ring import SR
from sage.symbolic.function import SFunction, sfunctions_funcs
from generalized_function import generalized_function_integrator
def maxima_integrator(*args, **kwds):
"""
Maxima Integrator
"""
from sage.calculus.calculus import integral as maxima_integral
# Ask maxima_integral to raise NotImplementedError if it fails
# to integrate, so that other algorithm can be tried.
kwds['raise_error'] = True
return maxima_integral(*args, **kwds)
def sympy_integrator(*args, **kwds):
"""
Sympy Integrator
"""
from sage.calculus.calculus import integral as sympy_integral
kwds['algorithm'] = "sympy"
# Check whether sympy could evaluate the integral
try:
return sympy_integral(*args, **kwds)
except:
raise NotImplementedError, "Sympy failed to integrate"
#####################################################
# List of integrator
#####################################################
integrator_list = \
[
generalized_function_integrator,
maxima_integrator,
sympy_integrator,
]
class SymbolicIntegration(SFunction):
def __init__(self, *args, **kwds):
"""
EXAMPLES::
"""
kwds['built_in_function'] = True
for name in sfunctions_funcs:
if hasattr(self, "_%s_"%name):
kwds['%s_func'%name] = getattr(self, "_%s_"%name)
SFunction.__init__(self, "integrate", *args, **kwds)
#from sage.symbolic.pynac import register_symbol
#register_symbol(self, self._conversions)
def _eval_(self, *args, **kwds):
"""
Returns the results of symvolic evaluation of the integral
EXAMPLES::
"""
for integrator in integrator_list:
try:
return integrator(*args, **kwds)
except NotImplementedError:
pass
return None
def _evalf_(self, f, x, a, b, **kwds):
"""
Returns the results of numerical evaluation of the inetgral
EXAMPLES::
"""
from sage.calculus.calculus import nintegral
return nintegral(f, x, a, b)[0]
def _print_latex_(self, f, x, a=None, b=None, **kwds):
r"""
Return LaTeX expression for integration of a symbolic function.
EXAMPLES::
sage: from sage.calculus.calculus import _integrate_latex_
sage: var('x,a,b')
(x, a, b)
sage: f(x) = function('f',x)
sage: _integrate_latex_(f(x),x)
'\\int f\\left(x\\right)\\,{d x}'
sage: _integrate_latex_(f(x),x,a,b)
'\\int_{a}^{b} f\\left(x\\right)\\,{d x}'
"""
from sage.misc.latex import latex
# Check whether its a definite integral
if a is not None:
return "\\int_{%s}^{%s} %s\\,{d %s}"%(latex(a), latex(b), latex(f), latex(x))
# Typeset as indefinite integral
return "\\int %s\\,{d %s}"%(latex(f), latex(x))
def _diff_derivative_(self, F, var):
"""
Return derivative of symbolic integration
EXAMPLES::
"""
op = F.operands()
f = op[0]
x = op[1]
a = None
b = None
if len(op) == 4:
a = op[2]
b = op[3]
if a is None:
if x._derivative(var) != 0:
return f*x._derivative(var)
else:
return f._derivative(arg).integral(x)
elif x._derivative(var) != 0: # dummy variable is same as arg
return 0
else:
return f.subs(x==b)*b._derivative(var) \
- f.subs(x==a)*a._derivative(var) \
+ f._derivative(var).integral(x,a,b)
integral = integrate = SymbolicIntegration()
"""
Symbolic Integration of Generalized Functions
This module performs symbolic integration when its integrand
contains at least one generalized function.
AUTHORS:
* Golam Mortuza Hossain (2009-06-25)
"""
##############################################################################
#
# Copyright (C) 2009 Golam Mortuza Hossain <[email protected]>
#
# Distributed under the terms of the GNU General Public License (GPL v2+)
# http://www.gnu.org/licenses/
#
##############################################################################
import math, operator
from sage.functions.generalized import dirac_delta, heaviside, unit_step
from sage.rings.complex_interval_field import ComplexIntervalField
from sage.symbolic.operators import FDerivativeOperator
from sage.symbolic.pynac import I
from sage.symbolic.ring import SR
from sage.misc.functional import integrate
##################################################################
#
# Symbolic integration algorithm where the integrand involves
# at least one generalized function such as dirac_delta,
# heaviside, or unit_step.
#
##################################################################
# List of generalized functions
gf_list = [dirac_delta, heaviside, unit_step]
# List of step functions
sf_list = [heaviside, unit_step]
def generalized_function_integrator(f, x, a=None, b=None, **kwds):
"""
Integrate expression involving generalized functions
"""
# Expand the integrand
f = f.expand()
# If f is in sum form "f1 + f2 + ..." then evalute them separately.
f_op = f.operator()
if f_op is operator.add:
ans = SR(0)
flist = list(f.iterator())
for subf in flist:
if a is None:
ans = ans + subf.integral(x, **kwds)
else:
ans = ans + subf.integral(x, a, b, **kwds)
return ans
return _evaluate_integral(f, x, a, b, **kwds)
gf_integral = generalized_function_integrator
def _evaluate_integral(f, x, a, b, **kwds):
"""
Computes integral when the integrand contains a
generalized function
"""
subf = _get_gf_from_integrand(f)
if subf is None:
raise NotImplementedError, "No generalized function " + \
"found in the integrand"
elif x not in subf.variables():
raise NotImplementedError, "Variable of integration " + \
"is not in args of generalized function"
# Remaining integrand
#fx = f.subs_expr(subf==1)
fx = f/subf
op = subf.operator()
# Indefinite integral
if a is None:
if op in sf_list:
return _heaviside_indefinite_integral(fx, subf, x)
else:
return _dirac_delta_indefinite_integral(fx, subf, x)
# Definite integral
else:
if op in sf_list:
return _heaviside_definite_integral(fx, subf, x, a, b)
else:
return _dirac_delta_definite_integral(fx, subf, x, a, b)
#############################################################
# Integration involving Dirac delta function
#############################################################
def _dirac_delta_indefinite_integral(fx, subf, x):
"""
Computes an indefinite integral when the integrand contains
dirac_delta function
EXAMPLES::
sage: from sage.symbolic.integration import gf_integral
sage: gf_integral(dirac_delta(x), x)
heaviside(x)
"""
op = subf.operator()
# Check for derivative operator
if isinstance(op, FDerivativeOperator):
ndiff = len(op.parameter_set())
else:
ndiff = 0
# The argument of dirac_delta(g(x))
gx = subf.operands()[0]
roots = _get_all_real_roots(gx,x)
# -=: Evaluation algorithm :=-
# Using integration by parts once
# "integrate(f(x)*D[n]dirac_delta(g(x)), x)" evaluates
# to "fx*D[n-1]dirac_delta(g(x))"
# - "integrate( f(x).diff(x)*D[n-1]dirac_delta(g(x)), x)"
# where "fx = f(x)/g(x).diff(x)"
# This process is repeated until it reduces to the form
# "integrate(fx*dirac_delta(g(x)), x)". Using the definition
# of Dirac delta this integration then evaluates to
# "sum(integrand, over roots of g(x)==0)" where
# "integrand = fx/g(x).diff(x).abs()|x=x_i"
ans = SR(0)
gderiv = gx.diff(x)
# if ndiff > 0 then use integration by parts
if ndiff > 0:
for j in range(1,ndiff+1):
fx = fx/gderiv
dd_factor = dirac_delta(x).diff(x,ndiff-j).subs(x==gx)
ans = ans + fx*dd_factor
# We apply derivative for next steps
fx = -fx.diff(x)
# the Dirac delta part
integrand = fx/gderiv.abs()
for xi in roots:
ans = ans + heaviside(x-xi)*integrand.subs(x==xi)
return ans
def _dirac_delta_definite_integral(fx, subf, x, a, b):
"""
Computes definite integral when integrand contains
dirac_delta function
"""
# Use indefinite integral
indef = _dirac_delta_indefinite_integral(fx, subf, x)
return indef.subs(x==b)-indef.subs(x==a)
#############################################################
# Integration involving Heaviside step functions
#############################################################
def _heaviside_indefinite_integral(fx, subf, x):
"""
Computes an indefinite integral when the integrand contains
heaviside step function
"""
# Using integration by parts
if fx == 1:
fxint = x
else:
fxint = fx.integral(x)
if "%s"%(fxint.operator()) == "integrate":
raise NotImplementedError
gx = subf.operands()[0]
newfx = gx.diff(x)*fxint
newsubf = dirac_delta(gx)
term2 = _dirac_delta_indefinite_integral(newfx, newsubf, x)
return heaviside(gx)*fxint - term2
def _heaviside_definite_integral(fx, subf, x, a, b):
"""
Computes definite integral when the integrand contains
heaviside step function
"""
gx = subf.operands()[0]
# If gx < 0 at both limits then certain heuristic can
# determine whether integration is zero
if bool(gx.subs(x==a) <= 0) and bool(gx.subs(x==b) <= 0):
# if gx is linear
if gx.diff(x,2) == 0:
return SR(0)
# Check whether any root lies in between the limits
roots = _get_all_real_roots(gx,x)
no_zeroes = True
for xi in roots:
if bool(xi > a) and bool(xi < b):
no_zeroes = False
break
if no_zeroes:
return SR(0)
# Now try using indefinite integral
indef = _heaviside_indefinite_integral(fx, subf, x)
return indef.subs(x==b)-indef.subs(x==a)
#############################################################
# Helper functions
#############################################################
def _get_gf_from_integrand(f):
"""
Get the generalized function from the integrand. If not
found then return None. Dirac delta gets the precedence.
"""
f_op = f.operator()
if f_op in gf_list:
return f
elif f_op is operator.mul:
flist = list(f.iterator())
elif isinstance(f_op, FDerivativeOperator):
flist = [f]
else:
return None
# Go through the list to see whether it has generalized functions in it
stepf = None
for subf in flist:
op = subf.operator()
# Check also for derivative operator of Dirac delta
if op is dirac_delta or isinstance(op, FDerivativeOperator) \
and op.function() is dirac_delta:
return subf
elif op in sf_list:
stepf = subf
return stepf
def _get_all_real_roots(gx, x):
"""
Get all real roots of the equation gx==0
"""
# If gx is linear (most common usage) then avoid
# using .roots() as it invokes maxima and that slows down
if gx.diff(x,2) == 0:
b = gx.subs(x==0)
a = gx.subs(x==1)-b
roots = [(-b/a, 1)]
else:
try:
roots = (gx==0).roots(x)
except RuntimeError:
raise NotImplementedError
valid_roots = []
for rt in roots:
n = rt[1]
# Multiplicity
if n > 1:
raise ValueError, "\'%s\' has root with multiplicity > 1"%(repr(gx))
xi = rt[0]
if _root_is_complex(xi):
pass
else:
valid_roots.append(xi)
return valid_roots
def _root_is_complex(x):
"""
Check whether the root is complex
"""
try:
approx_x = ComplexIntervalField()(x)
if bool(approx_x.imag() == 0): # x is real
return False
else:
return True
except:
pass
#if x.has(I) is True:
# return True
return False
