Sympy fails to integrate the product of a piecewise continuous f

ghz 8hours ago ⋅ 5 views

Sympy fails to integrate the product of a piecewise continuous function and a complex function across the discontinuity

If I do:

from sympy import *
x, L = symbols('x L', real=True)
f = Piecewise((1, x<=0), (-1, x<1), (0, True))
g = exp(-x * 1j)
integrate(f * g, (x, 0, L))

I get:

Piecewise((1.0*I*exp(-1.0*I*L) - 1.0*I, L <= 0), (-1.0*I*exp(-1.0*I*L) + 1.0*I, L < 1), (-1.0*I*exp(-1.0*I) + 1.0*I, True))

But if I change the last line to:

integrate(f*g, (x, L/2, L))

I get:

Integral(Piecewise((exp(-1.0*I*x), x <= 0), (-exp(-1.0*I*x), x < 1), (0, True)), (x, L/2, L))

Any insight would be appreciated!

Answer

The behavior you're seeing is due to the fact that sympy is not able to automatically handle the Piecewise function in the same way when you adjust the limits of integration. Let me explain what's going on in more detail.

First Scenario: integrate(f * g, (x, 0, L))

In the first case, where you're integrating from 0 to L, the Piecewise function f behaves predictably within that range. Since x is always positive in this range (i.e., x <= 0 doesn't hold for x in [0, L]), SymPy directly integrates the expression. The result involves a combination of terms, and the final expression is simplified correctly.

The output you see:

Piecewise((1.0*I*exp(-1.0*I*L) - 1.0*I, L <= 0), (-1.0*I*exp(-1.0*I*L) + 1.0*I, L < 1), (-1.0*I*exp(-1.0*I) + 1.0*I, True))

is a valid result, though it may be a bit more complex due to the conditional structure in Piecewise. The expression inside the Piecewise is evaluating the complex integrals, and the simplification happens appropriately.

Second Scenario: integrate(f * g, (x, L/2, L))

When you change the integration limits to (L/2, L), sympy has to consider how the Piecewise function behaves within this new range. The issue arises because the function f(x) involves three separate cases:

  • 1 when x <= 0
  • -1 when x < 1
  • 0 for all other values (True case)

The integral might be more difficult for SymPy to handle as it now needs to keep track of these conditions across the entire interval of integration. The resulting output keeps the Piecewise structure as-is:

Integral(Piecewise((exp(-1.0*I*x), x <= 0), (-exp(-1.0*I*x), x < 1), (0, True)), (x, L/2, L))

This suggests that sympy was unable to automatically compute the integral with the Piecewise structure and opted to return an unevaluated integral with the piecewise function still intact.

Possible Explanation:

  1. Piecewise Handling: SymPy integrates Piecewise functions by considering the cases defined in the function. However, when you change the integration limits, the function may span over different conditions in the piecewise structure, which complicates the symbolic integration.

  2. Simplification and Integration: When SymPy can't easily simplify the integral in the given range, it opts to return an unevaluated integral instead of computing it directly. This happens because f(x) has conditions dependent on x, and the integration limits in the second case may lead to discontinuities or mixed conditions that are harder to resolve symbolically.

Solutions and Approaches:

  1. Manual Case Analysis: If the integral can't be automatically solved symbolically, you can break down the integration into pieces manually by analyzing the piecewise function at different ranges (depending on the value of L/2 and L):

    # Define the different cases of f(x)
# You would need to manually compute the integrals for the different regions.

  2. Numerical Integration: If a symbolic result is not necessary, you can use numerical integration. You can use scipy.integrate.quad or SymPy’s integrate with numerical=True to perform numerical integration over the specified range:

    from sympy import integrate, symbols, exp, Piecewise
from scipy.integrate import quad
import numpy as np

x, L = symbols('x L', real=True)
f = Piecewise((1, x<=0), (-1, x<1), (0, True))
g = exp(-x * 1j)

# Define the function for numerical integration
def integrand(x_val, L_val):
    return f.subs(x, x_val) * g.subs(x, x_val).evalf()

result, error = quad(integrand, 0, float(L))

  3. Using subs() to evaluate at specific values: If you don't want to break things down manually but still need a specific result, you can substitute specific values into the function before integration, effectively turning the piecewise function into a regular function for integration.

    f_at_L = f.subs(x, L)

Conclusion:

The issue you're facing is related to how SymPy handles Piecewise functions when integration limits cross over different cases in the Piecewise. In your second case, SymPy keeps the integral unevaluated because it can't simplify it automatically. You can work around this by either doing case-based integration manually or using numerical integration methods.

#py