# runge.py
# Andrew Davison, Sept. 2025, ad@coe.psu.ac.th

'''
The Runge phenomenon is a problem that arises when you try to fit a high-degree polynomial through many points, especially when the points are evenly spaced over an interval.

https://en.wikipedia.org/wiki/Runge%27s_phenomenon

Here's the key idea:

Suppose you have n data points, and you construct the unique polynomial of degree n−1 that passes exactly through them (using interpolation, e.g., Newton's method).

Intuitively, you might expect that as n increases, the 
polynomial approximation gets better.

But in reality, for some functions (especially smooth ones like f(x)=1/(1+x^2), the polynomial interpolant oscillates more and more wildly near the edges of the interval as n grows.

High-degree polynomials are very flexible, but that flexibility causes them to overshoot between points.

The oscillations are worst at the interval boundaries (a 
manifestation of Gibbs-like behavior).


Workarounds

Use fewer points (lower-degree polynomial).

Piecewise interpolation (splines): instead of one huge 
polynomial, fit small polynomials over sub-intervals.
'''

import matplotlib.pyplot as plt


def evalPoly(a, xData, x):
  # Evaluate Newton's polynomial using nested multiplication
  n = len(xData) - 1
  poly = a[n]
  for k in range(1, n + 1):
    poly = a[n - k] + (x - xData[n - k]) * poly
  return poly


def coefs(xData, yData):
  # Compute divided difference coefficients
  m = len(xData)
  a = yData[:]  # make a copy
  for k in range(1, m):
    for j in range(m-1, k-1, -1):
      a[j] = (a[j] - a[j-1]) / (xData[j] - xData[j-k])
  return a


def f(x):
  return 1/(1 + x**2)


# Test different numbers of interpolation points
Ns = [5, 11, 15]

plt.figure(figsize=(10, 6))

# True curve
xsFn = [i/50 for i in range(-250, 251)]  
   # -5 to 5 in steps of 0.02
ysFn = [f(x) for x in xsFn]
plt.plot(xsFn, ysFn, "red", ls="--", lw=2, 
             label="True function 1/(1+x^2)")

# Interpolations
for N in Ns:
  # Equally spaced sample points
  xVals = [-5 + 10*i/(N-1) for i in range(N)]
  yVals = [f(x) for x in xVals]
  
  a = coefs(xVals, yVals)
  yInterp = [evalPoly(a, xVals, x) for x in xsFn]
  
  plt.plot(xsFn, yInterp, label=f"Interpolating poly, N={N}")
  plt.scatter(xVals, yVals, s=30, marker="o")

plt.xlabel("x")
plt.ylabel("y")
plt.title("Runge Phenomenon: Interpolation of 1/(1+x^2)")
plt.ylim(-1, 2)
plt.legend()
plt.grid(True)
plt.show()