I would like to solve a 1D equation using NonlinearSolve.
Appended below is a toy MNWE followed by the error message generated.
The issue I’m encountering is that in the actual problem I need to pass a workspace array ws_array
for the calculation. After reading the documentation and various threads on Julia Discourse, ComponentArrays
seemed to be the way to proceed.
# test_NonlinearSolve.jl
using ComponentArrays
using NonlinearSolve
@inline function a_function(
t::Float64,
cp_v::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(q = 1, p = 2, g = 3, l = 4, m = 5, w = ViewAxis(6:21, Shaped1DAxis((16,))))}}}
)::Float64
@static_unpack q, p, g, l, m = cp_v
cp_v.w[1] = sin(2.0*q)
cp_v.w[2] = l^(-1)
cp_v.w[3] = sin(q)
cp_v.w[4] = m^(-1)
cp_v.w[5] = cos(q)
cp_v.w[6] = cp_v.w[5]^2
cp_v.w[7] = cp_v.w[3]^2
cp_v.w[6] = - (1.0/2.0)*cp_v.w[6] + cp_v.w[7]
cp_v.w[7] = g^2
cp_v.w[6] = cp_v.w[4]*cp_v.w[6]*cp_v.w[7]*p^2
cp_v.w[8] = (1.0/2.0)*cp_v.w[5]
cp_v.w[9] = cp_v.w[2]^3
cp_v.w[8] = - g*p^4*cp_v.w[4]^3*cp_v.w[9]*cp_v.w[8]
cp_v.w[6] = cp_v.w[6] + cp_v.w[8]
cp_v.w[6] = cp_v.w[6]*cp_v.w[9]
cp_v.w[8] = cp_v.w[3]*cp_v.w[5]
cp_v.w[9] = - (1.0/3.0)*cp_v.w[1] + (1.0/2.0)*cp_v.w[8]
cp_v.w[9] = m*cp_v.w[3]*cp_v.w[9]*g^3
cp_v.w[6] = cp_v.w[9] + (1.0/3.0)*cp_v.w[6]
cp_v.w[6] = t*cp_v.w[2]*cp_v.w[6]
cp_v.w[8] = (1.0/2.0)*cp_v.w[1] - cp_v.w[8]
cp_v.w[7] = p*cp_v.w[8]*cp_v.w[7]*cp_v.w[2]^2
cp_v.w[6] = (1.0/2.0)*cp_v.w[7] + cp_v.w[6]
f_t = t^2*cp_v.w[6]
return f_t
end
function main()
# Workspace array
ws_array::Array{Float64, 1} = zeros(Float64, 16)
# Initial value of scalar t
t_o::Float64 = 0.5
cmpnt_vec =
ComponentVector(
q = 0.7552444383914849,
p = -0.31175068141804535,
g = 6.67408, l = 1.0, m = 0.1,
w = ws_array
)
@show typeof(cmpnt_vec)
println("")
# Test a_function
@show f_t = a_function(t_o, cmpnt_vec)
println("")
# Test NonlinearSolve
nl_prblm = NonlinearProblem(a_function, t_o, cmpnt_vec)
@show t_f = solve(nl_prblm, SimpleNewtonRaphson(), abstol = 1.0e-12)
println("")
end
begin
main()
end
However, this code generates the following error message:
julia> include("test_NonlinearSolve.jl")
. . .
ERROR: LoadError: MethodError: no method matching a_function(::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, ::ComponentVector{Float64, Vector{…}, Tuple{…}})
The function `a_function` exists, but no method is defined for this combination of argument types.
Closest candidates are:
a_function(::Float64, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(q = 1, p = 2, g = 3, l = 4, m = 5, w = ViewAxis(6:21, Shaped1DAxis((16,))))}}})
@ Main ~/projects/Hamiltonian/test/test_NonlinearSolve.jl:6
a_function(::Float64, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(q = 1, p = 2, g = 3, l = 4, m = 5, w = ViewAxis(6:69, Shaped1DAxis((16,))))}}})
@ Main ~/projects/Hamiltonian/test/test_NonlinearSolve.jl:6
Stacktrace:
[1] (::NonlinearFunction{…})(::ForwardDiff.Dual{…}, ::Vararg{…})
@ SciMLBase ~/.julia/packages/SciMLBase/dld4W/src/scimlfunctions.jl:2604
[2] pushforward(f::NonlinearFunction{…}, backend::AutoForwardDiff{…}, x::Float64, tx::Tuple{…}, contexts::DifferentiationInterface.Constant{…})
@ DifferentiationInterfaceForwardDiffExt ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl:40
[3] derivative(f::NonlinearFunction{…}, backend::AutoForwardDiff{…}, x::Float64, contexts::DifferentiationInterface.Constant{…})
@ DifferentiationInterfaceForwardDiffExt ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl:206
[4] compute_jacobian!!(::Nothing, prob::SciMLBase.ImmutableNonlinearProblem{…}, autodiff::AutoForwardDiff{…}, fx::Float64, x::Float64, ::SimpleNonlinearSolve.Utils.DINoPreparation)
@ SimpleNonlinearSolve.Utils ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/utils.jl:110
[5] #__solve#32
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/raphson.jl:55 [inlined]
[6] __solve
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/raphson.jl:34 [inlined]
[7] #simplenonlinearsolve_solve_up#40
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/SimpleNonlinearSolve.jl:123 [inlined]
[8] simplenonlinearsolve_solve_up
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/SimpleNonlinearSolve.jl:118 [inlined]
[9] #solve#39
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/SimpleNonlinearSolve.jl:109 [inlined]
[10] solve
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/SimpleNonlinearSolve.jl:92 [inlined]
[11] #solve#36
@ ~/.julia/packages/SimpleNonlinearSolve/GiCJO/src/SimpleNonlinearSolve.jl:64 [inlined]
which, if I understand correctly, ComponentArrays
is not compatible with the default ForwardDiff
.
How to solve this issue? Should I use SciMLStructs and, if so, then how?