Генерация функции в цикле для дифференциальных уравнений в Юлии

Question

Принимая пример Лоренца ( JuliaDiffEq GitHub)

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

Если бы мне пришлось генерировать уравнения в цикле, я бы попытался объединить строки в цикле, а затем преобразовать их в функцию. С Matlab это работает ( str2func). С Юлией я потерялся (это на самом деле самое первое, что я пытаюсь в Юлии). Моя первая попытка была:

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

Не работает (u не определено; за ним следуют строки ошибок, вероятно, возникшие в результате использования этой функции в prob = ODEProblem(lorenz,u0,tspan) Может ли кто-нибудь, пожалуйста, указать мне в правильном направлении. Первая проблема заключается в том, что я не смог найти, что происходит с функцией, когда она входит в ODEProblem, Кто-то, пожалуйста, направьте меня туда, где я могу узнать, что на самом деле есть функция в Юлии.

---

EDIT1

Я только что попробовал

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

Это работает. Так что я знаю, как это работает. Но я боюсь, что мог бы наткнуться на "путь"; не предпочтительный способ. Может ли кто-то более знаком с Julia Прокомментируйте, пожалуйста.

---

EDIT2

Пример системы уравнений, которую необходимо решить:

du[1] = u[1] + v[1] + U
dv[1] = v[1] + U

du[2] = u[2] + v[2] + U
dv[2] = v[2] + U

...

Вот U это сумма u[i] (я иду с 1 до M; должно быть возможно установить значение M в переменной).

Answer 1

В MATLAB вы должны попытаться избежать цикла. В Джулии циклы в таком контексте оптимизированы, поэтому нет причин избегать этого. Простой способ написать эту функцию:

function f(t,u,du)
  U = sum(u)
  for i in eachindex(u)
    du[i] = u[i] + v[i] + U
    dv[i] = v[i] + U
  end
end

Обратите внимание, что для этого даже не нужна ссылка на длину, поэтому он будет работать с любым размером u а также du, Так создай ODEProblem(f,u0,tspan) где u0 ваш размер M массив и тебе будет хорошо.

Мы можем использовать некоторые макросы производительности, чтобы оптимизировать это немного больше. @inbounds отключает проверку границ, так что это может быть приятно сделать

function f(t,u,du)
  U = sum(u)
  @inbounds for i in eachindex(u)
    du[i] = u[i] + v[i] + U
    dv[i] = v[i] + U
  end
end

В некоторых случаях, если цикл достаточно длинный, мы можем захотеть @simd в теме:

function f(t,u,du)
  U = sum(u)
  @simd for i in eachindex(u)
    @inbounds du[i] = u[i] + v[i] + U
    @inbounds dv[i] = v[i] + U
  end
end

или, если он действительно длинный, многопоточность:

function f(t,u,du)
  U = sum(u)
  Threads.@threads for i in eachindex(u)
    @inbounds du[i] = u[i] + v[i] + U
    @inbounds dv[i] = v[i] + U
  end
end

(и убедитесь, что вы включаете многопоточность, конечно).