Solving integral equation using DifferentialEquations

Specific Domains Numerics
1

Problem

I need to solve an equation in x which can be simplified as an MWE into

x = \int_{x}^\infty A(z) dz

Importantly, A(z) is decreasing, integrable, and goes to 0 at infinity.

Possible solution concept

I thought I would define

K(x) = \int_x^\infty A(z) dz - x

and evaluate eg K(0) using quadrature, then use an ODE solver on

K'(x) = -A(x) - 1

starting from eg 0 and see where it crosses 0.

I just don’t know how to implement this with the DifferentialEquations ecosystem. I looked at the integrator interface, but all I could think of is finding where I cross 0, then backtracking.

Other suggestions for solution methods are appreciated.

For an MWE, consider

A(z) = exp(-z)

for which the (numerically obtained) solution should satisfy z = e^{-z}.

2

Why not newton solver?

2 Likes
3

Sorry, I don’t know what that is. An example would be appreciated.

4

if you know how to compute the integral with quadrature, why do you want to use DiffEq? and chose a (non)linear solver?

5

Use Newton and a quadrature routine. That is, you are trying to find a root of

f(x) = x - \int_x^\infty A(z) dz

The derivative of this is f'(x) = 1 + A(x). So, you can just repeatedly take Newton steps

x \leftarrow x - \frac{f(x)}{f'(x)}

where f(x) is evaluated with a quadrature routine like QuadGK:

using QuadGK
f(x) = x - quadgk(A, x, Inf, rtol=1e-4)[1]

for example. (QuadGK handles the semi-infinite integral for you by a coordinate transformation.)

9 Likes
6

OK, I see that you and @stevengj suggest that I just use a univariate solver directly on K(x), and Newton should be a good choice since K'(x) is available at no cost. Got it, thanks!

7

Yes, for example:

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

converges to about 10 digits in only 4 steps for A(x) = e^{-x} / (x^2 + 1) starting at an initial guess of x=1:

julia> x = mysolve(t -> exp(-t) / (t^2 + 1), 1, 1e-4)
f(1) = 0.9033475191037287
f(0.23699872265724586) = -0.17704370803589597
f(0.3383383890915101) = -0.005483914249519217
f(0.3416828040022412) = -5.745324205275182e-6
0.34168631519454107

julia> x - quadgk(t -> exp(-t) / (t^2 + 1), x, Inf, rtol=1e-7)[1]
-4.3091752388590976e-11
5 Likes
8

Wouldn’t you have to use the full Leibniz integration rule here? I believe the correct derivative of K(x) should then be -2A(x) - 1.
Edit: Ah, never mind. The A inside the integral only depends on z, not x explicitly, so you were right.