Solve systems of ODEs in Julia with initial condition as 2D array

user_8583532 0
Viewed 905

http://diffeq.sciml.ai/latest/tutorials/ode_example.html

I am trying to use the ODE solver in Julia (DifferentialEquations.jl) to solve a system of n interacting particles. Let's say that the system is in 2D, and the equation of motion of each particle is described by a second-order ODE of its position with respect to time. Then four variables are needed for each particle, two for positions and two for velocities. Then 4n variables are needed to be declared. Is there a way to generalize, such that one does not need to list all 4n equations one by one?

For example:

http://diffeq.sciml.ai/latest/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1

I try to modify the Lorenz equation in the link above a little bit to n particles (which is a very very rough attempt since I actually have no idea how to do it) by trying to extend u and du to 2D arrays.

using DifferentialEquations
using Plots
n = 4
function lorenz(du,u,p,t,i)
    du[i,1] = 10.0*(u[i,2]-u[i,1])*sum(u[1:n,1])
    du[i,2] = (u[i,1]*(28.0-u[i,3]) - u[i,2])*sum(u[1:n,1])
    du[i,3] = (u[i,1]*u[i,2] - (8/3)*u[i,3])*sum(u[1:n,1])
end

u0 = hcat([1.0;0.0;0.0], [0.0;1.0;0.0], [0.0;0.0;1.0])
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob)

This, without surprise, does not work, but I hope that you get the idea what I am trying to do. Is there anyway for the ODE solver to solve u as a 2D array (or other ways that can achieve similar purposes?)

differential-equations julia arrays
1 Answers

The problem is not the 2D Array. For example replacing your lorenz definition with 

function lorenz(du,u,p,t)
    du[1,1] = 10.0*(u[1,2]-u[1,1])
    du[1,2] = (u[1,1]*(28.0-u[1,3]) - u[1,2])
    du[1,3] = (u[1,1]*u[1,2] - (8/3)*u[1,3])
end

will work.

The problem is the function signature, the additional i is not supported. If you want to solve a network of Lorenz oscillators you have to code it with a function with the same signature, e.g. lorenz_network!(du, u, p, t) for the inplace version. Put a loop over the individual oscillators in your function and you are almost there.

user_47170970
Related