# Funny behaviour when plotting a polynomial of high degree and large coefficients

I am trying to plot a polynomial of degree $29$ on the domain $[0,1]$, with fairly large coefficients:

poly[z_] = -1.1126829840302355 + 113.28783661058498 z -
9878.742379213338 z^2 + 584715.8646149524 z^3 -
2.280647160113914*^7 z^4 + 6.176933520283158*^8 z^5 -
1.217581378062843*^10 z^6 + 1.811582263559531*^11 z^7 -
2.0920133095023196*^12 z^8 + 1.9158620445090305*^13 z^9 -
1.414955836879161*^14 z^10 + 8.539129365981781*^14 z^11 -
4.254563022912091*^15 z^12 + 1.764182366816184*^16 z^13 -
6.125080435776876*^16 z^14 + 1.7883504482275766*^17 z^15 -
4.403320010637656*^17 z^16 + 9.154951756734264*^17 z^17 -
1.6068672087698447*^18 z^18 + 2.3765393965161196*^18 z^19 -
2.950846281328122*^18 z^20 + 3.0579497598096415*^18 z^21 -
2.6220913470110597*^18 z^22 + 1.837651151556163*^18 z^23 -
1.0344252684666292*^18 z^24 + 4.5602474296077024*^17 z^25 -
1.5155510563521117*^17 z^26 + 3.568596763872067*^16 z^27 -
5.304183668348243*^15 z^28 + 3.7404713997980006*^14 z^29


The problem is that when I try to plot the polynomial, all seems fine near 0, but then midway through the unit interval I observe some erratic behaviour as shown below. Can anyone tell me why this is and how I can avoid it?

Plot[{poly[z]}, {z, 0, 1}, WorkingPrecision -> precision]


• Have you considered that you're working with a 29th degree polynomial with insanely huge coefficients? It's not surprising to me that you get a plot like that. It may seem "nice" early on, but that's because z^29 (and other high-order terms) is very small near 0. Near approximately 0.85, it rapidly approaches 1. Commented Mar 18, 2012 at 2:14
• @Mithc, I will have to disagree with you. Just because the terms of the ploynomial get large as z goes to 1 does not mean you should expect this type of erratic behaviour. Indeed the polynomial can have no more than 29 real roots! So why you expect a plot like this is beyond me. You should not have voted the question down. I suspect a numerical issue here and was hoping somebody could help me pin it down.
– aukie
Commented Mar 18, 2012 at 3:09
• It's cancellation error. Per advice in the response, use higher precision and/or rationalize the coefficients. And be sure you know what you are doing in terms of what might comprise an "expected" outcome, because the result might not mean much if the actual coefficients really have error intervals associated with them. Commented Mar 19, 2012 at 1:16

If you Rationalize your real numbers you will be able to use Mathematica's arbitrary precision engine:

poly2 = Rationalize[poly[z], 0];

Plot[poly2, {z, 0, 1}, WorkingPrecision -> 50]


## Arbitrary and machine precision

Mathematica has two kinds of numeric calculations: machine precision, and arbitrary precision. Machine precision is fast but is limited to 53 binary (≈16 decimal) digits(1)(2) and may lose precision during a calculation. Mathematica also does not track the precision of a result.

Numbers entered as 1.234 or 1234 are taken to be machine precision.

Arbitrary precision is slower, but Mathematica will track the precision of calculations, often using more calculations as necessary to preserve precision, and print the precision of the result.

Exact values such as 1234 and 1/2 can be used in arbitrary precision calculations. Numbers can also be entered with e.g. 1.23420 specifying 20 digits of precision, and these will automatically use arbitrary precision if all other values are either exact or arbitrary.

Precision can be checked with the function Precision:

Precision /@ {1.234, 1234, 1.23420, 7}

{MachinePrecision, MachinePrecision, 20., ∞ }


Precision can only be preserved if all values in a calculation have at least that precision. Also, arbitrary precision arithmetic may be used with numbers having a precision less than MachinePrecision -- Mathematica will show the true precision of the result.

Precision[1.234 + 7]

Precision[1.23420 + 1.23412]

MachinePrecision

12.301


Precision can be set with SetPrecision. It is probably better to use this rather than Rationalize to put numbers into a form that the arbitrary precision engine will use, because the latter will be manufacturing false precision.

poly3 = SetPrecision[poly[z], 15];

Plot[Evaluate[poly3], {z, 0, 1}, WorkingPrecision -> 50]


During evaluation of In[91]:= Plot::precw: The precision of the argument function <<>> is less than WorkingPrecision (50.). >>

This is an important warning because it lets you know that your results may not be valid.

See this tutorial for more information about precision. Take time to understand the difference between Mathematica's meanings of Accuracy and Precision.

## Interval

Michael's answer shows the folly of simply doing a Rationalize as I did at the start of this answer. Since questions like this come up often it would be good to have a general solution that is easily applied. I propose this rule using his formula:

machineToInterval =
c_Real?MachineNumberQ :>
Interval[ {1 - 2^-54, 1 + 2^-54} SetPrecision[c, ∞] ];


This converts any machine numbers into explicit Interval form.

To extract an ordered pair of values from a numeric Interval we merely need First.

poly3 = poly[z] /. machineToInterval;

Plot[{First @ poly3, poly[z]}, {z, 0, 1}
, WorkingPrecision -> 50
, PlotStyle -> {{Thick, Red}, ColorData[97][1]}
]


• @Wizard, Thank you. If possible can you give me some insight into the problem and the solution.
– aukie
Commented Mar 18, 2012 at 3:22
• @aukie see the updated answer Commented Mar 18, 2012 at 12:05
• To complement this, calculating the polynomial value in Horner form (HornerForm) increases precision somewhat, but not enough. Commented Mar 18, 2012 at 22:26
• @Szabolcs But strangely enough including WorkingPrecision->20 does allow HornerForm to work: Plot[HornerForm[poly[z]] // Evaluate, {z, 0, 1}, WorkingPrecision -> 20] // Quiet.
– JimB
Commented Feb 11, 2022 at 23:42

Here is another way to look at the OP's plot, which came to me after reading Daniel Lichtblau's comment under the question. His comment is worth emphasizing, which I think the following will show. (Mr.Wizard alludes to the issue in his remark on validity near the end of his answer.)

Let's say that machine precision is binary64, which has 53 bits of precision, and the maximum relative rounding error is 2^-53. Let's be a little fuzzy and say that the average error is roughly half that or 2^-54. The we can get almost exactly the envelope of the OP's plot as follows. We can replace the coefficients with ones that are slightly larger and smaller (according with the relative error); since z >= 0, the two resulting polynomials form estimates for the upper and lower bounds on poly[z]. (I'll use a high working precision, as in Mr.Wizard's answer, to compute the bounds precisely.)

Let's show two sets of bounds: bounds is in accordance with the average relative error and will form the envelope of the OP's plot; and absbounds is in accordance with the maximum relative error.

boundsPlot[error_, opts : OptionsPattern[Plot]] := With[{
polyMin = poly[z] /. c_Real :> (1 - Sign[c] error) SetPrecision[c, Infinity],
polyMax = poly[z] /. c_Real :> (1 + Sign[c] error) SetPrecision[c, Infinity]},
Plot[{polyMin, polyMax}, {z, 0, 1}, WorkingPrecision -> 50, Filling -> {1 -> {2}}, opts]
];
bounds = boundsPlot[2^-54, PlotStyle -> {Darker@Blue, Darker@Green}, PlotRange -> 5];
absbounds = boundsPlot[2^-53, PlotRange -> 5];

plotOP = Plot[poly[z], {z, 0, 1}, PlotStyle -> {AbsoluteThickness[0.7], Red}, PlotRange -> 5];

Show[plotOP, absbounds, bounds]


A good way to think of floating-point numbers is that they represent an interval of real numbers. At best a floating-point number is an approximation obtained by rounding the true value in an interval to the binary representative. An average error of 2^-54 is about as precise as the coefficients could be if they are thought of as random variables uniformly distributed throughout each interval. The bounds above show the range of likely possible values of poly[z], although occasionally poly[z] exceeds bounds.

The upshot is this: Under the hypothesis that the coefficients are known to MachinePrecision, the OP's original plot more accurately represents what is known about poly[z], than the better-looking high-precision plot. The plot bounds seems to show the likely limits of the plot of the polynomial, under the same hypothesis (although calculating the probability is beyond me). Similarly, the plot absbounds shows absolute bounds on the poly[z]. Under the hypothesis, using absbounds or bounds seems the appropriate way to represent poly[z].

If the error is greater than 2^-53, say, because the coefficients are calculated from scientific measurements with less precision, then one can compute the corresponding bounds.

• Great analysis. +1 Commented Sep 14, 2014 at 0:28
• I forgot that I needed to apply the variation to the mantissa. Correcting that I get a result that matches yours. Commented Sep 14, 2014 at 6:11
• Actually it doesn't quite, and I have a guess why. Commented Sep 14, 2014 at 6:32
• I am using r_Real :> (Interval[Rationalize[#, 0] + {-1, 1} 2^-54] 10^#2 & @@ MantissaExponent[r]). Is this more or less accurate than your method? Commented Sep 14, 2014 at 6:53
• @Mr.Wizard Ah, I see. I missed that. SetPrecision directly converts the binary representation of a machine real into a rational number with an appropriate 2^n denominator. Rationalize replaces c by the rational number with the smallest denominator within roughly c * \$MachineEpsilon of c (e.g. SetPrecision[1./3, Infinity] vs. Rationalize[1./3, 0]). SetPrecision accurately corresponds to the process I described in my answer. Rationalize does something different but close; I don't think that it is the function to use in this analysis. Commented Sep 14, 2014 at 12:30

Surprisingly enough, converting the polynomial to the Bernstein basis yields a result that Plot[] can stably deal with:

pbconv[poly_, x_] /; PolynomialQ[poly, x] := Module[{bc, n},
n = Exponent[poly, x];
bc = CoefficientList[poly, x]/Binomial[n, Range[0, n]];
Do[bc[[j + 1]] += bc[[j]], {k, n, 1, -1}, {j, k, n}];
bc]

pb[z_] = With[{n = Exponent[poly[z], z]},
N[pbconv[SetPrecision[poly[z], ∞], z]] .
BernsteinBasis[n, Range[0, n], z]]


(The cheating with SetPrecision[] is necessary because, as mentioned in the answer I linked to, the conversion procedure is inherently ill-conditioned.)

A plot:

Plot[pb[z], {z, 0, 1}]


This suggests that whatever process that generated the original polynomial might have been better reformulated in the Bernstein basis.

Settled complain about my post being "off topic". (My post got complaints about not using a special language thats too expensive for me to buy:"this is a Wolfram Language, not a mathematics, site". I did not pay attention to the brand name in the domain name. I removed the offending post for the reader's benefit.)

Here is @fischertranscripts's idea in Mathematica code (from @MichaelE2), which idea is to center the expansion at the center of the interval to reduce catastrophic cancellation:

p4 = SetPrecision[poly[z], Infinity] /. z -> y + 4/10 // Expand // N
(*
1.036149520209514 + 2.6409755026207726 y -
0.26953631057648186 y^2 + 3.449324571015004 y^3 +
338.0294181630971 y^4 - 11.745703101835632 y^5 -
43728.76625187574 y^6 - 8023.6344953736 y^7 +
2.8981359277591794*^6 y^8 + 319591.7190295984 y^9 -
1.1561838608969387*^8 y^10 + 1.6950583018096235*^7 y^11 +
2.9997779821296453*^9 y^12 - 1.5058103605142357*^9 y^13 -
5.259500639613633*^10 y^14 + 4.791402548859846*^10 y^15 +
6.297914684319244*^11 y^16 - 8.57877729466455*^11 y^17 -
5.062935150456753*^12 y^18 + 9.554822737249043*^12 y^19 +
2.5656319594205586*^13 y^20 - 6.75013786580717*^13 y^21 -
6.66863314102059*^13 y^22 + 2.924066304682213*^14 y^23 -
1.0628093808994717*^13 y^24 - 6.905409809370671*^14 y^25 +
5.294845683864764*^14 y^26 + 5.772127663081625*^14 y^27 -
9.652368445825622*^14 y^28 + 3.7404713997980006*^14 (y^29)
*)

Plot[p4, {y, -0.4, 0.6}, AxesOrigin -> {-0.4, 0},
Ticks -> {ChartingScaledTicks[{# - 0.4 &, # + 0.4 &}], Automatic}]


• Perhaps, you could plot the curve for your transformed polynomial. Commented Feb 11, 2022 at 10:27
• And show the expansion that gave that particular set of coefficients. Commented Feb 11, 2022 at 13:33
• Hard to beat Plot[Evaluate@{poly[z]}, {z, 0, 1}, WorkingPrecision -> precision], for "easy", if you're starting from Plot[{poly[z]}, {z, 0, 1}, WorkingPrecision -> precision] (assuming precision is >20). As this is a Wolfram Language, not a mathematics, site, it would be better to post your solution in WL code. For instance, p4 = SetPrecision[poly[z], Infinity] /. z -> y + 4/10 // Expand // N; Plot[p4, {y, -0.4, 0.6}], which takes 0.159 sec. compared to 0.0167 sec for precision = 20`. Commented Feb 11, 2022 at 21:37
• It was a pretty good idea though, and it happened to work pretty well. But for someone who does not know or have the tools you used, it's a lot harder than using Mathematica to accomplish the equivalent thing. -- Roll back my edit if you'd rather your answer be deleted. Commented Feb 12, 2022 at 2:19
• @fischertranscripts What I meant by "happened" is that the coefficients after the transformation happened to have negligible cancellation error. I think your transformation will reduce it, but does it always make it negligible? Commented Feb 12, 2022 at 15:53