I want to write a code in Fortran to solve this integral numerically: $$\int_{1095}^\infty \frac{dx}{x\sqrt{(x+644.153)(4.17 \cdot 10^{-5} x+0.145)}}$$ What is the best method for it?

I tried Monte-Carlo and RK4 but I had to put a big number instead of infinity.

To avoid infinity, I changed the variable $x$ to $1/t$. Then the integral becomes: $$\int_0^{1/1095}\frac{dt}{(t\sqrt{(1/t+644.153)(4.17 \cdot 10^{-5}/t+0.145)}}$$

But this time the problem is that the lower limit of the integral is one of the roots of the denominator. Can you suggest another way of changing variables to avoid this problem, too?


You are almost there, just put $t$ under the square root, and it will become

$ \displaystyle\int_0^{1/1095} {\frac{dt}{\sqrt{(1+644.153t)(4.17 \cdot 10^{-5} + 0.145 t)}}}, $

which eliminates the singularity at $t=0$


There's a closed form solution. For $a,b,c,d,x_0 > 0$, $$ \int_{x_0}^\infty \frac{dx}{x\sqrt{(a+bx)(c+dx)}} = \\ \frac{\log(2 \sqrt{b d (a + b/x_0) (c + d/x_0)} + a d + b c + 2 b d/x_0) - \log(2 \sqrt{b d a c} + a d + b c)}{\sqrt{b d}} $$

a <- 1
b <- 644.153
c <- 4.17e-5
d <- 0.145
x <- 1/1095
(log(2*sqrt(b*d*(a+b*x)*(c+d*x)) + a*d + b*c + 2*b*d*x) - 
    log(2*sqrt(a*b*c*d) + a*d + b*c))/sqrt(b*d) 
# 0.08334367

Using $c = 4.15e^{-5}$ as @nicoguaro, we find the same as his result:

c <- 4.15e-5
(log(2*sqrt(b*d*(a+b*x)*(c+d*x)) + a*d + b*c + 2*b*d*x) - 
    log(2*sqrt(a*b*c*d) + a*d + b*c))/sqrt(b*d) 
# 0.08344436

I think that you don't need a change of variable (yourself) for this problem. Quadpack seems to work just fine for it. It uses a Gauss-Kronrod quadrature.

I tried it (in Python) and it seems to work.

import numpy as np
from scipy.integrate import quad

fun = lambda x: 1/(x*np.sqrt((x + 644.153)*(4.15e-5*x + 0.145)))
inte, err = quad(fun, 1095, np.inf)

This gives 0.08344435887373103 as result with an error of 2.7989753696121526e-10.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.