0

Let's say that I have a system of differential equations and I want to solve it with odeint. Some of the system's constants are time depended and I have their values stored in arrays (a,b,c and d with shape (8000, ) ). I want the system to use a different value of these constants at each time step. See the simplified code example:

t = np.linspace(0,100,8000)

a = np.array([0, 5, 6, 12, 1.254, ..., 0.145])     # shape (8000, )
b = np.array([1.45, 5.9125, 1.367, ..., 3.1458])
c = np.array([0.124, 0.258, 0.369, ..., 0.147])
d = np.array([7.145, 5.123, 6.321, ..., 0.125])

def system(k,t):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)

    return [vcx_i_dot, vcy_i_dot, psi_i_dot wz_i_dot]

k0 = [0.1257, 0, 0, 0]

k = odeint(system, k0, t)

vcx_i = k[:,0]
vcy_i = k[:,1]
psi_i = k[:,2]
wz_i = k[:,3]

psi_i = [system(t_i, k_i)[2] for k_i, t_i in zip(k,t)]
wz_i = [system(t_i, k_i)[3] for k_i, t_i in zip(k,t)]

The most relevant solution I found so far is this: Solving a system of odes (with changing constant!) using scipy.integrate.odeint? But since I have only values of my variables in arrays and not the equation of the variables that depend on time (e.g. a=f(t)), I tried to aply an interpolation between the values in my arrays, as shown here: ODEINT with multiple parameters (time-dependent) I managed to make the code running without errors, but the total time increased dramatically and the results of the system solved are completely wrong. I tried any possible type of interpolation that I found here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html but still wrong outcome. That means that my interpolation isn't the best possible, or my points in the arrays (8000 values) are too much to interpolate between them and solve the system correctly. What is the right way to solve a system like this? I am new to python and I use python 2.7.12 on Ubuntu 16.04 LTS. Thank you in advance.

Meanmachine
  • 97
  • 12
  • *"... and the results of the system solved are completely wrong."* How do you know the results are wrong? Also, if you don't show us the code that you tried, we can't help figure out what the problem might be. It is much easier for someone to help you if you provde a [minimal, complete and verifiable example](https://stackoverflow.com/help/mcve). – Warren Weckesser Feb 21 '19 at 22:48
  • @WarrenWeckesser I know how the correct results look like because I have solved a bigger system, which includes the system that I want to solve now. The results must be the same. I posted my full code here: https://stackoverflow.com/questions/54773543/use-numpy-arrays-as-arguments-in-odeint but the post was a little bit chaotic, so I tried to simplify the problem in the post above. – Meanmachine Feb 22 '19 at 11:21

1 Answers1

1

The interpolators are usually very fast, so probably there is something else in your function. You can however, attempt different interpolators (like InterpolatedUnivariateSpline), or reducing the interpolation nodes to increase speed. But I would aim at your integration instead.

Lately, ode and odeint have been replaced by other, more flexible functions (see here)

I would start with and explicit method instead of an implicit one (default for solve_ivp is runge kutta, while default for odeint is LSODA):

interp = scipy.interpolate.interp1d(t,(a,b,c,d))

def system(t,k):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]
    a, b, c, d = interp(t)

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)
    return [vcx_i_dot, vcy_i_dot, psi_i_dot, wz_i_dot]

k0 = [0.1257, 0, 0, 0]
steps = 1
method = 'RK23'
atol = 1e-3
s = solve_ivp(dydt, (0, 100), k0, method=method, t_eval=t, atol=atol, vectorized=True)

you can increase / reduce atol or change the method.

Tarifazo
  • 4,118
  • 1
  • 9
  • 22
  • I tried to implement the solve_ivp solution as you mentioned and I got the error: 'vcx_i = k[0] TypeError: 'float' object has no attribute '__getitem__' ' – Meanmachine Feb 22 '19 at 12:16
  • Reverse the order of t and k in `system`(see edit - my mistake in reading function docs, also different from `odeint`) – Tarifazo Feb 22 '19 at 12:22
  • The system looks to work without errors, but in order to see the results I tried to save the output, like I did with odeint (see updated code above) and I get the error ` vcx_i = k[:,0] TypeError: unhashable type` Is there any difference between odeint and solve_ivp in saving the results of the system? – Meanmachine Feb 22 '19 at 13:26
  • Yes. `odeint` only returns the results (`y` vector), `solve_ivp` returns an instance of the `OdeResult` object which also stores other things (e.g. initial and final time, steps, etc.). You want to use `s.y` to get the result vector. – Tarifazo Feb 22 '19 at 13:33
  • So far so good. I managed to print the results. Now, let's say that I want to repeat a part of the same system (see updated code above), just like you pointed out in : https://stackoverflow.com/questions/54365358/odeint-function-from-scipy-integrate-gives-wrong-result with odeint. The general problem is the same, but this time I am using solve_ivp. How can I do that? – Meanmachine Feb 22 '19 at 14:13
  • ... when I do it like above I get this error: vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i) TypeError: ufunc 'cos' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe'' – Meanmachine Feb 22 '19 at 14:27
  • in order to make it more clear, I want to solve the half of the system I had in [https://stackoverflow.com/questions/54365358/odeint-function-from-scipy-integrate-gives-wrong-result] (2nd half), so I put the results of the 1rst half in arrays and use them as inputs in the 2nd half (see this post). When I use solve_ivp, I have the same problem I had in the previous post (I get wrong theta_i). Is there any similar solution to the one you pointed out in the previous post, but with the use of solve_ivp ? – Meanmachine Feb 23 '19 at 12:31
  • I solved the problem of the previous comments in the post [https://stackoverflow.com/questions/54843250/converting-odeint-system-to-solve-ivp-dimensions-problem] – Meanmachine Feb 26 '19 at 12:56