Piecewise polynomial integration

[drmikecooke]

I have a Piecewise expression:

Ok=[Piecewise((r**(k + 6)/36, (r >= 0) & (r <= 1)), (r**k*(3*r**3 - 12*r**2 + 12*r - 4)**2/36, r <= 2), (r**k*(3*r**3 - 24*r**2 + 60*r - 44)**2/36, r <= 3), (r**k*(r**3 - 12*r**2 + 48*r - 64)**2/36, r <= 4), (0, True))]

or in latex/jax form (in case that is available, a preview option would be helpful):

$$\displaystyle \begin{cases} \frac{r^{k + 6}}{36} & \text{for}\: r \geq 0 \wedge r \leq 1 \\frac{r^{k} \left(3 r^{3} - 12 r^{2} + 12 r - 4\right)^{2}}{36} & \text{for}\: r \leq 2 \\frac{r^{k} \left(3 r^{3} - 24 r^{2} + 60 r - 44\right)^{2}}{36} & \text{for}\: r \leq 3 \\frac{r^{k} \left(r^{3} - 12 r^{2} + 48 r - 64\right)^{2}}{36} & \text{for}\: r \leq 4 \0 & \text{otherwise} \end{cases}$$

This is a simplified product of two cubic B-splines and r**k. (r and k are the obvious symbols)

I 'Ok.integrate((r,0,1))' and get:

Piecewise((oo, Eq(k, -7) | Eq(k, -6) | Eq(k, -5) | Eq(k, -4) | Eq(k, -3) | Eq(k, -2) | Eq(k, -1)), (-0**(k + 7)/(36*(k + 7)) + 1/(36*(k + 7)), Eq(k, -6) | Eq(k, -5) | Eq(k, -4) | Eq(k, -3) | Eq(k, -2) | Eq(k, -1) | ((k > -oo) & (k < oo) & Ne(k, -7))), (oo, True))

or

$$\displaystyle \begin{cases} \infty & \text{for}\: k = -7 \vee k = -6 \vee k = -5 \vee k = -4 \vee k = -3 \vee k = -2 \vee k = -1 \- \frac{0^{k + 7}}{36 \left(k + 7\right)} + \frac{1}{36 \left(k + 7\right)} & \text{for}\: \left(k > -\infty \wedge k < \infty \wedge k \neq -7\right) \vee k = -6 \vee k = -5 \vee k = -4 \vee k = -3 \vee k = -2 \vee k = -1 \\infty & \text{otherwise} \end{cases}$$

To my mind apart from k<-7, this should evaluate to 1/36/(k+6) or something similar. If I was being generous I could go around the singularity an infinite number of times and so on, but that would raise problems only for non-integer k. I assume the 0**(k+7) term is trying to express.

The problematic term for me is the first which immediately give $/infty$, since it excludes the negative integers I am interested in for using B-splines to evaluate simple QM radial energy levels.

[drmikecooke]

Didn't see the Preview-tab, and sorry the formulae are a mess. Ok is:

$$\begin{cases} \frac{r^{k + 6}}{36} & \text{for}\: r \geq 0 \wedge r \leq 1 \ \frac{r^{k} \left(3 r^{3} - 12 r^{2} + 12 r - 4\right)^{2}}{36} & \text{for}\: r \leq 2 \ \frac{r^{k} \left(3 r^{3} - 24 r^{2} + 60 r - 44\right)^{2}}{36} & \text{for}\: r \leq 3 \ \frac{r^{k} \left(r^{3} - 12 r^{2} + 48 r - 64\right)^{2}}{36} & \text{for}\: r \leq 4 \ 0 & \text{otherwise} \end{cases}$$

And the integrated result:

$$ \begin{cases} \infty & \text{for}\: k = -7 \vee k = -6 \vee k = -5 \vee k = -4 \vee k = -3 \vee k = -2 \vee k = -1 \ -\frac{0^{k + 7}}{36 \left(k + 7\right)} + \frac{1}{36 \left(k + 7\right)} & \text{for}\: \left(k > -\infty \wedge k < \infty \wedge k \neq -7\right) \vee k = -6 \vee k = -5 \vee k = -4 \vee k = -3 \vee k = -2 \vee k = -1 \ \infty & \text{otherwise} \end{cases}$$

[drmikecooke]

There seems to be a problem with simplifying piecewise expressions.

$$\begin{cases} x - 2 & \text{for}\: x \geq 2 \wedge x \leq 3 \ 4 - x & \text{for}\: x \geq 3 \wedge x \leq 4 \ 0 & \text{otherwise} \end{cases}$$

simplifies to:

$$\displaystyle \begin{cases} x - 2 & \text{for}\: x \geq 2 \wedge x \leq 3 \ 4 - x & \text{for}\: x \leq 4 \ 0 & \text{otherwise} \end{cases}$$

The former expression is continuous with support [2,4], while the latter is discontinuous with support [$-\infty$,4]. The second conditional of the simplified version is wrong. [Can't seem to get inline latex to render.]

[oscarbenjamin]

Please show complete code for any examples (including imports, defining symbols, etc). Also please say what SymPy version you are using.

I just manually typed out the code for the last example and guessed what functions you called and I don't see what you suggested with SymPy's current master branch:

In [16]: from sympy import *

In [17]: x = Symbol('x')

In [18]: p = Piecewise((x - 2, (x >= 2) & (x <= 3)), (4 - x, (x >= 3) & (x <= 4)), (0, True))

In [19]: p
Out[19]: 
⎧x - 2  for x ≥ 2 ∧ x ≤ 3
⎪                        
⎨4 - x  for x ≥ 3 ∧ x ≤ 4
⎪                        
⎩  0        otherwise    

In [20]: simplify(p)
Out[20]: 
⎧x - 2  for x ≥ 2 ∧ x ≤ 3
⎪                        
⎨4 - x  for x ≤ 4 ∧ x > 3
⎪                        
⎩  0        otherwise 

[drmikecooke]

I'm using jupyterlab on a raspberry pi. I hope this is clear enough (a slightly different problem, but similar result):

import sympy as sp
import sympy.functions.special.bsplines as spbs

d,N=1,4
sknots=tuple(n for n in range(-d,N+d+1))
x=sp.Symbol('x')
BS=sp.Matrix(spbs.bspline_basis_set(d,sknots,x))

BS[2]

$$\displaystyle \begin{cases} x - 1 & \text{for}\: x \geq 1 \wedge x \leq 2 \ 3 - x & \text{for}\: x \geq 2 \wedge x \leq 3 \ 0 & \text{otherwise} \end{cases}$$

sp.simplify(BS[2])

$$\displaystyle \begin{cases} x - 1 & \text{for}\: x \geq 1 \wedge x \leq 2 \ 3 - x & \text{for}\: x \leq 3 \ 0 & \text{otherwise} \end{cases}$$

sp.__version__

'1.12.1'

Rather than using simplify, I seem to have started getting reasonable results by using piecewise_fold to combine piecewise expressions.

[oscarbenjamin]

This is already fixed in SymPy 1.13. Please try updating SymPy to the latest version.

The fix was 3d7c5c8d6a742798e2b3b4ef5704ac645be0d3c4 from gh-25649

[drmikecooke]

That's better:

import sympy as sp
import sympy.functions.special.bsplines as spbs

d,N=1,4
sknots=tuple(n for n in range(-d,N+d+1))
x=sp.Symbol('x')
BS=sp.Matrix(spbs.bspline_basis_set(d,sknots,x))

BS[2]

$$\displaystyle \begin{cases} x - 1 & \text{for}\: x \geq 1 \wedge x \leq 2 \ 3 - x & \text{for}\: x \geq 2 \wedge x \leq 3 \ 0 & \text{otherwise} \end{cases}$$

sp.simplify(BS[2])

$$\displaystyle \begin{cases} x - 1 & \text{for}\: x \geq 1 \wedge x \leq 2 \ 3 - x & \text{for}\: x \leq 3 \wedge x > 2 \ 0 & \text{otherwise} \end{cases}$$

sp.__version__

'1.13.3'

[drmikecooke]

Thanks.