Translating the Euler code in scipy's solve_ivp

Question

My code is based on the similarity transformation X=VZ.I simulate the model for transformed equations involving Z by replacing the state space modelAX+BU and RX+SU with transformed equations where X is replaced with VZ and WZ resp then I apply transformation again to obtain X back. I want to execute the following Euler code using solve_ivp. Although I'm able to frame the basic equations but not able to access the solution array at each time step as in Euler to obtain X back from Z. X and Z are simply the state variable matrices of order 4*1. My euler code is :

  for i in range(0,500000):
           if (w== 0 and  X[1]> vdon) or (w==1 and X[0]> 0):
                      zdot=inv(V)*A*V*Z+inv(V)*B*U
                      Z(i+1)=Z(i)+ h*zdot
                      w=1
                      X(i+1)=V*Z(i+1)
           else:
                     zdot=inv(W)*R*W*Z+inv(W)*S*U
                     Z(i+1)=Z(i)+ h*zdot
                    w=0
                    X(i+1)=W*Z(i+1)

where V and W are the eigen vector matrices of A and R respectively obtained using

e1,V=LA.eig(A)
e2,W=LA.eig(R)

In scipy's solve_ivp I define the function as

def conv(t,Z):
   if (w==0 and X[1]>vdon) or (w==1 and X[0]>0):
      zdot=inv(V)*A*V*Z+inv(V)*B*U
      w=1
   else:
       zdot=inv(W)*R*W*Z+inv(W)*S*U
       w=0
return zdot

The if condition shown here doesn't work as expected and has been shown just for code understanding. and I define the solver equations as: w=0 #intially sol= solve_ivp(conv, tspan1,X0) aa1=sol.t bb1=sol.y But I'm unable to define X=W*Z or X=V*Z at each solver time step during run time.

Answer 1

Independent of how that makes sense or not, you have a system with two modes of operation or phases, and the phase change condition is different for each phase, enabling something like hysteresis if the system parameters are right.

So you have the differential systems

def f0(t,Z): return A.dot(Z)+B.dot(U(t))
def f1(t,Z): return R.dot(Z)+S.dot(U(t))

and the change conditions

def g0(t,Z): return vdon(t) - V[1].dot(Z)
def g1(t,Z): return W[0].dot(Z)

The logic of your Euler code then reads as

• if in phase 0: Integrate $Z'=f_0(t,Z)$ as long as $g_0(t,Z)>0$. Over that segment, store $X(t)=VZ(t)$. On $g_0(t,Z)=0$ switch to phase 1.
• if in phase 1: Integrate $Z'=f_1(t,Z)$ as long as $g_1(t,Z)>0$. Over that segment, store $X(t)=WZ(t)$. On $g_1(t,Z)=0$ switch to phase 0.

This recipe has an obvious failure in that if both $g_1<0$ and $g_2<0$ then no phase is admissible, both phases are left as soon as they are entered. One could change the condition from "while $g_k>0$" to "until $g_k$ changes from positive to negative", which can be realized using events with direction

g0.terminal=True; g0.direction=-1;
g1.terminal=True; g1.direction=-1;

and then the switching mechanism is

t,z=t0,z0
sw=0
while t<t_final:
    if sw==0:
        sol = solve_ivp(f0,[t,t_final],z,events=g0)
        Z = sol.y.T
        X = V.dot(Z)
    else:
        sol = solve_ivp(f1,[t,t_final],z,events=g1)
        Z = sol.y.T
        X = W.dot(Z)
    # end if
    t,z=sol.t[-1], Z[-1]
    # collect sol.t, Z, X segments in a list to be later concatenated, 
    # see https://stackoverflow.com/a/54796425 (a discontinuous ODE entering sliding mode)
    # or directly append to large lists
#end while