1

我相信这段代码有一个错误。为简洁起见，我将只编写定义 ODE 的函数

function clones(du,u,p,t)

    (Nmut,f) = p

    # average fitness
    phi = sum(f.*u)

    # constructing mutation kernel
    eps = 0.01
    Q = Qmatrix(Nmut,eps)

    # defining differential equations
    Nclones=2^Nmut;
    du = zeros(Nclones)

    ufQ = transpose(transpose(u.*f)*Q)
    du = ufQ .- phi*u

end

如果需要整个代码，我可以提供它，但它很乱，我不确定如何创建一个最小的示例。我在 Nmut = 2 时尝试过这个，所以我可以与硬编码版本进行比较。du 在第一个时间步的输出是相同的。但是这个版本似乎永远不会更新你，它保持在规定的 u0 上。

有谁知道为什么会这样？我也可以提供完整的脚本，但如果有人能看到你为什么不更新，我想避免这种情况。

编辑：

maxdim=4;
for i in 1:maxdim
    du[i] = 0.0;
    for j in 1:maxdim
        du[i] += u[j].*w[j].*Q[j,i] 
    end
    du[i] -= u[i].*phi
end

如果我们使用这个版本，du 会正确更新。为什么会这样？

4

1 回答 1

1

您正在使用就地表单。有了这个，你必须更新du. 在您使用的脚本中

du = ufQ .- phi*u

那就是用新数组替换名称du，但不会改变原始du数组中的值。一个快速的解决方法是使用变异等于：

du .= ufQ .- phi*u

请注意，这是.=.

要了解这在更多基于示例的格式中意味着什么，请考虑这一点。我们有一个数组：

a = [1,2,3,4]

现在我们将一个新变量指向同一个数组

a2 = a

当我们更改 的值时，a我们会看到反映的值，a2因为它们指向相同的内存：

a[1] = 5
println(a2) # [5,2,3,4]

但是现在如果我们替换a我们没有注意到任何事情发生，a2因为它们不再引用同一个数组

a = [1,2,3,4]
println(a2) # [5,2,3,4]

像DifferentialEquations.jl 这样的包使用了变异形式，这样用户就可以通过重复更改缓存数组中的值来摆脱数组分配。因此，这意味着您应该更新du而不是替换其指针的值。

如果您觉得不使用突变更舒服，您可以使用函数 syntax f(u,p,t)，尽管如果状态变量是（非静态）数组，这样做会影响性能。

于 2018-09-30T02:25:31.167 回答