Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have been waiting for Mathematica to give me something for the following integral, errors welcome, but it has been "running" for almost 30 minutes now.

     h[s_] := If[1 < s, 1, 0]
              Abs[1/(b1^2 + b2^2)
               2 E^(-b1 s) ((a1 b1 + a2 b2) E^(
              b1 s) - (a1 b1 + a2 b2) Cos[b2 s] + (-a2 b1 + a1 b2) Sin[
               b2 s]) - h[s]],  {s, 0, Infinity}, Assumptions -> b1 > 0]

Is there a way to simplify this so that Mathematica would not take this long, or is this an entirely lost case? Thank you in advance.

share|improve this question
Actually you can get rid of the function h by breaking the Integral to two, one in the range $(0,1)$ and one in the range $(1,\infty)$. Now if you have more information about the constants in order to get rid of the Abs things might get better. – Spawn1701D Apr 15 '13 at 4:32
Also it might be better instead of h to use UnitStep[s-1]. – Spawn1701D Apr 15 '13 at 4:39
No other assumptions on your other parameters? Note that h[s] is expressible as Boole[1 < s]. – J. M. Apr 15 '13 at 4:39
I am afraid that if you don't get rid of the Abs things will be difficult. – Spawn1701D Apr 15 '13 at 5:10
Actually, you need more than $b_1\gt 0$: you also must require $2(a_1 b_1 + a_2 b_2) = b_1^2 + b_2^2$ to avoid divergence as $s\to \infty$. Generically, this "doesn't happen" because it describes a hypersurface (of zero measure) in $\mathbb{R}^4$. – whuber Apr 15 '13 at 12:46
up vote 3 down vote accepted

Provided $2(a_1 b_1 + a_2 b_2) = b_1^2 + b_2^2$, there exists a closed-form solution that Mathematica can obtain quickly. But you have to do a little analysis along the way :-).

I won't carry out the full solution, but will show how to do the difficult part. The motivation for the first step is recognizing that expressions like $b_1^2+b_2^2$, $a_1 b_1 + a_2 b_2$, etc., have geometric meanings that can be expressed in terms of the lengths of the vectors $(a_1, a_2)$ and $(b_1, b_2)$ and the angles between them. Accordingly, let's introduce polar coordinates for both:

{a1, a2} = a {Cos[\[Alpha]], Sin[\[Alpha]]};
{b1, b2} = b {Cos[\[Beta]], Sin[\[Beta]]};

Remember, for future use, that $a\ge 0$ and $b \ge 0$.

Applying FullSimplify to the integrand yields an expression similar to this one, in which I have manually distributed the exponentials into the numerator:

$$\frac{ \left|2 a e^{-\Re(b s \cos (\beta ))}\cos (\alpha -\beta -b s \sin (\beta ))+ (b h(s)-2 a \cos (\alpha -\beta ))\right|}{|b|}.$$

For the integral to converge, the first term must decay--implying $b\cos(\beta)\gt 0$--and the second term simply cannot be present for arbitrarily large $s$. For such values, $h(s)=1$. We must therefore enforce

$$b - 2 a \cos(\alpha-\beta) = 0.$$

This is equivalent to the condition stated at the beginning of this answer:

(-2 (a1 b1 + a2 b2) + ( b1^2 + b2^2))/b // FullSimplify

$$b-2 a \text{Cos}[\alpha -\beta]$$

When this condition holds, and remembering $a$ and $b$ are nonnegative, the integrand is

$$\frac{2 a e^{-b s \cos (\beta )} |\cos (\alpha -\beta -b s \sin (\beta ))|}{b}$$

To simplify this further, let's re-parameterize it. It looks like it's in the form of a constant times and exponential times a cosine. Here are the new parameters that will make it so:

rules = {c -> 2 a / b, k -> Cot[\[Beta]], u -> b Sin[\[Beta]], \[Gamma] -> (\[Alpha] - \[Beta])};

As a check, the calculation

c Exp[-k u s] Abs[Cos[\[Gamma] - u s]] /. rules  // TraditionalForm

reproduces the previous formula.

Now for the analysis. First, it is clear that a change of variable from s to u s will help: it eliminates one free parameter u (at the cost merely of dividing the answer by $u$ afterwards). Some care will be needed to distinguish the three cases $u \lt 0$, $u=0$, and $u \gt 0$. I will illustrate the case $u \gt 0$. (The case $u \lt 0$ is almost the same and $u = 0$ is easy because the cosine disappears.) Let's illustrate:

With[{k = 1/8, g = 2}, 
 base = Plot[Exp[-k s] Cos[g + s], {s, 0, 20}, Filling -> Axis, 
  FillingStyle -> {RGBColor[.9, .6, .6], RGBColor[.6, .6, .9]}, PlotRange -> Full]];
first = Plot[Exp[-k s] Cos[g + s], {s, 0, 2 \[Pi] - \[Pi]/2 - u}, 
   Filling -> Axis, FillingStyle -> GrayLevel[0.9]]];
Show[base, first]]


Here I have plotted the argument without the absolute value. With the absolute value, we are requesting the total (unsigned) area of the gray, red, and blue parts. The strategy is to compute this area in three pieces according to the colors.

We therefore break the argument of the absolute value into its positive and negative parts. To find out where it is positive--if it's not plainly obvious already--apply Reduce:

Reduce[Cos[\[Gamma] - s] >= 0 , s, Reals]

$C[1]\in \text{Integers}\&\&\frac{1}{2} (-\pi +2 \gamma +4 \pi C[1])\leq s\leq \frac{1}{2} (\pi +2 \gamma +4 \pi C[1])$

This is a countable set of regions indexed by the integer C[1]. Let's perform the integration over a single such region and, afterwards, replace the cosines with $1$ and the sines with $0$:

  Exp[-k s] Cos[\[Gamma] - s], 
    {s, 1/2 (-\[Pi] + 2 \[Gamma] + 4 \[Pi] C[1]), 1/2 (\[Pi] + 2 \[Gamma] + 4 \[Pi] C[1])}, 
  Assumptions ->  C[1] \[Element] Integers && k > 0 && \[Gamma] \[Element] Reals] 
/. {Cos[2 \[Pi] C[1]] -> 1, Sin[2 \[Pi] C[1]] -> 0, C[1] -> n}

$$\frac{e^{-\frac{1}{2} k (\pi +4 n \pi +2 \gamma )} \left(1+e^{k \pi }\right)}{1+k^2}$$

The same treatment of the regions where the integrand is the absolute value of a negative number yields

$$\frac{e^{-k \left(\left(\frac{3}{2}+2 n\right) \pi +\gamma \right)} \left(1+e^{k \pi }\right)}{1+k^2}$$

It remains now only to

  1. Integrate over the first partial period of $\cos(\gamma - s)$ (the gray lobe in the plot) which has a closed form value that Integrate will handily produce,

  2. Add to that the sum of the preceding values starting at an appropriate initial value of $n$ (indexing precisely the red and blue lobes in the plot).

At this juncture I will only demonstrate that the sum has a closed form (it is, in fact, merely a geometric series) and leave it to interested (or more dedicated) readers to determine where to start the summation.

Sum[(E^(-(1/2) k (\[Pi] + 4 n \[Pi] + 2 \[Gamma])) (1 + E^(k \[Pi])))/(1 + k^2)
  + (E^(-k ((3/2 + 2 n) \[Pi] + \[Gamma])) (1 + E^(k \[Pi])))/(1 + k^2), 
  {n, 1, Infinity}] // FullSimplify

$$\frac{e^{-2 k (\pi +\gamma )} \left(e^{k \left(\frac{3 \pi }{2}+\gamma \right)}+e^{\frac{1}{2} k (\pi +2 \gamma )}\right)}{\left(-1+e^{k \pi }\right) \left(1+k^2\right)}$$

To account for the change of variable we must divide this by $u$ and then back-substitute:

%/u /. rules

$$\frac{e^{-2 (\pi +\alpha -\beta ) \text{Cot}[\beta ]} \left(e^{\frac{1}{2} (\pi +2 (\alpha -\beta )) \text{Cot}[\beta ]}+e^{\left(\frac{3 \pi }{2}+\alpha -\beta \right) \text{Cot}[\beta ]}\right) \text{Csc}[\beta ]}{b \left(-1+e^{\pi \text{Cot}[\beta ]}\right) \left(1+\text{Cot}[\beta ]^2\right)}$$

One may go further and re-express this in terms of the $a_i$ and $b_i$, but this is far enough.

share|improve this answer
big +1 for both patience and exposition! – gpap Apr 17 '13 at 16:26
@whuber, I was not expecting such a complete response. Even though I am not sure if I'll be able to understand all of it considering my novice mathematica skills, I want to say : eternal thanks. This is just amazing. – Korben Dallas Apr 20 '13 at 1:09

This integral diverges; you can see that by looking at the behaviour of the integrand for real values of the parameters.

Compile the integrand:

h[s_] := UnitStep[s - 1];
g = Compile[{{a1, _Real}, {a2, _Real}, {b1, _Real}, {b2, _Real}, {s,_Real}},
       Abs[1/(b1^2 + b2^2) 2 E^(-b1 s) ((a1 b1 + a2 b2) E^(b1 s) - (a1 b1 + 
          a2 b2) Cos[b2 s] + (-a2 b1 + a1 b2) Sin[b2 s]) - h[s]], 
  CompilationTarget -> "C"

And then Plot it for (arbitrary) ranges of your a1, a2, b1, b2:

     Plot[g[a1, a2, b1, b2, s], {s, 0, 100}, PlotRange -> All], 
 {a1, -5, 5}, {a2, -5, 5}, {b2, -5, 5}, {b1, 0, 5}]

for large s it asymptotically reaches a non-zero value:

asymptotic divergence

So you need to impose some other constraint if you want this to converge.


So given the additional constraint $ \frac{b_1^2 +b_2^2}{2} = a_1 b_1 + a_2 b_2 $ that @whuber spotted, your integrand becomes:

Abs[1/(b1^2 + 
   b2^2) 2 E^(-b1 s) ((a1 b1 + a2 b2) E^(b1 s) - (a1 b1 + 
     a2 b2) Cos[b2 s] + (-a2 b1 + a1 b2) Sin[b2 s]) - h[s]] /. (b1^2 + b2^2) -> 2 ( a1*b1 + a2*b2)

Abs[(E^(-b1 s) ((a1 b1 + a2 b2) E^(
  b1 s) - (a1 b1 + a2 b2) Cos[b2 s] + (-a2 b1 + a1 b2) Sin[b2 s]))/(a1 b1 + a2 b2) 
  - UnitStep[-1 + s]]

Now, integrating this will be hard: the integral oscillates a lot if $ |b_2|>>1 $ and it gets dampened very slowly if $b_1<<1.$ I don't think it has an analytic solution. Numerically you may get a reasonable result for a range of parameters with this:

intg[a1_?NumberQ, a2_?NumberQ, b1_?NumberQ /; (b1 > 0), b2_?NumberQ] :=
  Abs[(E^(-b1 s) ((a1 b1 + a2 b2) E^(
        b1 s) - (a1 b1 + a2 b2) Cos[b2 s] + (-a2 b1 + a1 b2) Sin[
         b2 s]))/(a1 b1 + a2 b2) - UnitStep[-1 + s]]
  , {s, 0, Infinity}, 
  Method -> {"DoubleExponentialOscillatory", 
    "SymbolicProcessing" -> 0}, WorkingPrecision -> 10, 
  MinRecursion -> 10, MaxRecursion -> 100, AccuracyGoal -> 8]

(which is an adapted example from the NIntegrate Integration Strategies) but, as I said, this breaks down around $ b_1 \sim .01 $ and/or $ b_2 \sim 100 $ or more, i.e.

intg[23, 3, .001, .12]

  Numerical integration converging too slowly; suspect one of the following: singularity, value of the 
     integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

For these, you'd need to do something more exotic like integrate period by period and see if the series of integrals converges or other such things. Sorry I can't help more.

share|improve this answer
Actually, this integral does converge for certain values of the $a_i$ and $b_i$. Your observation is useful in that it indicates Mathematica may be working hard to characterize when the integral is convergent. – whuber Apr 15 '13 at 12:29
@whuber, yes, I completely agree. It's just the OP implied in a comment that $ b_1>0 $ is a sufficient condition for this to converge whereas clearly it isn't. – gpap Apr 15 '13 at 12:44
(+1) I just noticed that comment and added a response describing a (minimal) extra condition that is needed. Imposing it still leaves a difficult-to-evaluate integrand of the form $\exp(-k s)|\cos(u + s)| ds$. – whuber Apr 15 '13 at 12:49
@whuber, you are absolutely right. I added the line to the top 2 (a1*b1 + a2*b2) == b1^2 + b2^2 yet, it has been 15 minutes, and it says "Running" – Korben Dallas Apr 17 '13 at 2:27

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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