# Unexpected solutions solving an ODE using odeint

I am trying to solve a system of 8 coupled differential equations using scipy's odeint. I have already written my code and it runs fine, but the solutions I get are completely different from what I expected. Originally, I would write all these equations within a loop, but since I was having problems with them I wrote them all separately.

The Zd argument that I need to pass is an exponentially decreasing variable the same size as my time vector. However, since odeint only calculates one sample at a time, I need to give it the right Zd sample for each time, so I use a copy of t to get it inside the function. The differential equations have four terms each and they need to be turned on/off at different times, so I use the step function there.

This system of equations is supposed to represent a sort of energy going through adjoining systems, so the eight Z functions are the evolution of energy within each one of these systems over time. At t=0, the energy enters the first of these systems and the others have no energy whatsoever. Zd represents the input energy into the system, which changes over time (it decays exponentially as a result of scattering) and needs to be turned on/off as energy enters/leaves each layer. By the end of my time window, I expected to have little or no energy at all. However, my resulting Zs grow almost exponentially and this doesn't make sense.

I have checked it over and over again but I can't see what is wrong in it. The signs and times in the equations are correct. Any idea on what is wrong here?

This is my code and a plot of the results:

def ODE_solver(Z,t,Q,Zd,Z0,dts,ts,cf,t_copy):

import numpy as np

# Define functions to solve for:
[Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8] = Z
# Define time intervals:
[dt1, dt2, dt3, dt4, dt5, dt6, dt7, dt8] = dts
# Define cumulative times:
[t1, t2, t3, t4, t5, t6, t7, t8] = ts
# Q factors:
[Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8] = Q

# Set some boundary conditions:
Z9=0; dt9=1/4 # (It doesn't matter which value I have here, since Z9 = 0)
t0=0; dt0=1/4;

# Find index of the time sample that is being used at the moment:
t_ind=(np.abs(t-t_copy)).argmin()
Zd=Zd[t_ind]

# j=1
dZ1dt= ((1/(4*dt0)) * Z0 * step_fun (t,t0)) + ((1/(4*dt2)) * Z2 * step_fun (t,t0))
- ((1/(4*dt1)) * Z1 * step_fun (t,t1)) - ((1/(4*dt1)) * Z1 * step_fun(t,t0))
+ (((2*np.pi*cf)*Q1) * Zd * step_fun(t,t0) * inv_step_fun(t,t1))

# j=2
dZ2dt= ((1/(4*dt1)) * Z1 * step_fun (t,t1)) + ((1/(4*dt3)) * Z3 * step_fun (t,t1))
- ((1/(4*dt2)) * Z2 * step_fun (t,t2)) - ((1/(4*dt2)) * Z2 * step_fun(t,t1))
+ (((2*np.pi*cf)*Q2) * Zd * step_fun(t,t1) * inv_step_fun(t,t2))

# j=3
dZ3dt= ((1/(4*dt2)) * Z2 * step_fun (t,t2)) + ((1/(4*dt4)) * Z4 * step_fun (t,t2))
- ((1/(4*dt3)) * Z3 * step_fun (t,t3)) - ((1/(4*dt3)) * Z3 * step_fun(t,t2))
+ (((2*np.pi*cf)*Q3) * Zd * step_fun(t,t2) * inv_step_fun(t,t3))

# j=4
dZ4dt= ((1/(4*dt3)) * Z3 * step_fun (t,t3)) + ((1/(4*dt5)) * Z5 * step_fun (t,t3))
- ((1/(4*dt4)) * Z4 * step_fun (t,t4)) - ((1/(4*dt4)) * Z4 * step_fun(t,t3))
+ (((2*np.pi*cf)*Q4) * Zd * step_fun(t,t3) * inv_step_fun(t,t4))

# j=5
dZ5dt= ((1/(4*dt4)) * Z4 * step_fun (t,t4)) + ((1/(4*dt6)) * Z6 * step_fun (t,t4))
- ((1/(4*dt5)) * Z5 * step_fun (t,t5)) - ((1/(4*dt5)) * Z5 * step_fun(t,t4))
+ (((2*np.pi*cf)*Q5) * Zd * step_fun(t,t4) * inv_step_fun(t,t5))

# j=6
dZ6dt= ((1/(4*dt5)) * Z5 * step_fun (t,t5)) + ((1/(4*dt7)) * Z7 * step_fun (t,t5))
- ((1/(4*dt6)) * Z6 * step_fun (t,t6)) - ((1/(4*dt6)) * Z6 * step_fun(t,t5))
+ (((2*np.pi*cf)*Q6) * Zd * step_fun(t,t5) * inv_step_fun(t,t6))

# j=7
dZ7dt= ((1/(4*dt6)) * Z6 * step_fun (t,t6)) + ((1/(4*dt8)) * Z8 * step_fun (t,t6))
- ((1/(4*dt7)) * Z7 * step_fun (t,t7)) - ((1/(4*dt7)) * Z7 * step_fun(t,t6))
+ (((2*np.pi*cf)*Q7) * Zd * step_fun(t,t6) * inv_step_fun(t,t7))

# j=8
dZ8dt= ((1/(4*dt7)) * Z7 * step_fun (t,t7)) + ((1/(4*dt9)) * Z9 * step_fun (t,t7))
- ((1/(4*dt8)) * Z8 * step_fun (t,t8)) - ((1/(4*dt8)) * Z8 * step_fun(t,t7))
+ (((2*np.pi*cf)*Q8) * Zd * step_fun(t,t7) * inv_step_fun(t,t8))

return [dZ1dt,dZ2dt,dZ3dt,dZ4dt,dZ5dt,dZ6dt,dZ7dt,dZ8dt]

def step_fun(x,a):

if x==a:  return 0.5
elif x<a: return 0
elif x>a: return 1

def inv_step_fun(x,a):

if x==a:  return 0.5
elif x<a: return 1
elif x>a: return 0

####################################################################################################

import copy
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Number of equations:
N=8

# Define t:
t=np.arange(0,50,0.025)
t_copy=copy.deepcopy(t)
Nt=len(t)

# Define other parameters/arguments:
Q=(0.008017307188069223, 0.008117153569161557, 0.008081825341069762, 0.00803759248462624, 0.00803759248462624, 0.008081825341069762, 0.008117153569161557, 0.008017307188069223)
Z0=np.array([22050, 0., 0., 0., 0., 0., 0., 0.])
Zd0=Z0[0]
dts=np.array([20.67120623,  1.43225437,  1.55738981,  1.63132137,  1.63132137, 1.55738981,  1.43225437, 20.67120623])
ts=np.array([20.67120623, 22.10346059, 23.66085041, 25.29217178, 26.92349315, 28.48088296, 29.91313733, 50.58434356])
cf=1.5

# Define Zd: this is a vector of size Nt, but I will need to give the function only the right sample
# at each time t.
Zd=1/(1+np.exp(t))*5e4

# Solve system of differential equations:
Z=odeint(ODE_solver,Z0,t,args=(Q,Zd,Zd0,dts,ts,cf,t_copy))

# Let's rearrange this results so it is easier to read in the future:
Zs=[]
for n in range(N):
Zs.append([])
for vec in Z:
Zs[n].append(vec[n])
Zs[n]=np.array(Zs[n])

plt.figure()
for vec in Zs:
plt.plot(t,vec)


• Could you describe the terms of dZ2dt in their physical meaning? From the flow of arguments I would expect that the 4th term has an inv, that is, Z2 is on for t<t1 and for t>t2 and thus off for t1<t<t2, and where it is off it gets replaced by the Zd contribution. Do you want Zd to be a finely grained step function? Why not directly use the formula to compute its value at t? Zd_func=lambda t: 1/(1+np.exp(t))*5e4 and then in the ODE function use Zd=Zd_func(t). – Lutz Lehmann Feb 12 '19 at 10:44
• As I said before, the Z functions represent the evolution of energy over time for each one of these systems. For t1<t, the "input energy" Zd is in the first system, but it moves forward, so that at t2, it enters the second system, at t3 the third one and so on. The fourth term of these equations is on only for the time interval for which this "input" energy remains within the corresponding system. – Xay Feb 12 '19 at 13:25
• The rest of the terms represent the leakage of energy from each system into the systems above/below, so they are turned on/off at different times (for t<t1, there is energy leaking from system 1 into system 2, and for t>t2, there is energy flowing from system 3 into system 2). I hope this makes it easier to understand. – Xay Feb 12 '19 at 13:27
• Why do you say Zd is a finely grained step function? In the example I posted here, Zd is calculated from a function, but in reality this Zd comes from data and can't be calculated with a function. I do have, though, a Zd sample for each t sample and it is, as the Zd I have here, a function that decays exponentially over time. – Xay Feb 12 '19 at 13:39
• The internal points that the ODE function is called are not the same ones as the points given to evaluate for the output. The internal points can be more or less dense than the output times, and usually do not coincide. The output is interpolated from the internal steps. In that view, it would probably be better to just use interpolation, Zd=np.interp(t,t_samples, Zd_samples). It would probably also a good idea to pass a maximal step size smaller than the compartment sizes. – Lutz Lehmann Feb 12 '19 at 14:00

For the interpretation where the left and right "virtual" compartments are infinite sinks at level Z0=Z9=0, and the walls between the compartments j and j+1 are removed at time ts[j] where also the external filling of compartment j starts, I get the shortened implementation

def ODE_func(Z,t,Q,Zd,Z0,dts,ts,cf,t_copy):

# Declare some boundary conditions:
Z9=0; dt9=1/4; t9=ts[-1]+dt9 # (It doesn't matter which value I have here, since Z9 = 0)
Z0=0; dt0=1/4; t0=0;
# Define functions to solve for:
Z = np.concatenate([[Z0],Z,[Z9]]);
# Define time intervals:
dt = np.concatenate([[dt0],dts,[dt9]])
# Define cumulative times:
ts = np.concatenate([[t0],ts,[t9]]);

# Find index of the time sample that is being used at the moment:
Zd=np.interp(t,t_copy,Zd);

dZdt = [ (Z[j-1]/dt[j-1]-Z[j]/dt[j])*step_fun(t,ts[j-1])/4
+ (Z[j+1]/dt[j+1]-Z[j]/dt[j])*step_fun(t,ts[j])/4
+ 2*np.pi*cf*Q[j-1]*Zd*step_fun(t,ts[j-1])*step_fun(ts[j],t) for j in range(1,9)]

return dZdt


and the plots of the "fill heights" Z[j]/dt[j], j=1,...,8 over the long term as

where indeed the levels start to fall after some time.

• This is exactly the kind of solution I expected to see, but my time vector is a lot shorter, it only goes up to ~50 seconds. Could you please elaborate on what you did there? – Xay Feb 12 '19 at 15:45
• Then you have to mentally cut the images at t=50. I just extended the computation by extending the t=linspace(0,250,0.25), using that 1/(1+exp(x))=(1-1/(1+exp(-x))) to avoid overflow in the exponential. In externally provided data, extend the t_copy and Zd arrays by 250 and 0 to have values available in the interpolation. – Lutz Lehmann Feb 12 '19 at 15:51
• Thanks a lot, but after about 50 second my "input" energy leaves the last system, so at this time the energy should not be growing. It doesn't make much sense, physically. – Xay Feb 12 '19 at 15:54
• That is also true numerically, as 1/(1+exp(50)) is essentially zero. The last compartments are rising because the first compartments have a much higher "height" value. So while the first compartments starts to fall early, it does so from a height of 1000, overpowering the compartments after it that have drastically lower heights, 500, 350, 250,... – Lutz Lehmann Feb 12 '19 at 16:00
• I have tried making my initial conditions equal to zero for all systems and changing the initial amplitude of my Zd, but I still get an increase in the Zj functions (the amplitudes are different, but they still grow). I have also tried solving it with ode, in case it made a difference, but the solutions are identical. – Xay Feb 12 '19 at 16:19