# Solving integral equation using DifferentialEquations

**URL:** <https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992>\
**Category:** Numerics\
**Tags:** question\
**Created:** [September 20, 2019, 9:39am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992 "2019-09-20T09:39:01Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 20, 2019, 9:39am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/1 "2019-09-20T09:39:01Z")

</div>

### 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

```julia
A(z) = exp(-z)

```

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

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 20, 2019, 9:59am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/2 "2019-09-20T09:59:42Z")

</div>

> [@Tamas\_Papp](#):
>
> using quadrature

Why not newton solver?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 20, 2019, 10:12am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/3 "2019-09-20T10:12:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 20, 2019, 10:20am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/4 "2019-09-20T10:20:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 20, 2019, 11:26am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/5 "2019-09-20T11:26:24Z")

</div>

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:

```julia
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.)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 20, 2019, 11:34am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/6 "2019-09-20T11:34:22Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 20, 2019, 11:42am UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/7 "2019-09-20T11:42:16Z")

</div>

Yes, for example:

```julia
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
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

```

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [September 20, 2019, 6:26pm UTC](https://discourse.julialang.org/t/solving-integral-equation-using-differentialequations/28992/8 "2019-09-20T18:26:21Z")

</div>

Wouldn’t you have to use the full [Leibniz integration rule](https://en.m.wikipedia.org/wiki/Leibniz_integral_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.