pgfplot renders jagged polynomial

Question

I keep getting a very jagged plot of this 9th degree polynomial, regardless of using smooth or increasing the sample points.

Ideas?

\documentclass{standalone}

\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\usetikzlibrary{fpu}

\begin{document}
\begin{tikzpicture}[scale=.55]
\begin{axis}[
xlabel={$x$}, ylabel={$y$}
,restrict y to domain=-10^9:10^9
,axis lines=middle
,samples=2000,
%,grid
,thick
,domain=0:4
,xtick={0,1,2,3,4}
,ytick={0,9}
,xmin=0
,xmax=4.5
,ymin=-1
,ymax=1.8
,xlabel shift={1in}
,y tick label style={yshift={5pt}}
%,legend pos=outer north east
]
\addplot+[black,no marks,smooth]{1.5*x + 3.83333*x^2 - 15.0104*x^3 + 22.7454*x^4 - 19.2028*x^5 +
 9.37731*x^6 - 2.5978*x^7 + 0.377315*x^8 - 0.0222801*x^9};
\end{axis}
\end{tikzpicture}
\end{document}

If you use gnuplot, \addplot+[black,no marks,smooth] gnuplot {1.5*x + 3.83333*x^2 - 15.0104*x^3 + 22.7454*x^4 - 19.2028*x^5 + 9.37731*x^6 - 2.5978*x^7 + 0.377315*x^8 - 0.0222801*x^9};, the problem disappears. You need to run with -shell-escape and have gnuplot installed. — , Nov 24 '19 at 04:07
I installed gnuplot, but how do I do the '-shell-escape'? I am using miktex on windows 8 rather than compiling my latex on the command line. — JohnD, Nov 24 '19 at 04:41
Unfortunately I do not know that but I would be surprised if no one had asked this. (I cannot test proposals because I do not use MikTeX.) — , Nov 24 '19 at 04:45
Shell-escape on miktex: https://tex.stackexchange.com/questions/37489/how-can-i-enable-write-18-on-a-miktex-installation — Rmano, Nov 24 '19 at 11:48
One can gain some improvement using ((((((((-0.0222801*x + 0.377315)*x - 2.5978)*x + 9.37731)*x - 19.2028)*x + 22.7454)*x - 15.0104)*x + 3.83333)*x +1.5)*x, but the round off errors will still bite you in the end. I do all my calculations using C++ and write the values into a table for pgfplots. — John Kormylo, Nov 24 '19 at 17:49
I said windows 8 above, but I meant windows 10. I still can't get even a basic gnuplot to compile with my miktex set up even using -shell-escape or -enable-write18. It seems miktex is not interacting with gnuplot properly. — JohnD, Nov 24 '19 at 17:51
@JohnKormylo --- second order segments are another possibility... ;-) — Rmano, Nov 24 '19 at 19:38
@Schrödinger'scat you can avoid gnuplot if you do not stress poor TeX math with x^9 ;-) — Rmano, Nov 24 '19 at 19:39
@Rmano Yes, that's a good point. gnuplot offers a solution that almost always works out of the box, but of course there are reasons for not wanting to depend on external software. — , Nov 24 '19 at 19:43

Rmano · Answer 1 · 2022-06-12T09:20:32.553

As commented by John Kormylo, the problem here is precision --- adding samples will not help.

My proposal is to split your polynomial into second-order-section terms; that will avoid the problematic high powers (like x^9). I have written a small SymPy thingy to do this and results that your polynomial is equivalent to:

-0.0222801*x*(x - 4.04307547777127)*(x + 0.193642687986076)*(x**2 - 7.81990533113625*x + 15.3415387466162)*(x**2 - 4.09164256641596*x + 4.41554633036612)*(x**2 - 1.1740869290557*x + 1.26942721414423)

Which is normally a much better-behaved function to compute. You can see it:

\documentclass[border=10pt]{standalone}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\usetikzlibrary{fpu}
\begin{document}
\begin{tikzpicture}[scale=.55]
\begin{axis}[
xlabel={$x$}, ylabel={$y$}
,restrict y to domain=-10^9:10^9
,axis lines=middle
,samples=200,
%,grid
,thick
,domain=0:4
,xtick={0,1,2,3,4}
,ytick={0,9}
,xmin=0
,xmax=4.5
,ymin=-1
,ymax=1.8
,xlabel shift={1in}
,y tick label style={yshift={5pt}}
%,legend pos=outer north east
]
% \addplot+[black,no marks,]{1.5x + 3.83333x^2 - 15.0104x^3 + 22.7454x^4 - 19.2028x^5 +
 % 9.37731x^6 - 2.5978x^7 + 0.377315x^8 - 0.0222801x^9};
\addplot+[red,no marks,]{-0.0222801x(x - 4.04307547777127)(x + 0.193642687986076)
    (xx - 7.81990533113625x + 15.3415387466162)
    (xx - 4.09164256641596x + 4.41554633036612)
    (xx - 1.1740869290557*x + 1.26942721414423)};
\end{axis}
\end{tikzpicture}
\end{document}

(I can share the python code if you want --- but it's just a hack, and probably everybody out here can write something more elegant).

SymPi code to obtain the second-order-terms form of the polynomial

(After request, here is the python code); notice that I left the mathjax code for the formula, because the PNG are really small...

import sys
import os
# BEWARE only for command line interactive usage
from math import *
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
sys.path.append("/home/romano/lib/python")
import sympy as sym
from sympy import init_printing
init_printing()
sym.__version__

'1.1.1'

x, y, z = sym.symbols("x y z")

Input the polynomial here

Y=1.5*x + 3.83333*x**2 - 15.0104*x**3 + 22.7454*x**4 - 19.2028*x**5 +9.37731*x**6 - 2.5978*x**7 + 0.377315*x**8 - 0.0222801*x**9

P=sym.Poly(Y); P

$$\operatorname{Poly}{\left( - 0.0222801 x^{9} + 0.377315 x^{8} - 2.5978 x^{7} + 9.37731 x^{6} - 19.2028 x^{5} + 22.7454 x^{4} - 15.0104 x^{3} + 3.83333 x^{2} + 1.5 x, x, domain=\mathbb{R} \right)}$$

Find roots of the polynomial

rs=sym.polys.polyroots.roots(P, multiple=True); rs

$$\left [ -0.193642687986076, \quad 0, \quad 4.04307547777127, \quad 0.587043464527849 - 0.961668958061643 i, \quad 0.587043464527849 + 0.961668958061643 i, \quad 2.04582128320798 - 0.479751610251982 i, \quad 2.04582128320798 + 0.479751610251982 i, \quad 3.90995266556812 - 0.23196745382247 i, \quad 3.90995266556812 + 0.23196745382247 i\right ]$$

Dirty hack to obtain the sos (second order segments) form. It relies on the fact that the order of the roots is fixed, and I am not sure if it's true: the complex conjugate terms must come one after the other.

tot=1
skip = False
for r in rs:
    rr = complex(r)
    if rr.imag==0:
        term = (x-rr)
        tot = tot*term
    else:
        if skip==False: #first complex root 
            skip=True
            a = rr.real
            b = rr.imag
            term = (x*x-2*a*x+(a*a+b*b))
            tot = term * tot
        else:
            skip=False
tot

$$x \left(x - 4.04307547777127\right) \left(x + 0.193642687986076\right) \left(x^{2} - 7.81990533113625 x + 15.3415387466162\right) \left(x^{2} - 4.09164256641596 x + 4.41554633036612\right) \left(x^{2} - 1.1740869290557 x + 1.26942721414423\right)$$

f=sym.LC(P); print(f); print (f*tot)

-0.0222801000000000
-0.0222801*x*(x - 4.04307547777127)*(x + 0.193642687986076)*(x**2 - 7.81990533113625*x + 15.3415387466162)*(x**2 - 4.09164256641596*x + 4.41554633036612)*(x**2 - 1.1740869290557*x + 1.26942721414423)

check that the two things are the same... you never know

Z=sym.expand(f*tot); Z

$$- 0.0222801 x^{9} + 0.377315 x^{8} - 2.5978 x^{7} + 9.37731 x^{6} - 19.2028 x^{5} + 22.7454 x^{4} - 15.0104 x^{3} + 3.83333 x^{2} + 1.5 x$$

Could you share the python code for approximating the curve? — antshar, Jun 11 '22 at 08:34
@antshar I'll see if I can find it again, but basically I just looked for the roots of then original polinomial and the joined the complex conjugate ones. — Rmano, Jun 11 '22 at 15:49