Hi Burcin, On Wed, Jul 8, 2009 at 11:02 AM, Burcin Erocal<bur...@erocal.org> 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 sage-devel@googlegroups.com To unsubscribe from this group, send email to sage-devel-unsubscr...@googlegroups.com 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 <gmhoss...@gmail.com> # # 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 <gmhoss...@gmail.com> # # 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