# Roots of a two variables equation involving a numerical integration

**URL:** <https://discourse.julialang.org/t/roots-of-a-two-variables-equation-involving-a-numerical-integration/37544>\
**Category:** New to Julia\
**Tags:** question, plotting, roots\
**Created:** [April 14, 2020, 4:08am UTC](https://discourse.julialang.org/t/roots-of-a-two-variables-equation-involving-a-numerical-integration/37544 "2020-04-14T04:08:39Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Daniel\_Tamayo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daniel_tamayo/32/10191_2.png) [@Daniel\_Tamayo](https://discourse.julialang.org/u/Daniel_Tamayo)\
**Post date:** [April 14, 2020, 4:08am UTC](https://discourse.julialang.org/t/roots-of-a-two-variables-equation-involving-a-numerical-integration/37544/1 "2020-04-14T04:08:39Z")

</div>

I’m in troubles to find the zeros of the next equation of two variables:  
(w - k)A + B = 0  
The issue is that the terms ‘_A_’ and ‘_B_’ are not numbers but are integral of certain functions ‘_f(k,w,p)_’ and ‘_g(k,w,p)_’ involving the three variables in a very complicate way (since the integration is only on p, it eventually must disappear),

A(k,w) = \int\_0^{\infty} f(k,w,p) dp  
B(k,w) = \int\_0^{\infty} g(k,w,p) dp

I have tried to use a double subsequent _for loop_ to the variables k and w, so that I can integrate numerically over p only. Then I check which pairs (k, w) satisfy the equation by setting a range close to zero. The next is a simple case that works, to illustrate my method

```julia
using QuadGK

eps = 0.01
T = 100
c = 0.0008/(8*pi^2)

nF(p) = 1/(exp(p/T) + 1) 
nB(p) = 1/(exp(p/T) - 1)

for w in 0.01:0.01:5
    for k in 0.01:0.01:w
        if w > k
            L1(p) = log( (p + 0.5*(w+k))/(p + 0.5*(w-k)) ) - log( (p - 0.5*(w+k))/(p - 0.5*(w-k)) )  
            L2(p) = log( (p + 0.5*(w+k))/(p + 0.5*(w-k)) ) + log( (p - 0.5*(w+k))/(p - 0.5*(w-k)) )
            
            
            A(p) = (c/k^2)*( (4*p + (w^2 - k^2)*L1(p)/(2*k) - 2*p*w*log((w+k)/(w-k))/k + w*p*L2(p)/k)*nF(p) + (4*p + (w^2 - k^2)*L1(p)/(2*k) - 2*p*w*log((w+k)/(w-k))/k + p*w*L2(p)/k )*nB(p) )
            B(p) = (c/k^2)*( (2*p*(w^2 - k^2)*log((w+k)/(w-k))/k - p*(w^2-k^2)*L2(p)/k - 4*w*p - w*(w^2-k^2)*L1(p)/(2*k))*nF(p) + (2*p*(w^2-k^2)*log((w+k)/(w-k))/k - p*(w^2-k^2)*L2(p)/k - w*(w^2-k^2)*L1(p)/(2*k) - 4*w*p)*nB(p) )

            IA = quadgk(A, 5, Inf, rtol=1e-3)
            IB = quadgk(B, 5, Inf, rtol=1e-3)
            
            if -eps < ((w-k)*(1+IA[1]) + IB[1]) < eps
            println(k, w)
            
            end
        end
    end
end

```

However, it has not been efficient. Do you know of any method to perform this calculation in a better way?  
Another problem I have is, for slightly more complex A and B functions, a certain configurations of values of k,w and p, I get indeterminacies by evaluating logarithms of negative numbers (I’m only interested in real solutions).  
How can I evaluate the sign of the expression inside the logarithms to discard the pairs (k, w) that lead to these unwanted ones?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 14, 2020, 4:12am UTC](https://discourse.julialang.org/t/roots-of-a-two-variables-equation-involving-a-numerical-integration/37544/2 "2020-04-14T04:12:49Z")

</div>

One option is to use a ML package to try to find a minimum of `((y - x)A + B )^2`. Since global minima of the function happen when the original equals 0, that should work (ish).

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [April 14, 2020, 5:25am UTC](https://discourse.julialang.org/t/roots-of-a-two-variables-equation-involving-a-numerical-integration/37544/3 "2020-04-14T05:25:02Z")

</div>

> [@Daniel\_Tamayo](#):
>
> However, it has not been efficient.

First, put the whole thing in a function. The globals may be extremely slow (or may not change anything! You never know…)

> [@Daniel\_Tamayo](#):
>
> quadgk(A, 5, Inf, rtol=1e-3)

Since `quadgk` is adaptive quadrature, I have always found it to be extremely slow for this sort of thing. It is usually better to use gaussian quadrature or something along those lines if possible when solving within an optimizer. You can get the nodes manually with someting like [GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl), and [GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects](https://github.com/QuantEcon/Expectations.jl) can do it a little more automatically, as it was built with this usecase in mind.