# NonlinearSolve: how to pass parameters including a mutable workspace array

**URL:** <https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [August 6, 2025, 6:45pm UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421 "2025-08-06T18:45:01Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Audrius-St](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/audrius-st/32/24175_2.png) [@Audrius-St](https://discourse.julialang.org/u/Audrius-St)\
**Post date:** [August 6, 2025, 6:45pm UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421/1 "2025-08-06T18:45:01Z")

</div>

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.

```julia-auto
# 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-auto
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?

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [August 7, 2025, 6:24am UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421/2 "2025-08-07T06:24:11Z")

</div>

Hi 🙂 First of all, thanks for providing a minimal example and the full error message.

Regarding your issue: Try removing all type annotations from the function signature, i.e. turn this

```julia
@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

```

into this

```julia
function a_function(t, cp_v)

```

What happens here is not that `ComponentArrays` is incompatible with `ForwardDiff`, but rather `a_function` is incompatible. `ForwardDiff` passes `Dual` numbers through the functions it wants to differentiate (see e.g. [here](https://juliadiff.org/ForwardDiff.jl/stable/dev/how_it_works/) for an explanation). If you restrict them to only take `Float64` as input you will get the error message you posted (note the type of the first argument is not `Float64`, that’s the main issue

```julia-auto
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

```

The second issue that will appear after fixing the signature is that the `ComponentVector` you use as a “buffer” for storing intermediate values is also fixed to be of type `ComponentVector{Float64, …}` and `ws_array` is a `Vector{Float64}`. So as soon as you try to do something like `cp_v.w[…] = t * …` it won’t work, because `t` will be a `Dual` and not a `Float64` (for `ForwardDiff` to work). But a `Vector{Float64}` can only accept `Float64` values. The solution is to either not use the array at all (you don’t need it, as far as I can see) or to use a “diff cache” which contains/creates arrays of all the different types you might need.

[EDIT: Of course, if you really need the array for some reason the following doesn’t really apply. But perhaps then you could clarify why you need it – there might still be ways around using the array.] Why you don’t need a component array is because you just compute a single number from a fixed set of inputs and operations. In `a_function` you are basically maintaining your own “function stack” by re-using all the entries of the vector `w`. This is quite confusing and most likely slower than just using the real function stack (i.e. just make regular variables inside the function)… Just use a `NamedTuple` for all the parameters and get rid of the vector all together (you can search\replace all `cp_v[i]` by some variable `w_i` if you want, or, ideally, try to make it more transparent what happens. Something like this:

```julia
# test_NonlinearSolve.jl

using NonlinearSolve

function b_function(
			t,
			parameters
		)

	(; q, p, g, l, m) = parameters
	w_1 = sin(2.0*q)
	w_2 = l^(-1)
	w_3 = sin(q)
	w_4 = m^(-1)
	w_5 = cos(q)
	w_6 = w_5^2
	w_7 = w_3^2
	w_6 = - (1.0/2.0)*w_6 + w_7
	w_7 = g^2
	w_6 = w_4*w_6*w_7*p^2
	w_8 = (1.0/2.0)*w_5
	w_9 = w_2^3
	w_8 = - g*p^4*w_4^3*w_9*w_8
	w_6 = w_6 + w_8
	w_6 = w_6*w_9
	w_8 = w_3*w_5
	w_9 = - (1.0/3.0)*w_1 + (1.0/2.0)*w_8
	w_9 = m*w_3*w_9*g^3
	w_6 = w_9 + (1.0/3.0)*w_6
	w_6 = t*w_2*w_6
	w_8 = (1.0/2.0)*w_1 - w_8
	w_7 = p*w_8*w_7*w_2^2
	w_6 = (1.0/2.0)*w_7 + w_6
	f_t = t^2*w_6
	return f_t
end

function main()
    t_o = 0.5

    parameters = (;
            q = 0.7552444383914849, 
            p = -0.31175068141804535, 
            g = 6.67408, l = 1.0, m = 0.1, 
        )

    # Test b_function
    @show f_t = b_function(t_o, parameters)
    println("")

    # Test NonlinearSolve
    nl_prblm = NonlinearProblem(b_function, t_o, parameters)
    @show t_f = solve(nl_prblm, SimpleNewtonRaphson(), abstol = 1.0e-12)
    println("")
end

begin
    main()
end

```

Note that removing the type annotations in the function signature has no impact on performance, but makes all the difference for writing generic code (working with regular floats and also duals). Unless you need dispatch (multiple methods which do different things based on the input type) type annotations of function arguments are usually unnecessary and often harmful when composing code (like here).

Also note that I removed `@inline`. No idea if it will have an impact here, but I doubt it (in my experience it was rarely necessary to speed things up, but you can test the speed of the solve with and without it if you are unsure).

---

<div class="post-metadata">

**Author:** ![Audrius-St](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/audrius-st/32/24175_2.png) [@Audrius-St](https://discourse.julialang.org/u/Audrius-St)\
**Post date:** [August 7, 2025, 7:07pm UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421/3 "2025-08-07T19:07:36Z")

</div>

Hello. Thank for your detailed explanation and example solution.  
Much appreciated.

I was wondering why the type annotations were not specified for functions in the `NonlinearSolve` documentation, now I understand that it’s due to the passing of `Dual` numbers as arguments.

Also, you’re quite right that I don’t need the workspace array `ws_array`. In the full project, the code is first computed using `Symbolics.jl`. It is then processed by a non-Julia symbolic application called [FORM](https://arxiv.org/abs/1409.5223) which reduces the number of floating point operations by about a factor of 8. The output is generated as Fortran code [which is converted back to Julia]. The Fortran compiler prefers arrays to scalars, but I can switch this to generate scalars instead as in your solution.

As an aside, I’ve yet to find a Julia package that can reduce floating point ops via symbolic code optimization as effectively as [FORM](https://www.nikhef.nl/~form/aboutform/aboutform.html). Are you aware of any such Julia package?

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [August 7, 2025, 9:13pm UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421/4 "2025-08-07T21:13:37Z")

</div>

Great!

> [@Audrius-St](#):
>
> As an aside, I’ve yet to find a Julia package that can reduce floating point ops via symbolic code optimization as effectively as [FORM](https://www.nikhef.nl/~form/aboutform/aboutform.html). Are you aware of any such Julia package?

Hmm, I don’t know any package like that, no. Symbolics.jl has the `simplify` function with various options and also allows to define custom term-rewriting rules. I don’t think the default simplifications are quite as powerful as what you are describing though. But there might be packages that implement something like Horner’s method (like here). There is also DynamicExpressions.jl which might be useful in implementing something like this directly in Julia. But if FOAM already exists for that purpose, just calling it like you do now is probably the easiest…

Just out of curiosity: How large of expressions are we talking about and does the rewriting really have a significant impact on the runtime in your final application? It might be best to open a new question with some more details on these details. Other people with more experience on this kind of optimization might also chime in 🙂

---

<div class="post-metadata">

**Author:** ![Audrius-St](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/audrius-st/32/24175_2.png) [@Audrius-St](https://discourse.julialang.org/u/Audrius-St)\
**Post date:** [August 10, 2025, 5:30am UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421/5 "2025-08-10T05:30:48Z")

</div>

> [@Sevi](#):
>
> Hmm, I don’t know any package like that, no. [Symbolics.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Symbolics) has the `simplify` function with various options and also allows to define custom term-rewriting rules. I don’t think the default simplifications are quite as powerful as what you are describing though. But there might be packages that implement something like Horner’s method (like here). There is also [DynamicExpressions.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/DynamicExpressions) which might be useful in implementing something like this directly in Julia. But if FORM already exists for that purpose, just calling it like you do now is probably the easiest…

For now, I’ll stay with using FORM at I’ve invested the time learning to use it and it does the job.

> [@](#):
>
> Just out of curiosity: How large of expressions are we talking about

Well, in one of my functions initially there are 16581 floating point operations. After processing with FORM the number is reduced to 2270. So about a 7.3 reduction in the number of floating point ops. All done in about 1 second on my system. Although it’s a bit like using a sledgehammer to crack a walnut as FORM was initially designed for high energy particle physics where the computations typically involve a far larger number of terms / floating point ops.

> [@](#):
>
> and does the rewriting really have a significant impact on the runtime in your final application?

Yes, it does. Not only runtime, but accuracy at the first step is to convert polynomial and polynomial-like expressions into Horner form followed by extensive common subexpression elimination, greedy optimization, and partial factorization.

> [@](#):
>
> It might be best to open a new question with some more details on these details. Other people with more experience on this kind of optimization might also chime in 🙂

Thanks for your suggestion. Once I’m done I’ll open such a new question. Or I may write an generic Julia wrapper for FORM to generalize my current code.

---

<div class="post-metadata">

**Author:** ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)\
**Post date:** [August 10, 2025, 8:00am UTC](https://discourse.julialang.org/t/nonlinearsolve-how-to-pass-parameters-including-a-mutable-workspace-array/131421/6 "2025-08-10T08:00:06Z")

</div>

Hi, the équivalent for simplifying is ModellingToolKit.jl and it’s structural\_simplify I think.