1

以洛伦兹为例（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）。对于 Julia，我迷路了（这实际上是我在 Julia 中尝试的第一件事）。我的第一次尝试是：

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

不起作用（你未定义；随后出现的错误行可能是由于使用此函数prob = ODEProblem(lorenz,u0,tspan) 引起的进入ODEProblem. 有人还请指导我到哪里可以找到 Julia 中的真正功能。

编辑1

我刚试过

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评论一下。

编辑2

需要求解的方程组示例：

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在变量中设置 的值）。

4

1 回答 1

1

在 MATLAB 中，您应该尽量避免循环。在 Julia 中，这种上下文中的循环已经过优化，因此没有理由避免它。编写此函数的一种简单方法是：

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)在哪里，你会很好。u0M

我们可以使用一些性能宏来进一步优化它。@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

（并确保您当然启用线程）。

于 2017-12-29T14:35:00.403 回答