1
$\begingroup$

I want to plot the difference between $sin(x)$ and its Maclaurin series of order $200$ and more on $[0,1]$.

When I run the code

error[x_] = Sin[x] - Normal[Series[Sin[x], {x, 0, 200}]];
N[error[1], 500]

I get an error message

Internal precision limit $MaxExtraPrecision = 50.` reached while
evaluating
-((1963410391502797893390796845628341196181950954203345943657372226021
7880597655701332158979418609723560351303854634053732931517193406757011
5<<94>>608012596588421747004213863754902836820901759185701726945750964
7034137503676042794159391336019551182168002237136382962149659462965975
315433)/(
2333307299895829892166755070805281249393893425109980068919356607231832
69219941817432042687577510243758632311283704035156880959405397229912<<
94>>483799533750218648640003720134850794852739739361450662906365978872
1780441485881349813724774400000000000000000000000000000000000000000000
000))<<1>>Sin<<1>><<1>>].

also when I try to plot the difference I get another error message

Plot[error[x], {x, 0, 1}, PlotRange -> All]

(2.77086*10^-249
(-1))/4274883284060025564298013753389399649690343788366813724672000000
000000 is too small to represent as a normalized machine number;
precision may be lost.

The question is, how to plot and calculate very small values such as $10^{-1000}$?

$\endgroup$
2
  • $\begingroup$ MMA informs you that you have numbers too small to be represented as normalized real number. A normalized real number is 1.xxxx times 2^exponent, where x and exponent are all binary digits. A nonnormalized number reads: 0.xxxx 2^exponent. This has lower pecision than 1.xxx 2^exponent. $\endgroup$ Commented Dec 22, 2023 at 17:59
  • $\begingroup$ You can eval each individually at high precision, subtract, and live with the cancellation error that loses you substantial precision. In[56]:= sinser[x_] = Normal[Series[Sin[x], {x, 0, 200}]]; In[58]:= N[Sin[1], 500] - N[sinser[1], 500] Out[58]= \ 6.30818921631005374655044903054825814974963992472620376751602664814354\ 07711329570715726165870370939006583396742390256305715*10^-378 For plotting, I'd use high WorkingPrecision and a log scale, like so. ListPlot[Table[Log[10, numdiff[x, 500]], {x, 0, 1, 1/100}]] $\endgroup$ Commented Dec 22, 2023 at 20:58

1 Answer 1

7
$\begingroup$

High-order polynomials with alternating signs are notoriously difficult to evaluate numerically. Such is the case for the sine function:

Sum[(-1)^n x^(2 n + 1)/(2 n + 1)!, {n, 0, ∞}]
(*    Sin[x]    *)

But Mathematica contains a lot of special functions that are implemented with much more stable algorithms. For example, if we take this sum only to order $m$ (in your case, $m=200$):

approx[m_, x_] = Sum[(-1)^n x^(2 n + 1)/(2 n + 1)!, {n, 0, m}]
(*    Sin[x] + (-1)^m x^(2 m + 3) HypergeometricPFQ[{1}, {m + 2, m + 5/2}, -x^2/4]/(2 m + 3)!    *)

the second term in this formula, $\frac{(-1)^m x^{2 m+3}}{(2 m+3)!}{}_1F_2\left(1;m+2,m+\frac52;-\frac{x^2}{4}\right)$, is exactly the finite sum you're looking for:

error[m_, x_] = Sin[x] - approx[m, x]
(*    -(-1)^m x^(2 m + 3) HypergeometricPFQ[{1}, {m + 2, m + 5/2}, -x^2/4]/(2 m + 3)!    *)

Try it out for $m=3$:

error[3, x] // Expand
(*    -x + x^3/6 - x^5/120 + x^7/5040 + Sin[x]    *)

Its evaluation is much more stable, though:

N[error[200, 1], 20]
(*    -2.4038518371210338825*10^-877    *)

ListLinePlot[Table[{x, N[Log[-error[200, x]], 20]}, {x, 0, 1, 1/20}]]

enter image description here

We can even write the ratio of a hypergeometric function and a factorial as a "regularized" hypergeometric:

error[m_, x_] = -(-1)^m Sqrt[π] (x/2)^(2m+3) HypergeometricPFQRegularized[{1}, {m+2, m+5/2}, -x^2/4]

This makes it very simple to look up its properties.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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