Hi! I’m new to using NonlinearSolve.jl and trying to understand how to provide an analytical jacobian to the problem. I realize that this is unnecessary if we have autodiff in the mix since that will be enough - but in cases where I have say registered functions (and corresponding registered jacobians) I want to understand how to explicitly supply the jacobian to the problem.
Here’s what I have as an example:
# Starting with a simple example
function f(u, p)
@. u*u - p
end
u0 = [1.0, 1.0, 1.0]
p = [2.0, 1.0, 4.0]
prob = NonlinearProblem(f, u0, p)
julia> sol = solve(prob, NewtonRaphson())
retcode: Success
u: 3-element Vector{Float64}:
1.4142135623730951
1.0
2.000000000000002
Now if I try to give it the jacobian and setup a NonlinearFunction (was not very sure how else to setup a NonlinearProblem with a jacobian than doing this):
df(u,p) = 2.0 .* u
h = NonlinearFunction(f, jac = df)
prob = NonlinearProblem(h, u0, p)
If I try solving this problem now I get this:
julia> sol = solve(prob, NewtonRaphson())
retcode: Unstable
u: 3-element Vector{Float64}:
1.0
1.0
1.0
And looking at this I tried specifying to the NewtonRaphson algorithm concrete_jac = true but that still does not seem to help.
julia> sol = solve(prob, NewtonRaphson( ;concrete_jac = true))
retcode: Unstable
u: 3-element Vector{Float64}:
1.0
1.0
1.0
All this is using the latest NonlinearSolve version:
[8913a72c] NonlinearSolve v3.9.1
Any thoughts or suggestions would be very appreciated!