Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 2020-06-19 15:35:27 +0100

Florentin Jaffredo gravatar image

Hello,

I don't know why this is the case, but here is a dirty workaround:

from operator import eq
def integrate_delta(E, x, a, b):
    """
    Integrate E over [a, b] with respect to x, but tries to find a dirac delta first.
    """

    t = var('temporary_variable')
    E *= t
    # try to find a pattern of the shape f(x)*delta(g(x))
    if not E.find(SR.wild(0)*dirac_delta(SR.wild(1))):
        return E.integrate(x, a, b) # if not found, integrate normaly

    match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
    f, g = match[SR.wild(0)], match[SR.wild(1)] # extract the functions f and g

    roots = list(g.solve(x, solution_dict=True, to_poly_solve="force")) # Solve g(x)==0
    roots_with_multiplicity = []

    for r in roots:
        # sometimes there is an unknown integer in the solutions
        extra_variables = set(r[x].variables()).difference(g.variables())
        if len(extra_variables)==0: # if not, keep the root
            roots_with_multiplicity.append(r[x])
            continue
        if len(extra_variables)>1: # if there are more, we don't know how to solve, fall back.
            return E.integrate(x, a, b)
        z = extra_variables.pop()

        # if there is, we must make sure to keep only the values that stay inside the bounds
        a0, b0 = [t1 for t1 in solve_ineq([r[x]>=a, r[x]<=b], (z,))
               if all(t2.operator() is not eq for t2 in t1)][0]
        roots_with_multiplicity += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

    # returns sum(f(r)/abs(g'(r))) where the sum is over the roots of g.
    return sum(f.subs({x: r, t: 1})/abs(g.diff(x).subs({x: r, t: 1})) for r in roots_with_multiplicity if r>=a and r<=b)

This function will in most case call the usual integrate method, but will first try to find a Dirac delta. The code is absolutely ugly, and there is no guarantee whatsoever that it always works. But I tried it with some trigonometric functions, and it seems fine. Use it like this:

sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2

Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
2/3*sqrt(3)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
4*sqrt(3)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
0
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
1/2
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
1

Hello,

I don't know why this is the case, but here is a dirty workaround:

from operator import eq
def integrate_delta(E, x, a, b):
    """
    Integrate E over [a, b] with respect to x, but tries to find a dirac delta first.
    """

    t = var('temporary_variable')
    E *= t
    # try to find a pattern of the shape f(x)*delta(g(x))
    if not E.find(SR.wild(0)*dirac_delta(SR.wild(1))):
(E.find(SR.wild(0)*dirac_delta(SR.wild(1))) or E.find(dirac_delta(SR.wild(1)))):
        return E.integrate(x, a, b) # if not found, integrate normaly

    match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
    if match is None:
        match = E.match(dirac_delta(SR.wild(1)))
        match[SR.wild(0)] = 1

    f, g = match[SR.wild(0)], match[SR.wild(1)] # extract the functions f and g

    roots = list(g.solve(x, solution_dict=True, to_poly_solve="force")) # Solve g(x)==0
    roots_with_multiplicity = []

    for r in roots:
        # sometimes there is an unknown integer in the solutions
        extra_variables = set(r[x].variables()).difference(g.variables())
        if len(extra_variables)==0: # if not, keep the root
            roots_with_multiplicity.append(r[x])
            continue
        if len(extra_variables)>1: # if there are more, we don't know how to solve, fall back.
            return E.integrate(x, a, b)
        z = extra_variables.pop()

        # if there is, we must make sure to keep only the values that stay inside the bounds
        a0, b0 = [t1 for t1 in solve_ineq([r[x]>=a, r[x]<=b], (z,))
               if all(t2.operator() is not eq for t2 in t1)][0]
        roots_with_multiplicity += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

    # returns sum(f(r)/abs(g'(r))) sum(f(r)*1/abs(g'(r))) where the sum is over the roots of g.
    return sum(f.subs({x: r, t: 1})/abs(g.diff(x).subs({x: r, t: 1})) r})/abs(g.diff(x).subs({x: r})) for r in roots_with_multiplicity if r>=a and r<=b)

This function will in most case call the usual integrate method, but will first try to find a Dirac delta. The code is absolutely ugly, and there is no guarantee whatsoever that it always works. But I tried it with some trigonometric functions, and it seems fine. Use it like this:

sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2

Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
2/3*sqrt(3)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
4*sqrt(3)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
0
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
1/2
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
1

Hello,

I don't really know why this is the case, but here simplifying Dirac deltas is a dirty workaround:hard problem, since solving an equation is always needed.

Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.

from operator import eq
def integrate_delta(E, x, a, b):
    """
    Integrate E over [a, b] with respect to x, but tries to find a dirac delta first.
    """
    # try to find a pattern of the shape f(x)*delta(g(x))
    if not (E.find(SR.wild(0)*dirac_delta(SR.wild(1))) or E.find(dirac_delta(SR.wild(1)))):
        return E.integrate(x, a, b) # if not found, integrate normaly

    match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
    if match is None:
        match = E.match(dirac_delta(SR.wild(1)))
        match[SR.wild(0)] = 1

    f, g = match[SR.wild(0)], match[SR.wild(1)] # extract the functions f and g

    roots = list(g.solve(x, solution_dict=True, to_poly_solve="force")) # Solve g(x)==0
    roots_with_multiplicity = []

    for r in roots:
        # sometimes there is an unknown integer in the solutions
        extra_variables = set(r[x].variables()).difference(g.variables())
        if len(extra_variables)==0: # if not, keep the root
            roots_with_multiplicity.append(r[x])
            continue
        if len(extra_variables)>1: # if there are more, we don't know how to solve, fall back.
            return E.integrate(x, a, b)
        z = extra_variables.pop()

        # if there is, we must make sure to keep only the values that stay inside the bounds
        a0, b0 = [t1 for t1 in solve_ineq([r[x]>=a, r[x]<=b], (z,))
               if all(t2.operator() is not eq for t2 in t1)][0]
        roots_with_multiplicity += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

    # returns sum(f(r)*1/abs(g'(r))) where the sum is over the roots of g.
    return sum(f.subs({x: r})/abs(g.diff(x).subs({x: r})) for r in roots_with_multiplicity if r>=a and r<=b)

This function will in most case call the usual integrate method, but will first try to find a Dirac delta. The code is absolutely ugly, and there is no guarantee whatsoever that it always works. But I tried it with some trigonometric functions, and it seems fine. Use it like this:

sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2

Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
2/3*sqrt(3)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
4*sqrt(3)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
0
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
1/2
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
1

Hello,

I don't really know why this is the case, but simplifying Dirac deltas is a hard problem, since solving an equation is always needed.

Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.

from operator import eq
def integrate_delta(E, x, a, b):
    """
    Integrate E over [a, b] with respect to x, but tries to find a dirac delta first.
    """
    # try to find a pattern of the shape f(x)*delta(g(x))
    if not (E.find(SR.wild(0)*dirac_delta(SR.wild(1))) or E.find(dirac_delta(SR.wild(1)))):
        return E.integrate(x, a, b) # if not found, integrate normaly

    # If expression is a sum, call recursively
    if E.operator() is sage.symbolic.operators.add_vararg:
        return sum(integrate_delta(op, x, a, b) for op in E.operands())

    match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
    if match is None:
        match = E.match(dirac_delta(SR.wild(1)))
        match[SR.wild(0)] = 1

    f, g = match[SR.wild(0)], match[SR.wild(1)] # extract the functions f and g

    roots = list(g.solve(x, solution_dict=True, to_poly_solve="force")) # Solve g(x)==0
    roots_with_multiplicity = []

    for r in roots:
        # sometimes there is an unknown integer in the solutions
        extra_variables = set(r[x].variables()).difference(g.variables())
        if len(extra_variables)==0: # if not, keep the root
            roots_with_multiplicity.append(r[x])
            continue
        if len(extra_variables)>1: # if there are more, we don't know how to solve, fall back.
            return E.integrate(x, a, b)
        z = extra_variables.pop()

        # if there is, we must make sure to keep only the values that stay inside the bounds
        a0, b0 = [t1 for t1 in solve_ineq([r[x]>=a, r[x]<=b], (z,))
               if all(t2.operator() is not eq for t2 in t1)][0]
        roots_with_multiplicity += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

    # returns sum(f(r)*1/abs(g'(r))) where the sum is over the roots of g.
    return sum(f.subs({x: r})/abs(g.diff(x).subs({x: r})) for r in roots_with_multiplicity if r>=a and r<=b)

This function will in most case call the usual integrate method, but will first try to find a Dirac delta. The code is absolutely ugly, and there is no guarantee whatsoever that it always works. Use it like this:

sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2

Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
2/3*sqrt(3)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
4*sqrt(3)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
0
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
1/2
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
1

Hello,

I don't really know why this is the case, but simplifying Dirac deltas is a hard problem, since solving an equation is always needed.

Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.

from operator import eq
def integrate_delta(E, x, a, b):
    """
    Integrate E over [a, b] with respect to x, but tries to find a dirac delta first.
    """
    # try to find a pattern of the shape f(x)*delta(g(x))
    if not (E.find(SR.wild(0)*dirac_delta(SR.wild(1))) or E.find(dirac_delta(SR.wild(1)))):
        return E.integrate(x, a, b) # if not found, integrate normaly

    # If expression is a sum, call recursively
    if E.operator() is sage.symbolic.operators.add_vararg:
        return sum(integrate_delta(op, x, a, b) for op in E.operands())

    match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
    if match is None:
        match = E.match(dirac_delta(SR.wild(1)))
        match[SR.wild(0)] = 1

    f, g = match[SR.wild(0)], match[SR.wild(1)] # extract the functions f and g

    roots = list(g.solve(x, solution_dict=True, to_poly_solve="force")) # Solve g(x)==0
    roots_with_multiplicity explicit_roots = []

    for r in roots:
        # sometimes there is an unknown integer in the solutions
        extra_variables = set(r[x].variables()).difference(g.variables())
        if len(extra_variables)==0: # if not, keep the root
            roots_with_multiplicity.append(r[x])
explicit_roots.append(r[x])
            continue
        if len(extra_variables)>1: # if there are more, we don't know how to solve, fall back.
            return E.integrate(x, a, b)
        z = extra_variables.pop()

        # if there is, we must make sure to keep only the values that stay inside the bounds
        a0, b0 = [t1 for t1 in solve_ineq([r[x]>=a, r[x]<=b], (z,))
               if all(t2.operator() is not eq for t2 in t1)][0]
        roots_with_multiplicity explicit_roots += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

    # returns sum(f(r)*1/abs(g'(r))) where the sum is over the roots of g.
    return sum(f.subs({x: r})/abs(g.diff(x).subs({x: r})) for r in roots_with_multiplicity explicit_roots if r>=a and r<=b)

This function will in most case call the usual integrate method, but will first try to find a Dirac delta. The code is absolutely ugly, and there is no guarantee whatsoever that it always works. Use it like this:

sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2

Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
2/3*sqrt(3)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
4*sqrt(3)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
0
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
1/2
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
1
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59