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 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]]