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
我正在尝试使用 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