From bd7e39ca9ceda95001fa1e4c73b0a7700952045d Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@gmail.com>
Date: Wed, 12 Feb 2020 12:30:56 +0100
Subject: [PATCH] Finished Exercise 3

---
 week1/ex3.py | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)
 create mode 100755 week1/ex3.py

diff --git a/week1/ex3.py b/week1/ex3.py
new file mode 100755
index 0000000..3f2ccb6
--- /dev/null
+++ b/week1/ex3.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+
+# Runge's Phenomenon
+#----------------------------------
+
+import numpy as np
+
+def func_to_interpolate(x):
+    return 1/( 1 + 25 * x**2 )
+
+def equidistant_nodes(n):
+    return 2*np.arange(0,n+1)/n -1
+
+def chebychev_nodes(n):
+    return np.cos( (2*np.arange(1,n+1) -1) * np.pi / (2*n) )
+
+def random_nodes(n):
+    return 2*np.random.rand(n) -1
+
+def main():
+    import matplotlib.pyplot as plt
+    from scipy.interpolate import lagrange as scipy_lagrange
+    
+    np.random.seed(0)
+
+
+    n_special = 20
+    n_cheby = n_special
+    n_equi = n_special
+    n_random = n_special
+
+    n_func = 1e4
+    x_func = np.linspace(-1, 1, n_func, endpoint = True )
+
+    fig, ax = plt.subplots()
+
+    ax.set_title("Lagrange Interpolations with different types of nodes")
+    ax.grid()
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+
+    func = func_to_interpolate
+    # Plot the formula
+    ax.plot(x_func, func(x_func), '--', label="Function $1/(1+25x^2)$")
+
+    # Using equidistant nodes
+    x_equi = equidistant_nodes(n_equi)
+    ax.plot(x_equi, scipy_lagrange(x_equi, func(x_equi))(x_equi), '^', label="Equidistant (n={})".format(n_equi))
+
+    # Using chebychev nodes
+    x_cheby = chebychev_nodes(n_cheby)
+    ax.plot(x_cheby, scipy_lagrange(x_cheby, func(x_cheby))(x_cheby), '>',label="Chebyshev (n={})".format(n_cheby))
+
+    # Using random nodes
+    x_random = random_nodes(n_random)
+    ax.plot(x_random, scipy_lagrange(x_random, func(x_random))(x_random), '.', label="Random (n={})".format(n_random))
+
+
+    ax.legend()
+
+
+    print("""
+    Runge's Phenomenon
+    ------------------
+
+    According to the docs [1], the scipy Lagrange
+    interpolation is unstable and there should not be more
+    than 20 points.
+
+
+    The region where the equidistant nodes break is near
+    x = -1, where the interpolation is fluctuating between
+    large negative and large positive values. From n > 25,
+    this starts to get visible.
+
+    The chebyshev nodes break near x = 1 when n > 30, where
+    it interpolates the function to be negative.
+
+    Funnily enough, the equidistant nodes vary much more
+    than the chebyshev nodes.
+
+    For the random nodes, you just have be lucky that not
+    too many points are close together at the extremes of
+    the interval.
+
+    [1] = https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.lagrange.html
+""")
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()