Julia BVP 求解器中边界条件代码的令人费解的结果

Puzzling result from boundary condition code in Julia BVP solver

我正在尝试使用 BoundaryValueDiffEq 包解决 Julia 中的边界值问题,按照发现的示例 here。在边界条件函数中,示例需要一个for循环来单独更新每个索引,à la

function bc1!(residual, u, p, t)
    for i in 1:n
        residual[i] = u[end][i] - 10
    end
end

我想用下面的代码,应该效率更高:

function bc1!(residual, u, p, t)
    residual = u[end] .- 10
end

虽然两个版本的代码的 residual 的结果值相同,但求解器在第一种情况下给出了正确的结果,而在第二种情况下给出了错误的结果。

我能想到的就是更新residual有一些区别 逐个索引并为其分配一个新向量,即使结果在值和类型上都相同。为什么会这样,有没有可能在保持正确结果的同时让代码更高效?

这是完整的代码,希望对您有帮助。

using BoundaryValueDiffEq, Plots

n = 3
f(t) = .1
F(t) = .1*t

function du!(du,u,p,t)
    fn(i) = 1/(u[i]-t)
    for i in 1:n
        du[i] = 1/(n-1)*F(u[i])/f(u[i])*((2-n)/(u[i]-t)+sum(map(fn,
        vcat(1:i-1,i+1:n))))
    end
end

function bc1!(residual, u, p, t)
    #residual = u[end] .- 10
    for i in 1:n
        residual[i] = u[end][i]-10
    end
end

# exact solution
xvals = LinRange(0,20/3,200)
yvals = 1.5*xvals

# solving BVP
tspan = (0.0,20/3)
bvp1 = BVProblem(du!, bc1!, 10*ones(Int8,n), tspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt=.2)

# plotting computed solution vs actual solution
plot(sol1,vars=(0,1))
plot!(xvals,yvals,label="Exact solution")

您覆盖了数组而不是改变它。您需要使用 .= 来更新它 in-place。

function bc1!(residual, u, p, t)
    residual .= u[end] .- 10
end

或更安全:

function bc1!(residual, u, p, t)
    @. residual = u[end] .- 10
end