# Solve for upper/lower bound of a numerical definite integral
**URL:** https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164
**Category:** General Usage
**Tags:** question, math
**Created:** [March 29, 2021, 2:14pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164 "2021-03-29T14:14:21Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@carloslesmes](https://discourse.julialang.org/u/carloslesmes)
#### Post date: [March 29, 2021, 2:14pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/1 "2021-03-29T14:14:21Z")
Hi, I’m trying to solve for b, given f, a and k, please help, Is there a way to do it in julia?
```julia
quadgk(f,a,b)=k
```
Thanks
---
### Author:  [@klaff](https://discourse.julialang.org/u/klaff)
#### Post date: [March 29, 2021, 2:35pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/2 "2021-03-29T14:35:55Z")
Possibly [https://github.com/JuliaMath/Roots.jl](https://github.com/JuliaMath/Roots.jl)
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 29, 2021, 2:59pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/3 "2021-03-29T14:59:57Z")
There’s always a Newton iteration:
```julia
b = some_initial_guess
tol = sqrt(eps())
while true
δb = (quadgk(f,a,b, rtol=tol*0.1)[1] - k) / f(b)
b -= δb
abs(δb) ≤ tol * abs(b) && break
end
```
Note that I’m using the fact that \frac{d}{db} \int\_a^b f(x) dx = f(b).
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 29, 2021, 5:00pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/4 "2021-03-29T17:00:53Z")
@stevengj, this is excellent.
Below a summary of different solutions, including yours and your advice to run Optim properly (_even if not recommended for this problem kept it for sake of completeness_).
```julia
f(x) = exp(x)*sin(x)
# forward modelling
a = 0
# b = π : is UNKNOWN to be found given f, a and k
k = 0.5*(1.0 + exp(1)^π) # k = 12.0703463163896
# SOLUTION-1: Newton solution by @stevengj
using QuadGK
b = 1.0 # initial guess
tol = sqrt(eps())
while true
δb = (quadgk(f,a,b, rtol=tol*0.1)[1] - k) / f(b)
b -= δb
abs(δb) ≤ tol * abs(b) && break
end
println("Solution = ", b) # 3.14159268 / solution π= 3.14159265
# SOLUTION-2: use Optim to minimizing g(x)^2
using Optim
g(x) = (quadgk.(f,a,x...))[1] - k
h(x) = ((g.(x)).^2)[1]
dh(x) = (2*g.(x).*f.(x))[1] # derivative
# Checks:
g(pi) ≈ 0 # true
h(pi) ≈ 0 # true
dh(pi) ≈ 0 # true
b = [1.0] # initial value / solution π= 3.14159265
result = optimize(h, b, Newton(), Optim.Options(g_tol = 1e-16))
Optim.minimizer(result) # 3.14159265
result = optimize(h, dh, b, Newton(); inplace=false) # using derivative
Optim.minimizer(result) # 3.1415795 # does not improve accuracy
# SOLUTION-3: use NLsolve to find roots
using NLsolve
g(x) = (quadgk.(f,a,x)...)[1] - k
dg(x) = f.(x) # derivative
# Check:
g([pi]) ≈ 0 # true
b = [1.0] # initial value / solution π= 3.14159265
nlsolve(g, b; ftol=1e-16) # 3.14159264
nlsolve(g, dg, b; ftol=1e-16) # 3.14159263 # same accuracy...
# SOLUTION-4: use DifferentialEquations and interpolation
using DifferentialEquations, Plots, Dierckx
u0 = 0; tspan = (a, 4.) # initial conditions & timespan
du(u,p,t) = f(t) # define the system
prob = ODEProblem(du, u0, tspan) # define the problem
sol0 = solve(prob, dtmax=0.01) # solve it
i = argmin(abs.(sol0.u .- k)) # index of time-step the closest to solution
tsol = LinRange(sol0.t[i-1], sol0.t[i+1],20) # spline it around the solution
spl = Spline1D(tsol, sol0.(tsol) .- k)
b = mean(roots(spl)) # 3.1415926536
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 29, 2021, 6:25pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/5 "2021-03-29T18:25:33Z")
> [@rafael.guerra](#):
>
> Do you know if it would be possible to use optim in such case?
Yes; it looks like the problem with your code is that `Optim.optimize` is expecting a vector of initial guesses, not a scalar. Try passing `[b]` instead of `b`.
**However** , this is a root-finding problem `g(b) = 0`. It is _almost always_ a mistake to try to solve a root-finding problem by using an optimization algorithm to minimize `|g|²`. You should use a root-finding algorithm like those in [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl).
Moreover, in this case you know the derivative analytically, as I mentioned above, in which case you should certainly provide it to the root-finding algorithm.
However, for a 1d root-finding problem (`b` is a scalar) with an analytically known derivative, all of the generic root-finding packages will basically boil down to a Newton iteration like the one I wrote. The clever algorithms are mainly for the case where you don’t know the derivative, or have lots of derivatives (a big Jacobian) and can’t affort to compute them all or to invert the Jacobian.
Alternatively, in 1d, especially for smooth functions, there are often more clever algorithms available. For example, you could use ApproxFun.jl to construct a high-accuracy polynomial approximation of your integrand f, call cumsum to compute its integral, and use the `ApproxFun.roots` function to find all of the places where the integral equals `k`.
---
### Author:  [@apo383](https://discourse.julialang.org/u/apo383)
#### Post date: [March 29, 2021, 7:17pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/6 "2021-03-29T19:17:12Z")
> [@stevengj](#):
>
> It is _almost always_ a mistake to try to solve a root-finding problem by using an optimization algorithm to minimize `|g|²` .
Rather than minimize `|g|²`, one could apply `g(b)=0` as an equality constraint to an optimization algorithm that supports such constraints, e.g. `NLopt.jl` or `JuMP.jl`. I often find it convenient to deal with roots the same way I deal with optimization. Sometimes you can leave out an objective function altogether, or just specify a constant like 0, to be “minimized” subject to the constraints (to some tolerance). However, I’ve never compared performance to a root finder.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 29, 2021, 9:10pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/7 "2021-03-29T21:10:08Z")
@stevengj, thank you for your top quality advice.
Tried to follow it the best possible and have updated code above for 3 different methods: Newton by yourself, Optim and NLsolve.
Including the derivative did not seem to help much, but any blame should surely be on me. BR
\*_PS: code above now runs but better not shaking it too much…_
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 29, 2021, 9:47pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/8 "2021-03-29T21:47:48Z")
> [@apo383](#):
>
> Rather than minimize `|g|²` , one could apply `g(b)=0` as an equality constraint to an optimization algorithm that supports such constraints, e.g. `NLopt.jl` or `JuMP.jl` . I often find it convenient to deal with roots the same way I deal with optimization. Sometimes you can leave out an objective function altogether, or just specify a constant like 0, to be “minimized” subject to the constraints (to some tolerance). However, I’ve never compared performance to a root finder.
Yes, solving the optimization problem \min\_x 1 subject to the constraint g(x)=0 is, of course, exactly equivalent to the root-finding problem g(x)=0, and a sufficiently clever NLP algorithm will perform steps essentially equivalent to Newton steps (perhaps limited by a trust region). They almost certainly aren’t going to be better than Newton steps, however, and in many cases will be worse (because they may spend a lot of effort doing computational steps that are only needed for more general NLPs).
In a single-variable problem where the derivative is known, if you aren’t going to do something more clever (e.g. Chebyshev fits, complex analysis, continued fractions, etc.) you might as well just use Newton. (Either implementing it yourself — it’s just a few lines — or calling a canned implementation like the one in [Roots.jl](https://github.com/JuliaMath/Roots.jl)).
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 30, 2021, 7:18am UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/9 "2021-03-30T07:18:16Z")
There is another way. Define R as
$$\int\_a^{R(x)}f=x$$
you get
$$\frac{d}{dx}R(x) = \frac{1}{f(R(x))},\quad R(a)=0 $$
You can use an ODE solver to get b=R(k). This is the main trick of my [paper](https://arxiv.org/abs/1504.06873) on PDMP.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 30, 2021, 9:09am UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/10 "2021-03-30T09:09:49Z")
@rveltz, it is only me, or your Latex is not displaying properly?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 30, 2021, 9:18am UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/11 "2021-03-30T09:18:40Z")
It gives:
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 30, 2021, 10:49am UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/12 "2021-03-30T10:49:11Z")
@rveltz, mathematically this looks pretty much like @stevengj’s solution above. The difference is that the differential equations solver will not use Newton’s method?
Would like to try this when time allows, but would be surprised if it beats Newton in terms of accuracy and fast convergence.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 30, 2021, 11:02am UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/13 "2021-03-30T11:02:18Z")
I would say it is better in that you dont need to recompute \int\_a^bf if you know that the solution b\_{sol} is greater (smaller) than b. You “waste” less computation
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [March 30, 2021, 11:46am UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/14 "2021-03-30T11:46:28Z")
Yeah, it seems to me that there is a lot of redundancy in re-computing the full integral at each step. It should be enough to compute the integral only over the ‘correction interval’. (One should probably take some care with accumulated errors, however.)
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [April 1, 2021, 6:05pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/15 "2021-04-01T18:05:40Z")
@rveltz, in your nice solution how do you handle the cases with singularities where the function `f()` has zeros and `1/f()` explodes to `∞`?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 1, 2021, 6:53pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/16 "2021-04-01T18:53:07Z")
Hum, I dont know! I guess the Newton solver will sufer as well. I’d say you have to use a stiff ODE solver that can handle finite time explosion.
---
### Author:  [@apo383](https://discourse.julialang.org/u/apo383)
#### Post date: [April 1, 2021, 6:59pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/17 "2021-04-01T18:59:31Z")
The topic “Solve equation” doesn’t do justice to the interesting things discussed here. @carloslesmes can I suggest editing the topic title to something more specific? For example, “Solve for upper/lower bound of a numerical definite integral” or some such.
---
### Author:  [@jacobadenbaum](https://discourse.julialang.org/u/jacobadenbaum)
#### Post date: [April 1, 2021, 7:37pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/18 "2021-04-01T19:37:29Z")
Am I looking at this crazy, or wouldn’t the correct boundary condition be R(0) = a instead of R(a) = 0?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 1, 2021, 8:36pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/20 "2021-04-01T20:36:46Z")
nope you are right! Thank you
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [April 1, 2021, 9:27pm UTC](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164/21 "2021-04-01T21:27:37Z")
> [@rveltz](#):
>
> I guess the Newton solver will sufer as well
Newton did not suffer for `f(x)=exp(x)sin(x)` in the example above with a zero at the solution `pi`. However, someone else is suffering to have your method working for that `f(x)`… 😓
[Next page](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164.md?page=2)