# Calculating and plotting very small values gives errors

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}$$?

• 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. Commented Dec 22, 2023 at 17:59
• 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}]] Commented Dec 22, 2023 at 20:58

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}]]


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.