9
$\begingroup$

I'm trying to calculate the function:

CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x]

Unfortunately Mathematica seems to hit a point (x~20) and things become unstable (see plot below), there shouldn't be any infinities in the area, so I am rather confused as to what is going on!

plot of above function

Does anyone know why this is happening and how to fix it?

Improve this question
$\endgroup$
0

4 Answers 4

Reset to default
12
$\begingroup$

You are calculating the small difference between numbers which are getting quite large and the default WorkingPrecision of Plot (which is MachinePrecision, usually about 16) is just not high enough.

So simply increase the WorkingPrecision of Plot, e.g.

Plot[CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x], {x, 0, 40}, 
 WorkingPrecision -> 40, PlotRange -> All]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Perfect! Thanks for your help. $\endgroup$
    – jake
    Commented Dec 13, 2013 at 12:37
8
$\begingroup$

I verified this behavior. I tried to increase $MaxExtraPrecision but I didn't get any improvement.

A solution is to increase WorkingPrecision in Plot:

Plot[CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x], {x, 0, 30}, 
 WorkingPrecision -> 50, PlotRange -> {{0, 30}, {-1, 1}}]

enter image description here

At last I created Cosh and Sinh integrals below with increased WorkingPrecision inside NIntegrate.

chi[z_] := N[EulerGamma + Log[z] + 
             NIntegrate[(Cosh[t] - 1)/t, {t, 0, z}, WorkingPrecision -> 50], 50]

shi[z_] := NIntegrate[Sinh[t]/t, {t, 0, z}, WorkingPrecision -> 50]

Then we can define res[x] as:

res[x_] := N[chi[x] Sinh[x] - shi[x] Cosh[x],50]

And Plot it:

Plot[res[x], {x, 0, 30}, 
 WorkingPrecision -> 50, PlotRange -> {{0, 30}, {-1, 1}}]]

It is a little bit slow but it works.

To check special values you can also use CoshIntegral Evaluation to verify the correctness for chi function above.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for the comment, unfortunately when I tried using the above method I still had the wacky behaviour. I think the Logarithm should be a natural log rather than base 10, and the integral should go to 0 as x->infty $\endgroup$
    – jake
    Commented Dec 13, 2013 at 11:36
  • $\begingroup$ @jake I searched on functions.wolfram.com/GammaBetaErf/CoshIntegral and it seems that log has base 10. Also there you can download a PDF file (after completing a small form) with all information about the function. What version of Mathematica are you using ? $\endgroup$
    – tchronis
    Commented Dec 13, 2013 at 11:55
  • $\begingroup$ I'm using mathematica 9. Thanks for that I will have a flick through. I was under the impression though, that the function Log[] in mathematica was by default a natural logarithm? $\endgroup$
    – jake
    Commented Dec 13, 2013 at 12:04
  • $\begingroup$ @Jake yes you are right. I got confused - it is the natural logarithm. And for this change I get very strange results... I will update shortly. $\endgroup$
    – tchronis
    Commented Dec 13, 2013 at 12:11
  • $\begingroup$ @tchronis You updated your answer at the same time I posted mine! $\endgroup$
    – MikeLimaOscar
    Commented Dec 13, 2013 at 12:34
5
$\begingroup$

An alternative to setting the input precision is to ask for enough output accuracy. The main idea is to set the precision of the input to Infinity and use N to get the desired accuracy -- something like this:

N[f[SetPrecision[x, Infinity], 6]

Here's such an wrapper for a function:

Options[nEval] := {AccuracyGoal -> Automatic, "MaxExtraPrecision" -> Automatic};
nEval[f_, x_?NumericQ, OptionsPattern[]] := 
  Block[{$MaxExtraPrecision = 
     OptionValue["MaxExtraPrecision"] /. Automatic -> $MaxExtraPrecision},
   N[f[SetPrecision[x, Infinity]], OptionValue[AccuracyGoal] /. Automatic -> 6]
   ];

Then 

myF[x_] := CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x];

Plot[nEval[myF, x], {x, 0, 30}, PlotRange -> {{0, 30}, {-1, 1}}]

Mathematica graphics

The precision is adapted to what is needed. 

Plot[nEval[myF, x], {x, 30, 60}, PlotRange -> {Automatic, {-0.35, 0}}]

Mathematica graphics

But one can still get a catastrophic loss of precision, due to the limit $MaxExtraPrecision.

Plot[nEval[myF, x], {x, 30, 100}, PlotRange -> {Automatic, {-0.35, 0}}]

Mathematica graphics

In that case, one can raise the amount of extra precision allowed, even to Infinity. (In some cases, it might take an exorbitant amount of time to finish, but in this case, it is still quick.)

Plot[nEval[myF, x, "MaxExtraPrecision" -> Infinity], {x, 30, 100}, 
 PlotRange -> {Automatic, {-0.35, 0}}]

Mathematica graphics

Improve this answer
$\endgroup$
2
$\begingroup$

A bit of an aside, but the limit behavior of this expression has a readily found simple form:

Simplify[Re@ 
 (*the series expansion has an imaginary part that asymptotically vanishes*)
         Normal@Series[
              CoshIntegral[x] Sinh[x] - SinhIntegral[x] Cosh[x],
                 {x, Infinity, 3}], Assumptions -> {x > 0}]

(* -(2+x^2)/x^3 *)

Now you can plot w/o using excessive precision:

Plot[Piecewise[ {
   {CoshIntegral[x] Sinh[x] - SinhIntegral[x] Cosh[x], x < 10},
   {-((2 + x^2)/x^3), True}
      }], {x, 0, 30},  PlotRange -> All ]
Improve this answer
$\endgroup$

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