# Vectorised second order ode solving in python

I am trying to write a python program that simulates the motion of a large number of particles by numerically integrating a second order ordinary differential equation. I first split the ODE into two coupled first order ODEs and solve using scipy.integrate.solve_ivp following a similar method to what is described in the answer to this question. However, my problem is that I want to solve this system for a large number of particles, each with different initial conditions. Naively I could do this with a for loop, but I'm sure numpy and scipy must have a way of vectorising this operation that would be much faster.

I have had a look at the documentation for scipy.integrate.solve_ivp and it talks about vectorisation, but not in the way I want. I would expect for a system with two coupled first order ODEs and n particles you could input the initial conditions as an array with size (2,n) but this is not the case.

Is there a way to solve for multiple initial conditions without resorting to a slow python for loop?

For reference the system of ODEs I want to solve looks like

dx/dt = v
dv/dt = F(x,v)


With initial conditions in an array like

initialConditions = [[x0,v0],[x1,v1],...,[xN,vN]]

• In order to get a compatible input for the scipy function, you need to write a large ODE system $\pmb{\dot{y}} = \pmb{f}(\pmb{y})$ with 2N equations for the 2N unknowns $\pmb{y} := (x_1,x_2,\dots,x_N,v_1,v_2,\dots,v_N)^{\top}$. You may use a different ordering of the unknowns if you like. – Christoph Apr 26 '20 at 8:19
• @Christoph thanks for that. As I expected, solving the ODE system using 2N equations as you suggested does provide some speed advantage, but only for relatively small N. My code was about five times faster using 2N equations than looping N times, up to about N=100. For larger N, the advantage became less and less, so for a very large N parallelization seems to be the way to go. – spacenut1 May 4 '20 at 5:45
• You're welcome. Parallelization is easy as long as the particles do not interact with each other. Otherwise, it becomes more difficult. – Christoph May 4 '20 at 7:09

Is there a way to solve for multiple initial conditions without resorting to a slow python for loop?
I think the solution to your problem might be parallelization, not vectorization, unfortunately. While Python is going to be slower than C or C++, the for loop in question isn't actually a part of the computation, it's a different problem numerically because of the different initial conditions. Besides, there's so much data involved in the computation there wouldn't be any spare room on the cache lines anyhow.