Changing variables in integral to avoid infinity

Question

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?

Answer 1

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$

Answer 2

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

Answer 3

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.