How to evaluate an integral not to get infinity?

Question

This is the integral I am trying to solve, and I need to use the answer in a loop which n1a , n2a, n1b and n2b are varying from 1 to 10 for instance. But it gives me infinity although I know that I should have finite answer for my integral. Does any one know how to simplify the answer or solve the integral not to get poles in the denominator?

sol = Integrate[4 A Exp[-(s2 - s1)] Sin[n1a π s1] Sin[n1b π s1] Sin[
n2a π s2] Sin[n2b π s2], {s2, s1, 1}] // Expand

sol2 = Integrate[#, {s1, 0, 1}] & /@ (sol)

Just to give you an example, I have previously solved this integral:

sol1 = Integrate[4 A (s2 - s1)^2 Sin[n1a π s1] Sin[n1b π s1] Sin[
n2a π s2] Sin[n2b π s2], {s2, s1, 1}]

which was giving me an answer having n1a-n1b in the denominator, to overcome that, I have calculated the integral this time for n1a=n1b

sol1 = Integrate[4 A (s2 - s1)^2 Sin[n1b π s1] Sin[n1b π s1] Sin[ n2a π s2] Sin[n2b π s2], {s2, s1, 1}]

and then defined the following function,

fxn[n1a_, n1b_, n2a_, n2b_] := N[ Which[
  n1a == n1b && n2a == n2b,
  1/6 (1 - 3/(n1a^2 \[Pi]^2) - 3/(n2a^2 \[Pi]^2)),
  n1a == n1b && n2a != n2b,
  (4 (1 + (-1)^(n2a + n2b)) n2a n2b)/((n2a^2 - n2b^2)^2 \[Pi]^2),
  n1a != n1b && n2a == n2b,
  (4 (1 + (-1)^(n1a + n1b)) n1a n1b)/((n1a^2 - n1b^2)^2 \[Pi]^2),
  True, -((
  32 (-1 + (-1)^(n1a + n1b)) (-1 + (-1)^(
  n2a + n2b)) n1a n1b n2a n2b)/((n1a^2 - n1b^2)^2 (n2a^2 - 
  n2b^2)^2 \[Pi]^4))]]

to use in the loop, so that I would never get infinity for any value of my variables.

In the new case, the answer is so complicated that applying the same solution is very complicated. I was wondering if anyone have any idea about this.

Answer 1

Work with Limit to get good results.

First do the two integrations with Assumptions (results too long to show here).

sol = Assuming[
   A \[Element] Reals && 0 <= s1 <= 1 && 1 <= n1a <= 10 && 
     1 <= n2a <= 10 && 1 <= n1b <= 10 && 1 <= n2b <= 10, 
     Expand@FullSimplify[
   Integrate[
     4 A Exp[-(s2 - s1)] Sin[n1a \[Pi] s1] Sin[n1b \[Pi] s1] Sin[
     n2a \[Pi] s2] Sin[n2b \[Pi] s2], {s2, s1, 1}]]]

(int[n1a_, n2a_, n1b_, n2b_] = 
     Integrate[#, {s1, 0, 1}, 
       Assumptions -> 
        A \[Element] Reals && 1 <= n1a <= 10 && 1 <= n2a <= 10 && 
        1 <= n1b <= 10 && 1 <= n2b <= 10] & /@ sol) // Timing

This takes about 500 seconds and yields a very long expression.

Whenever there is division by zero in calculating int, add a small parameter d and take the Limit d->0. This workes fine and is very fast. for example:

Limit[int[3 + d, 3, 4, 4], d -> 0] // FullSimplify

(*     (A (E + 75 E \[Pi]^2 - 3 (768 + 335 E) \[Pi]^4 + 
        1225 E \[Pi]^6))/(E (1 + 50 \[Pi]^2 + 49 \[Pi]^4)^2)     *)

Improvement:

Use Limit, only if neccesary.

sol2[n1a_, n2a_, n1b_, n2b_] := 
    Check[int[n1a, n2a, n1b, n2b], 
       Limit[int[n1a + d, n2a, n1b, n2b], d -> 0]] //Quiet // FullSimplify

Answer 2

Not an Answer, but difficult to explain this as a comment.

I did the integral and took its derivative. We get the difference to be 0 as expected.

In[11]:= expr = 
  4 A Exp[-(s2 - s1)] Sin[n1a π s1] Sin[n1b π s1] Sin[
    n2a π s2] Sin[n2b π s2];

In[12]:= sol = Integrate[expr, s2];

In[13]:= D[sol, s2] - expr // FullSimplify

Out[13]= 0

So we can be pretty sure the solution is correct.

Answer 3

Here is a solution inspired by the solution of @Akku14 . This solution is essentially the same as the previous solution, the only real difference being that here we do not evaluate the function once and if we get an infinity we evaluate again and take a limit. We start by defining our integrand, calculating the antiderivative g and defining the function h as follows:

With[
 {f = Exp[-(s2 - s1)] Sin[n1A π s1] Sin[n1B π s1] *
      Sin[n2A π s2] Sin[n2B π s2]},

 With[{g = Integrate[f, {s1, 0, 1}, {s2, s1, 1},
      Assumptions -> {{n1A, n1B, n2A, n2B} ∈ Integers,
        n1A > 0, n1B > 0, n2A > 0, n2B > 0}] // Expand},

  h[n1a_, n1b_, n2a_, n2b_] := 
   Limit[g /. {n1A -> n1a, n1B -> n1b, n2A -> n2a}, n2B -> n2b]
  ]
 ]

Note that we use With to avoid creating global variables; we do a double integration instead of 2 singles; we use {n1A, n1B, n2A, n2B} with capital A,B,C,D instead of {n1a, n1b, n2a, n2b} to define g and we have left off the 4*A multiplier.

The double integration does take a few minutes, but once it is complete our function is defined and when we evaluate h[n1a, n1b, n2a, n2b] it evaluates as a limit, whether direct substitution would give an infinity or not. The advantage is that we never see any infinity.

So, that's it, except for the 4*A multiplier.

Combinations that give infinities

Infinities can be caused by a zero denominator. Mathematica stores "fractions" as powers: 1/n is stored as Power[n,-1]. We can find all of the powers in our antiderivative g that have negative exponents like this:

powers = Cases[g, _Power, ∞];
negPowers = Select[powers, Last[#] < 0 &] // Union;

Now we can find the denominators of the sub-expressions that are negative powers, pick only the denominators that can produce zero for positive integer values of n1a, n1b, n2a, n2b, and normalize the denominators so the coefficient of n1a is positive:

dlist = Denominator /@ negPowers;
dlist = Select[dlist, Reduce[Exists[{n1A, n1B, n2A, n2B},
      And @@ Join[{# == 0, {n1A, n1B, n2A, n2B} ∈ Integers},
        {n1A > 0, n1B > 0, n2A > 0, n2B > 0}]]] &];
dlist = Sort[(Sign[Coefficient[#, n1A]] #) & /@ dlist];
dlist // TableForm

(* {
 {n1A - n1B - n2A - n2B},
 {n1A + n1B - n2A - n2B},
 {n1A - n1B + n2A - n2B},
 {n1A + n1B + n2A - n2B},
 {n1A - n1B - n2A + n2B},
 {n1A + n1B - n2A + n2B},
 {n1A - n1B + n2A + n2B}
}  *)

We find there are 7 denominators that can be zero for some combination of our parameters. Each is a linear combination of our 4 parameters. When we pick a combination of n1A, n1B, n2A, n2B that makes one of the above 7 linear combinations zero, we will get an infinity in our antiderivative g. Any other combination should be okay. More importantly, we can use direct substitution of any 3 of the integers and then evaluate the limit as the other one approaches a critical value. These 7 denominators also tell us that testing with combinations of integers greater than 3 will not produce any new infinities.