Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

4 Answers 4

up vote 10 down vote accepted

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]
share|improve this answer
    
Perfect! Thanks for your help. –  jake Dec 13 '13 at 12:37

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.

share|improve this answer
    
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 –  jake Dec 13 '13 at 11:36
    
@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 ? –  tchronis Dec 13 '13 at 11:55
    
Thanks for the edit! –  tchronis Dec 13 '13 at 12:01
    
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? –  jake Dec 13 '13 at 12:04
    
@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. –  tchronis Dec 13 '13 at 12:11

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

share|improve this answer

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 ]
share|improve this answer

Your Answer

 
discard

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.