Runge's phenomenon

Lagrange's interpolation appears in every first course in numerical analysis. The proof of the interpolation error formula is beautiful and it is useful to introduce several algorithms. In my opinion in some of these courses Runge's phenomenon, showing the failure of a naive approach to interpolation, it is not sufficiently emphasized.

The plots

We always consider n+1 nodes -1=x_0<...<x_n=1 evenly distributed in [-1,1].

A nice function like f(x)=sin(pi*x) show that even with few nodes the approximation is surprisingly good. This is the graph for 7 nodes. The interpolation polynomial (in blue) barely show up under the function (in red).

If we increase the number of nodes, nothing strange happens.

If we now try f(x)=1/(8x^2+1), that is also a nice infinitely differentiable function, the results are strange at the extremes and they worsen exponentially with the number of nodes.

7 nodes

14 nodes

21 nodes

28 nodes

A first guess is that this is related to the size of the derivative. Let us try f(x)=sin(5*pi*x) that has a derivative much larger that the previous example (even the first function is a counterexample). The first experiments seem to confirm the guess but there is a big difference: In this case the strange peaks eventually disappear while for f(x)=1/(8x^2+1), as we mentioned before, they grow exponentially.

7 nodes

14 nodes

21 nodes

28 nodes

One can explain in a satisfactory theoretical way these examples using interpolation error formula. The function f(x)=1/(8x^2+1) has a finite radius of convergence that forces the successive derivatives to grow exponentially. On the other hand, sine is an entire function.

The code

The plots were done with the following SAGE code. I know that it is highly inefficient. Even worse, I took the three lines that probably grind your gears from somewhere on the internet. My idea was to write quickly a brief program easy to understand for students.

#######################

#### RUNGE

#######################

      

      

def rung(N):

    #f = sin(5*pi*x)

    #f = sin(pi*x)

    #f = 1/(8*x^2+1)

    h = 2/(N-1)

    xs = srange(-1,1+h/2,h)

    ys = [ f(x=xx).n() for xx in xs]

      

    p_i = prod(x-ii for ii in xs); 

    pid = diff(p_i,x);

    g = sum(p_i/(x-ii)/pid(x=ii)*jj for (ii,jj) in zip(xs,ys))

      

    P = plot(g,x,-1, 1, thickness=3)

    P += plot(f,x,-1, 1, color='red', thickness=3)

    P += list_plot(zip(xs,ys), size=80, zorder=200)

      

    #P.set_aspect_ratio(1)

    #P.axes(False)

    P.fontsize(25)

    return P

      

P = rung(20)

      

P.show()