Finding roots of an integral equation

Specific Domains Numerics
1

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)?

2

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,)))

image
image620×505 15.3 KB

1 Like
3

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)
   @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/weightedgauss.jl:90
  kronrod(::Type{T}, ::Integer) where T<:AbstractFloat
   @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:316
  kronrod(::AbstractMatrix{<:Real}, ::Integer, ::Real, ::Pair{<:Tuple{Real, Real}, <:Tuple{Real, Real}})
   @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:390
  ...

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:564 [inlined]
  [2] _cachedrule
    @ ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:564 [inlined]
  [3] cachedrule
    @ ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:569 [inlined]
  [4] do_quadgk(f::typeof(f), s::Tuple{…}, n::Int64, atol::Nothing, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:7
  [5] #50
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:253 [inlined]
  [6] handle_infinities(workfunc::QuadGK.var"#50#51"{…}, f::typeof(f), s::Tuple{…})
    @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:145
  [7] #quadgk#49
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:252 [inlined]
  [8] quadgk
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:250 [inlined]
  [9] quadgk
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:247 [inlined]
 [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.

4

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.

5 
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.