2
$\begingroup$

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.

Improve this question
$\endgroup$
1

3 Answers 3

Reset to default
1
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you so much. $\endgroup$ Commented Sep 29, 2017 at 20:15
  • $\begingroup$ I replaced n1a with n1a+0.0001 in (int), and I think it is doing the same thing as the limit will do. (so that we do not create 0/0, and it can simplify the nominator with the denominator ) $\endgroup$ Commented Sep 30, 2017 at 18:46
1
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Are you sure that the denominator poles should get cancelled ? $\endgroup$ Commented Sep 29, 2017 at 15:23
  • $\begingroup$ Thank you Lotus, I have added an example to clarify my question,Could you please take a look at it again? $\endgroup$ Commented Sep 29, 2017 at 15:48
1
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ LouisB, my function "int" is only evaluated once, since I used "Set" instead of "SetDelayed".If it yields infinity, it it not calculated again (which would take 500 seconds), but the once calculated result is called with the slightly changed parameter na1 + d. Your Limit is called every time, which can take much time for some parameter combinations. My parameters are not restricted to Integers. $\endgroup$ Commented Sep 30, 2017 at 13:08

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.