Roots of a two variables equation involving a numerical integration

New to Julia
, ,
1

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

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?

2

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

3

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

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, and GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects can do it a little more automatically, as it was built with this usecase in mind.

1 Like