SciML / SimpleNonlinearSolve.jl

Fast and simple nonlinear solvers for the SciML common interface. Newton, Broyden, Bisection, Falsi, and more rootfinders on a standard interface.
https://docs.sciml.ai/NonlinearSolve/stable/
MIT License
68 stars 27 forks source link

SimpleNewtonRaphson not using provided Jacobian #63

Closed bclyons12 closed 1 year ago

bclyons12 commented 1 year ago

The documentation for SimpleNewtonRaphson says autodiff is ignored if a user provides a Jacobian. SimpleNewtonRaphson never uses a provided Jacobian, however, while NewtonRaphson does. For example:

using NonlinearSolve, StaticArrays

const R0 = 3.0
const Z0 = 0.0
const a  = 1.0

function f(u, p)
    r = u[1] * a
    st, ct = sincos(u[2])
    dR = R0 + r * ct - p[1]
    dZ = Z0 - r * st - p[2]
    return @SVector[dR, dZ]
end

u0 = @SVector[0.5, 0.0]
p = @SVector[3.3, 0.5]
probN = NonlinearProblem{false}(f, u0, p)

sol = NonlinearSolve.solve(probN, SimpleNewtonRaphson())
println("No Jac:", sol)
println()

function J(u, p)
    println("Using Jacobian")
    r = a * u[1]
    st, ct = sincos(u[2])
    J11 =   a * ct 
    J12 = - r * st
    J21 = - a * st
    J22 = - r * ct
    return @SMatrix[J11 J12; J21 J22]
end
fJ = NonlinearFunction(f; jac=J)
probN = NonlinearProblem{false}(fJ, u0, p)
sol = NonlinearSolve.solve(probN, SimpleNewtonRaphson())
println("Fails to use Jac: ", sol)
println()

sol = NonlinearSolve.solve(probN, NewtonRaphson())
println("Uses Jac: ", sol)

Gives

No Jac:[0.5830951894845299, -1.030376826524313]

Fails to use Jac: [0.5830951894845299, -1.030376826524313]

Using Jacobian
Using Jacobian
Using Jacobian
Using Jacobian
Uses Jac: [0.5830951894592991, -1.030376830345174]
bclyons12 commented 1 year ago

@utkarsh530 Thanks. Any update on when that fix would be released?

utkarsh530 commented 1 year ago

I think it's complete from my side. Waiting for review and merge.