Problems using numpy.piecewise

Question

1. The core problem and question

I will provide an executable example below, but let me first walk you through the problem first.

I am using solve_ivp from scipy.integrate to solve an initial value problem (see documentation). In fact I have to call the solver twice, to once integrate forward and once backward in time. (I would have to go unnecessarily deep into my concrete problem to explain why this is necessary, but please trust me here--it is!)

sol0 = solve_ivp(rhs,[0,-1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)
sol1 = solve_ivp(rhs,[0, 1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)

Here rhs is the right hand side function of the initial value problem y(t) = rhs(t,y). In my case, y has six components y[0] to y[5]. y0=y(0) is the initial condition. [0,±1e8] are the respective integration ranges, one forward and the other backward in time. rtol and atol are tolerances.

Importantly, you see that I flagged dense_output=True, which means that the solver does not only return the solutions on the numerical grids, but also as interpolation functions sol0.sol(t) and sol1.sol(t).

My main goal now is to define a piecewise function, say sol(t) which takes the value sol0.sol(t) for t<0 and the value sol1.sol(t) for t>=0. So the main question is: How do I do that?

I thought that numpy.piecewise should be tool of choice to do this for me. But I am having trouble using it, as you will see below, where I show you what I tried so far.

---

2. Example code

The code in the box below solves the initial value problem of my example. Most of the code is the definition of the rhs function, the details of which are not important to the question.

import numpy as np
from scipy.integrate import solve_ivp

# aux definitions and constants
sin=np.sin; cos=np.cos; tan=np.tan; sqrt=np.sqrt; pi=np.pi;  
c  = 299792458
Gm = 5.655090674872875e26    

# define right hand side function of initial value problem, y'(t) = rhs(t,y)
def rhs(t,y):
    p,e,i,Om,om,f = y
    sinf=np.sin(f); cosf=np.cos(f); Q=sqrt(p/Gm); opecf=1+e*cosf;        

    R = Gm**2/(c**2*p**3)*opecf**2*(3*(e**2 + 1) + 2*e*cosf - 4*e**2*cosf**2)
    S = Gm**2/(c**2*p**3)*4*opecf**3*e*sinf         

    rhs    = np.zeros(6)
    rhs[0] = 2*sqrt(p**3/Gm)/opecf*S
    rhs[1] = Q*(sinf*R + (2*cosf + e*(1 + cosf**2))/opecf*S)
    rhs[2] = 0
    rhs[3] = 0
    rhs[4] = Q/e*(-cosf*R + (2 + e*cosf)/opecf*sinf*S)
    rhs[5] = sqrt(Gm/p**3)*opecf**2 + Q/e*(cosf*R - (2 + e*cosf)/opecf*sinf*S)

    return rhs

# define initial values, y0
y0=[3.3578528933149297e13,0.8846,2.34921,3.98284,1.15715,0]

# integrate twice from t = 0, once backward in time (sol0) and once forward in time (sol1)
sol0 = solve_ivp(rhs,[0,-1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)
sol1 = solve_ivp(rhs,[0, 1e8],y0,rtol=10e-12,atol=10e-12,dense_output=True)

The solution functions can be addressed from here by sol0.sol and sol1.sol respectively. As an example, let's plot the 4th component:

from matplotlib import pyplot as plt

t0 = np.linspace(-1,0,500)*1e8
t1 = np.linspace( 0,1,500)*1e8
plt.plot(t0,sol0.sol(t0)[4])
plt.plot(t1,sol1.sol(t1)[4])
plt.title('plot 1')
plt.show()

---

3. Failing attempts to build piecewise function

3.1 Build vector valued piecewise function directly out of sol0.sol and sol1.sol

def sol(t): return np.piecewise(t,[t<0,t>=0],[sol0.sol,sol1.sol])
t = np.linspace(-1,1,1000)*1e8
print(sol(t))

This leads to the following error in piecewise in line 628 of .../numpy/lib/function_base.py:

TypeError: NumPy boolean array indexing assignment requires a 0 or 1-dimensional input, input has 2 dimensions

I am not sure, but I do think this is because of the following: In the documentation of piecewise it says about the third argument:

funclistlist of callables, f(x,*args,**kw), or scalars

[...]. It should take a 1d array as input and give an 1d array or a scalar value as output. [...].

I suppose the problem is, that the solution in my case has six components. Hence, evaluated on a time grid the output would be a 2d array. Can someone confirm, that this is indeed the problem? Since I think this really limits the usefulness of piecewiseby a lot.

3.2 Try the same, but just for one component (e.g. for the 4th)

def sol4(t): return np.piecewise(t,[t<0,t>=0],[sol0.sol(t)[4],sol1.sol(t)[4]])
t = np.linspace(-1,1,1000)*1e8
print(sol4(t))

This results in this error in line 624 of the same file as above:

ValueError: NumPy boolean array indexing assignment cannot assign 1000 input values to the 500 output values where the mask is true

Contrary to the previous error, unfortunately here I have so far no idea why it is not working.

3.3 Similar attempt, however first defining functions for the 4th components

def sol40(t): return sol0.sol(t)[4]
def sol41(t): return sol1.sol(t)[4]
def sol4(t): return np.piecewise(t,[t<0,t>=0],[sol40,sol41])
t = np.linspace(-1,1,1000)
plt.plot(t,sol4(t))
plt.title('plot 2')
plt.show()

Now this does not result in an error, and I can produce a plot, however this plot doesn't look like it should. It should look like plot 1 above. Also here, I so far have no clue what is going on.

Am thankful for help!

Answer 1

You can take a look to numpy.piecewise source code. There is nothing special in this function so I suggest to do everything manually.

def sol(t):
    ans = np.empty((6, len(t)))
    ans[:, t<0] = sol0.sol(t[t<0])
    ans[:, t>=0] = sol1.sol(t[t>=0])
    return ans

Regarding your failed attempts. Yes, piecewise excpect functions return 1d array. Your second attempt failed because documentation says that funclist argument should be list of functions or scalars but you send the list of arrays. Contrary to the documentation it works even with arrays, you just should use the arrays of the same size as t < 0 and t >= 0 like:

def sol4(t): return np.piecewise(t,[t<0,t>=0],[sol0.sol(t[t<0])[4],sol1.sol(t[t>=0])[4]])