Runge's Phenomenon

Author's Name: Elena Schirò

Date: 23/12/2022

Assignment

Explain in your words what Runge's Phenomenon is, its implications in polynomial interpolation and, in general, in data over-fitting.

Demonstrate Runge's Phenomenon using Runge's function :

What is a possible remedy to avoid spurious oscillations?

Demonstrate the effectiveness of this remedy though some

Background/Motivations:

In scientific and technological studies, when dealing with sampling or measuring equipment, one obtains a limited number of data points and often wants to create a function that passes through all these points and use it to find new points in the collected data set. This operation is called interpolation. There are several methods of interpolation, including linear interpolation. This is quick and easy but not very accurate. This is why we consider linear interpolation, in which a polynomial is used to interpolate all data points. The method of Lagrange polynomials is used to find the unique interpolating polynomial of degree at most n for a given set of n+1 points. In this field of numerical analysis, however, there is a problem that arises, namely the Runge phenomenon.

Problem statement:

The Runge phenomenon concerns the field of numerical analysis and is a problem that occurs in the context of polynomial interpolation, specifically Lagrange interpolation. This consists of an increase in the amplitude of the error near the extremes of the interval. Considering the function:

The mathematician Carl David Tolmè Runge discovered that by interpolating this function in a set of equidistant x_i points in the interval [-1;1], with a polynomial P_n (x) of degree ≤n, the resulting interpolation oscillates in amplitude towards the extremes of the interval. This consists of the amplitude increase of the error near the extremes of the interval.

from __future__ import division from numpy import (zeros_like, pi, cos, zeros, amax, abs, array, linspace, prod) import matplotlib.pyplot as plt def lagrange(x_int, y_int, x_new): y_new = zeros_like(x_new) for xi, yi in zip(x_int, y_int): y_new += yi*prod([(x_new - xj)/(xi - xj) for xj in x_int if xi!=xj], axis=0) return y_new def gauss_lobatto(N, tol=1e-15): x = -cos(linspace(0, pi, N)) P = zeros((N, N)) # Vandermonde Matrix x_old = 2 while amax(abs(x - x_old)) > tol: x_old = x P[:, 0] = 1 P[:, 1] = x for k in range(2, N): P[:, k] = ((2 * k - 1) * x * P[:, k - 1] - (k - 1) * P[:, k - 2]) / k x = x_old - (x * P[:, N - 1] - P[:, N - 2]) / (N * P[:, N - 1]) print(x) return array(x) runge = lambda x: 1/(1 + 25*x**2) x = linspace(-1, 1, 100) x_int = linspace(-1, 1, 11) x_int2 = gauss_lobatto(11) x_new = linspace(-1, 1, 1000) y_new = lagrange(x_int, runge(x_int), x_new) y_new2 = lagrange(x_int2, runge(x_int2), x_new) plt.plot(x, runge(x), "k") plt.plot(x_new, y_new) plt.plot(x_new, y_new2) plt.legend(["Runge function", "Uniform interpolation", "Lobatto-sampling interpolation"]) plt.xlabel("x") plt.ylabel("y") plt.show()

[-1.         -0.93199267 -0.7828368  -0.56457705 -0.29566305  0.
  0.29566305  0.56457705  0.7828368   0.93199267  1.        ]
[-1.         -0.93400371 -0.7844839  -0.56523534 -0.29575814  0.
  0.29575814  0.56523534  0.7844839   0.93400371  1.        ]
[-1.         -0.93400143 -0.78448347 -0.56523533 -0.29575814  0.
  0.29575814  0.56523533  0.78448347  0.93400143  1.        ]
[-1.         -0.93400143 -0.78448347 -0.56523533 -0.29575814  0.
  0.29575814  0.56523533  0.78448347  0.93400143  1.        ]
[-1.         -0.93400143 -0.78448347 -0.56523533 -0.29575814  0.
  0.29575814  0.56523533  0.78448347  0.93400143  1.        ]

The error expression of the Lagrange polynomial is:

It is wrong to think that increasing n, that is the number of nodes, decreases the error because often increasing n increases the complexity of the derivative (n+1) and in addition this causes the Runge phenomenon.

Methods:

To solve the problem of Runge's phenomenon, one can change the polynomial representing the error to one that minimises this oscillatory trend and a better result can be found. This is done by taking the zeros of the Chebyshev polynomial as nodes. In mathematics, the Čebyšëv knots, Čebyšëv-Gauss-Lobatto knots, or Čebyšëv roots, are the roots of the Čebyšëv polynomials. For each natural integer n, the n-th polynomial has n simple roots within the interval [-1;1]. Such an n-uple is a good choice for an interpolation on n points in the said interval, as it allows for an a priori increase in the interpolation error.

The Čebyšëv nodes of the n-th polynomial are given by

Results:

We have Tn that is the n-th polynomial of Čebyšëv:

The cosine function has periodic roots:

for each integer i giving:

Then the roots of the n-th Čebyšëv polynomial are found when:

which can be solved for xi obtaining:

Conclusions:

At the end of this analysis we can say that adding more terms to regular polynomials is useful for some functions, such as the sine, but fails with others, such as the Runge function. This can be seen in the following graphs:

import numpy as np import matplotlib.pyplot as plt def runge(x): """In some places the x range is expanded and the formula give as 1/(1+x^2) """ return 1 / (1 + 25 * x ** 2) def polyfit(x, y): """Return a polynomial fit of order n+1 to n points""" z = np.polyfit(x, y, x.size + 1) return np.poly1d(z) def cheby_nodes(x0, x1, n): return (x1 - x0) * (np.cos((2*np.arange(1, n + 1) - 1)/(2*n)*np.pi) + 1) / 2 + x0 def overlay_plot_fit(f, f_fit, x, y, x0=-1, x1=1): """x, y are the original points fitted to """ _x = np.linspace(x0, x1, 500) plt.plot(_x, f_fit(_x), label='polyfit') plt.plot(_x, f(_x), label='actual') plt.plot(x, y, 'k.') plt.title(x.size) plt.legend(loc='best', fontsize=9) def compare_fit(f, x, x0=-3, x1=3): y = f(x) f_fit = polyfit(x, y) overlay_plot_fit(f, f_fit, x, y, x0=x0, x1=x1) def compare_fit_panel(f, x0=-3, x1=3, zoom=False): n_list = [2, 4, 8, 16, 32, 64] for n in range(6): # x = np.linspace(x0, x1, n_list[n]) x = cheby_nodes(x0, x1, n_list[n]) ax = plt.subplot(2, 3, n + 1) compare_fit(f, x) if n != 3: plt.setp(ax, xticklabels=[], yticklabels=[]) if zoom: plt.setp(ax, ylim=[-1, 1]) plt.suptitle(f) plt.figure(figsize=(10, 6)) compare_fit_panel(np.sin) plt.figure(figsize=(10, 6)) compare_fit_panel(runge) plt.figure(figsize=(10, 6)) compare_fit_panel(runge, zoom=True)

/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)
/tmp/ipykernel_89/2993306733.py:35: RankWarning: Polyfit may be poorly conditioned
  f_fit = polyfit(x, y)

Čebyšëv's method proves to be the most effective in solving the problem of the runge phenomenon that occurs in the case of polynomial interpolation.

Bibliography:

https://www.wikiwand.com/it/Nodi_di_%C4%8Ceby%C5%A1%C3%ABv https://nicoguaro.github.io/posts/numerical_challenge/numerical-10/ https://gist.github.com/kghose/c6580def53156c92eb8e1025cc4c7444 http://www.mat.unimi.it/users/alzati/Geometria_Computazionale_98-99/apps/interpolanti/teoria.html https://it.wikipedia.org/wiki/Fenomeno_di_Runge https://www.wikiwand.com/it/Nodi_di_%C4%8Ceby%C5%A1%C3%ABv https://tex.unica.it/~gppe/did/ca/tesine/2003/03gallo.pdf

.css-ytumd6{-webkit-text-decoration:none;text-decoration:none;}Runge's Phenomenon