Решить системы ОДУ в Юлии с начальным условием в виде двумерного массива - PullRequest
Новая
       21

Решить системы ОДУ в Юлии с начальным условием в виде двумерного массива

1 голос
/ 14 октября 2019

http://docs.juliadiffeq.org/latest/tutorials/ode_example.html

Я пытаюсь использовать решатель ODE в Юлии (DifferentialEquations.jl) для решения системы из n взаимодействующих частиц. Скажем, система находится в 2D, а уравнение движения каждой частицы описывается ОДУ второго порядка ее положения по времени. Тогда для каждой частицы необходимы четыре переменные, две для положения и две для скорости. Затем необходимо объявить 4n переменных. Есть ли способ обобщения, чтобы не нужно было перечислять все 4n уравнений одно за другим?

Например:

http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1

Я пытаюсь изменитьуравнение Лоренца в приведенной выше ссылке немного на n частиц (что является очень и очень грубой попыткой, поскольку я на самом деле не знаю, как это сделать), пытаясь расширить u и du до двумерных массивов. 

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)

Это, без удивления, не работает, но я надеюсь, что вы поняли, что я пытаюсь сделать. В любом случае решатель ODE может решить u как двумерный массив (или другими способами, которые могут достичь аналогичных целей?)

1 Ответ

3 голосов
/ 14 октября 2019

Проблема не в 2D-массиве. Например, замена определения lorenz на 

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

будет работать.

Проблема в сигнатуре функции, дополнительная i не поддерживается. Если вы хотите решить сеть осцилляторов Лоренца, вы должны закодировать ее с помощью функции с такой же сигнатурой, например, lorenz_network!(du, u, p, t) для версии на месте. Поместите петлю над отдельными генераторами в вашей функции, и вы почти там.

...