# How to prove the following integration identity?

I have the following integration that I want to evaluate it using Mathematica.

Evaluate[
Sqrt[θ]/(Sqrt[p] BesselK[1, Sqrt[p θ]])
Integrate[
x BesselK[0, x Sqrt[p θ]]
(Sqrt[θ]/Sqrt[s]
(rTD BesselK[1, rTD Sqrt[s θ]]
(BesselI[0, x Sqrt[s θ]] - BesselK[0, x Sqrt[s θ]]
BesselI[0, Sqrt[s θ]]/BesselK[0, Sqrt[s θ]]) +
BesselI[0, x Sqrt[s θ]]
(x BesselK[1, x Sqrt[s θ]] - rTD BesselK[1, rTD Sqrt[s θ]])
HeavisideTheta[x - rTD] + BesselK[0, x Sqrt[s θ]]
(x BesselI[1, x Sqrt[s θ]] -
rTD BesselI[1, rTD Sqrt[s θ]]) HeavisideTheta[x - rTD])),
{x, 1, Infinity}]]


I am failing to obtain a solution using Mathematica, so I made a manual evaluation. I want to make sure that my solution is correct. It is as follows:

## 2 Answers

I suppose, all parameter p,s,rTD are positive?

Then with not so much effort you can use indefinite integration and Newton-Leibnitz rule.

First you split the integrand into two pieces, first integrated from x=1, second --- from x=rTD:

In[1]={integrand1, integrand2}=Factor@CoefficientList[Sqrt[θ]/(Sqrt[p] BesselK[1,Sqrt[p θ]])(x BesselK[0,x Sqrt[p θ]] (Sqrt[θ]/Sqrt[s] (rTD BesselK[1,rTD Sqrt[s θ]] (BesselI[0,x Sqrt[s θ]]-BesselK[0,x Sqrt[s θ]] BesselI[0,Sqrt[s θ]]/BesselK[0,Sqrt[s θ]])+BesselI[0,x Sqrt[s θ]] (x BesselK[1,x Sqrt[s θ]]-rTD BesselK[1,rTD Sqrt[s θ]]) HeavisideTheta[x-rTD]+BesselK[0,x Sqrt[s θ]] (x BesselI[1,x Sqrt[s θ]]-rTD BesselI[1,rTD Sqrt[s θ]]) HeavisideTheta[x-rTD]))),HeavisideTheta[-rTD+x]]


int1 should be integrated from x=1 to infinity, int2 --- from x=rTD to infinity. We take the indefinite integrals:

In[2]={integral1, integral2} = Integrate[{integrand1, integrand2}, x];


and then use the Newton-Leibnitz rule. Note that both integral1, integral2 decay exponentially at infinity, so we only account for the lower limits:

In[3]=FullSimplify[-Normal[Series[integral1,{x,1,0}]]-Normal[Series[integral2,{x,rTD,0}]],p>0&&s>0&&rTD>0]
Out[3]=(rTD (BesselK[1, rTD Sqrt[p θ]] - (p BesselK[0, Sqrt[p θ]] BesselK[1, rTD Sqrt[s θ]])/(Sqrt[p s] BesselK[0, Sqrt[s θ]])))/(p (-p + s) BesselK[1,Sqrt[p θ]])


This one coincides with your result.

NB Of course, not all definite integrals can be done via Newton-Leibnitz rule, but it is always a good idea to check if the indefinite integral can be calculated. The reason is that, in my opinion, nowadays Mathematica is much stronger (and faster) in taking the indefinite integrals then in taking the definite ones.

• "Mathematica is much stronger (and faster) in taking the indefinite integrals" - I would say instead that Mathematica is more cautious in evaluating definite integrals, and will refuse to return a result if it cannot ascertain things like continuity of the antiderivative. As you have said, applying Newton-Leibniz can sometimes lead to trouble. Jan 11, 2021 at 7:08
• Actually, my statement was different. There are many examples where the indefinite integral can not be expressed in terms of any known functions, while the definite one can still be taken. Mostly via complex analysis methods. On the other hand, if the indefinite integral can be taken, with due attention to crossing the singularities, one can "always" (with some reservations to the complexity of the appearing functions) apply the Newton-Leibnitz rule. Jan 11, 2021 at 7:41

This might get you most of the way there. First break up the integrand into 6 pieces:

integrand = Sqrt[θ]/(Sqrt[p] BesselK[1, Sqrt[p θ]]) x BesselK[0, x Sqrt[p θ]] (Sqrt[θ]/
Sqrt[s] (rTD BesselK[1, rTD Sqrt[s θ]] (BesselI[0, x Sqrt[s θ]] -
BesselK[0, x Sqrt[s θ]] BesselI[0, Sqrt[s θ]]/BesselK[0, Sqrt[s θ]]) +
BesselI[0, x Sqrt[s θ]] (x BesselK[1, x Sqrt[s θ]] - rTD BesselK[1, rTD Sqrt[s θ]]) HeavisideTheta[x - rTD] +
BesselK[0, x Sqrt[s θ]] (x BesselI[1, x Sqrt[s θ]] - rTD BesselI[1, rTD Sqrt[s θ]]) HeavisideTheta[x - rTD])) // Expand


Integrate each piece (with the integrands with HeavysideTheta functions just have the limits of integration changed):

i1 = Integrate[(rTD x θ BesselI[0, x Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[1, rTD Sqrt[s θ]])/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]),
{x, 1, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]

i2 = -Integrate[(rTD x θ BesselI[0, Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[0, x Sqrt[s θ]] BesselK[1, rTD Sqrt[s θ]])/(Sqrt[p] Sqrt[s] BesselK[0, Sqrt[s θ]] BesselK[1, Sqrt[p θ]]),
{x, 1, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]

i3 = -Integrate[(rTD x θ BesselI[1, rTD Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[0, x Sqrt[s θ]] )/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]),
{x, rTD, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]

i4 = Integrate[((x^2) θ BesselI[1, x Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[0, x Sqrt[s θ]] )/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]),
{x, rTD, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]

i5 = -Integrate[(rTD x θ BesselI[0, x Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[1, rTD Sqrt[s θ]])/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]),
{x, rTD, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]

i6 = Integrate[((x^2) θ BesselI[0, x Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[1, x Sqrt[s θ]] )/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]),
{x, rTD, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]


The integrations work for i1, i2, i3, and i5 but not for i4 and i6. We add the integrands for i4 and i6 and try again:

i46 = Integrate[(x^2 θ BesselI[1, x Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[0, x Sqrt[s θ]])/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]) +
(x^2 θ BesselI[0, x Sqrt[s θ]] BesselK[0, x Sqrt[p θ]] BesselK[1, x Sqrt[s θ]])/(Sqrt[p] Sqrt[s] BesselK[1, Sqrt[p θ]]) // FullSimplify,
{x, rTD, ∞}, Assumptions -> {p > s > 0, θ > 0, rTD > 0}]

(* (rTD BesselK[1, rTD Sqrt[p θ]])/(p s BesselK[1, Sqrt[p θ]]) *)


So FullSimplify helps with i46 because of the identity that Mathematica knows:

BesselI[ν, z] BesselK[ν + 1, z] + BesselI[ν + 1, z] BesselK[ν, z] == 1/z


Now add up all of the terms (and convert the hypergeometric functions in i1):

integral = i1 + i2 + i3 + i5 + i46 /.
Hypergeometric0F1Regularized[1, (s θ)/4] -> BesselI[0, Sqrt[s θ]] /.
Hypergeometric0F1Regularized[2, (s θ)/4] -> (2 BesselI[1, Sqrt[s θ]])/Sqrt[s θ]
integral = FullSimplify[integral, {p > s > 0, \[Theta] > 0}]

(* (rTD (-s BesselK[0, Sqrt[s θ]] BesselK[1, rTD Sqrt[p θ]] + Sqrt[p s] BesselK[0, Sqrt[p θ]] BesselK[1, rTD Sqrt[s θ]]))/(p (p - s) s BesselK[0, Sqrt[s θ]] BesselK[1, Sqrt[p θ]]) *)


$$\frac{r_{TD} \left(\sqrt{p s} K_0\left(\sqrt{p \theta }\right) K_1\left(r_{TD} \sqrt{s \theta }\right)-s K_0\left(\sqrt{s \theta }\right) K_1\left(r_{TD} \sqrt{p \theta }\right)\right)}{p s (p-s) K_1\left(\sqrt{p \theta }\right) K_0\left(\sqrt{s \theta }\right)}$$

Note that

Hypergeometric0F1Regularized[1, (s θ)/4] // FunctionExpand
(* BesselI[0, Sqrt[s] Sqrt[θ]] *)
Hypergeometric0F1Regularized[2, (s θ)/4] // FunctionExpand
(* (2 BesselI[1, Sqrt[s] Sqrt[θ]])/(Sqrt[s] Sqrt[θ]) *)


The answer matches what you have:

-((rTD BesselK[1, rTD Sqrt[p θ]])/(p (p - s) BesselK[1, Sqrt[p θ]])) +
(rTD Sqrt[p s] BesselK[0, Sqrt[p θ]] BesselK[1, rTD Sqrt[s θ]])/
(p (p - s) s BesselK[0, Sqrt[s θ]] BesselK[1, Sqrt[p θ]]) ==
rTD BesselK[0, Sqrt[p θ]]/((BesselK[1, Sqrt[p θ]]) (s - p)) (BesselK[1, rTD Sqrt[p θ]]/(p BesselK[0, Sqrt[p θ]]) -
BesselK[1, rTD Sqrt[s θ]]/(Sqrt[p s] BesselK[0, Sqrt[s θ]])) // FullSimplify

(* True *)


There is certainly a simpler way to do this but this is just an example as to what one might need to go through to solve a complicated integral.