# Finding roots of an integral equation

Hey out there!

I am trying to solve to following problem.

g(x) = \int_0^x f(y) dy - a

I would like to find the roots. This should be easily solveable by Newtons Method, like the following:

g'(x) = f(x)

    using QuadGK
function mysolve(f, x, a, rtol)
while true
g = quadgk(f, 0, x, rtol=rtol)[1] - a
println("g($x) =$g")
Δx = -g / f(x) # Newton step
x += Δx
abs(Δx) ≤ abs(x)*rtol && break
end
return x
end


This code I took from here and modified it slightly. (thanks to @stevengj)

I was wondering how a problem like this could be solved in the SciML ecosystem (Integrals.jl, NonlinearSolve.jl)?

Something like the following, I suppose?

julia> using QuadGK, NonlinearSolve, Integrals

julia> f(y, p) = y^2 - p.b
f (generic function with 2 methods)

julia> g(x, p) = solve(IntegralProblem(f, (zero(x), x), p), QuadGKJL()).u - p.a
g (generic function with 1 method)

julia> p = (; rtol = 1e-2, a = 0.1, b = 0.5)
(rtol = 0.01, a = 0.1, b = 0.5)

julia> prob = IntervalNonlinearProblem(g, (-1.0, 1.0), p)
IntervalNonlinearProblem with uType nothing. In-place: false
Interval: (-1.0, 1.0)

julia> solve(prob)
retcode: Success
u: -0.2058119300266939

julia> g(ans.u, p)
0.0

julia> using GLMakie

julia> plot(-1:0.1:1, g.(-1:0.1:1, (p,)))


1 Like

You would just have to define the system as f(u, p)=0 and you pass this function to the NonlinearSolve solvers.

Actually I thought it was easier but I have stumbled upon some undesired blocks when trying to code this. This actually works

using NonlinearSolve

f(x) = x^3

function int_equation(u, p)
a, rtol = p
g = quadgk(f, 0, u, rtol=rtol)[1] - a
return g
end

function jacobian(u, p)
return f(u)
end

F = NonlinearFunction(int_equation, jac=jacobian)

p = (1.0, 1e-10)
prob = NonlinearProblem(F, 1.0, p)
sol = solve(prob, Broyden())


but I had to use Broyden since with Newton automatic differentiation was being used and throws an error in the interaction with quadgk. Actually, this is surprising, since I am providing the analytical jacobian so I would expect that this is used instead of doing AD.

Indeed, the following code

using NonlinearSolve

f(x) = x^3

function int_equation(u, p)
a, rtol = p
g = quadgk(f, 0, u, rtol=rtol)[1] - a
return g
end

function jacobian(u, p)
return f(u)
end

F = NonlinearFunction(int_equation, jac=jacobian)

p = (1.0, 1e-10)
prob = NonlinearProblem(F, 1.0, p)
sol = solve(prob, NewtonRaphson(autodiff=false))


throws an error even if I set autodiff to false inside the solver (maybe this keyword argument is being ignored?). The error is the following

ERROR: MethodError: no method matching kronrod(::Type{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1}}, ::Int64)

Closest candidates are:
kronrod(::Any, ::Integer, ::Real, ::Real; rtol, quad)
kronrod(::Type{T}, ::Integer) where T<:AbstractFloat
kronrod(::AbstractMatrix{<:Real}, ::Integer, ::Real, ::Pair{<:Tuple{Real, Real}, <:Tuple{Real, Real}})
...

Stacktrace:
[1] macro expansion
[2] _cachedrule
[3] cachedrule
[4] do_quadgk(f::typeof(f), s::Tuple{…}, n::Int64, atol::Nothing, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
[5] #50
[10] int_equation
@ ~/Desktop/Test/test2.jl:45 [inlined]
[11] NonlinearFunction
@ ~/.julia/packages/SciMLBase/JUp1I/src/scimlfunctions.jl:2297 [inlined]
[12] JacobianWrapper
@ ~/.julia/packages/SciMLBase/JUp1I/src/function_wrappers.jl:106 [inlined]
[13] __value_derivative
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/internal/jacobian.jl:197 [inlined]
[14] JacobianCache
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/internal/jacobian.jl:128 [inlined]
[15] JacobianCache
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/internal/jacobian.jl:116 [inlined]
[16] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generalized_first_order.jl:214
[17] __step!
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generalized_first_order.jl:210 [inlined]
[18] #step!#213
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generic.jl:55 [inlined]
[19] step!
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generic.jl:50 [inlined]
[20] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generic.jl:13
[21] #__solve#212
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generic.jl:4 [inlined]
[22] __solve
@ ~/.julia/packages/NonlinearSolve/FLKyk/src/core/generic.jl:1 [inlined]
[23] #solve_call#44
@ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:612 [inlined]
[24] solve_call
@ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:569 [inlined]
[25] #solve_up#53
@ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1072 [inlined]
[26] solve_up
@ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1066 [inlined]
[27] #solve#52
@ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1060 [inlined]
[28] solve(prob::NonlinearProblem{…}, args::GeneralizedFirstOrderAlgorithm{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1050
[29] top-level scope
@ ~/Desktop/Test/test2.jl:57
Some type information was truncated. Use show(err) to see complete types.


In summary, you can use it with this SciML ecosystem but I would like to ask if anybody knows why it is behaving as I stated in this example because it is a little bit surprising to me. To be honest, I do not know how to turn off AD if that keyword argument does not work and I am not even sure how to enforce the use of the analytical jacobian, since I do not see a clear response in the documentation.

Thanks! Thats awesome.

julia> g(x, p) = solve(IntegralProblem(f, (zero(x), x), p), QuadGKJL()).u - p.a
g (generic function with 1 method)


Now that I see that line it is totally clear, but to be honest I did not come up with the idea to use the IntegralProblem in such a way to set up my function.

using QuadGK
g = quadgk(f, 0, u, rtol=rtol)[1] - a


Thats also something I tried and then feeding this to NonlinearSolve.jl, this also gave rise to errors which I did not completely understood.
And yes giving the analytical Jacobian would be nice!

But the solution from @mkitti is exactly what I was looking for.