使用Julia的NonlinearSolve求解稳态时函数调用错误

我正在用NonlinearSolve求解稳态问题，代码如下：

prob = SteadyStateProblem(eoms!, σ₀, p)
sol = solve(
    prob,
    SSRootfind(),
    abstol = slv.abstol,
    reltol = slv.reltol,
)

eoms!函数定义如下：

function eoms!(
    du::Vector{Float64},
    u::Vector{Float64},
    p::Tuple{Any,Any,Any,Any,Any},
    t::Float64 = 0.0,
)

第一次调用eoms!时参数类型正确，但第二次调用时，du的类型变成了ForwardDiff.Dual相关类型（具体类型如下）：

ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}, Float64, 8}[Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#41#50"{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.Solve.eoms!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{Main.Solve.Instance, Main.SymmetricLattice.Instance, Main.Laser.Instance, Float64, Int64}}, Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)]

注：用ODEProblem做时间演化时，求解器能得到合理结果。

问题出在eoms!函数的类型注解上：你把du和u的类型固定成了Vector{Float64}，但稳态求解器SSRootfind()会自动用自动微分（ForwardDiff）计算雅可比矩阵，这时候会传入ForwardDiff.Dual类型的数组，和你注解的类型不匹配，导致报错。

而ODE求解器默认可能没启用自动微分，或者用了其他方式，所以没触发这个问题。

修复方法很简单：去掉du、u和t的具体类型注解，让函数支持任意数值类型：

function eoms!(
    du::AbstractVector,
    u::AbstractVector,
    p::Tuple{Any,Any,Any,Any,Any},
    t = 0.0,
)
    # 你的函数实现代码
end

如果需要保留一定的类型约束，可以用参数化类型：

function eoms!(
    du::AbstractVector{T},
    u::AbstractVector{T},
    p::Tuple{Any,Any,Any,Any,Any},
    t::T = zero(T),
) where T<:Real
    # 你的函数实现代码
end

这样函数就能兼容Float64和ForwardDiff.Dual类型，稳态求解器就能正常运行了。

内容的提问来源于stack exchange，提问作者mrW