69 lines
2.0 KiB
Python
69 lines
2.0 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
def runge(x):
|
|
return 1 / (1 + 25 * x**2)
|
|
|
|
def divided_differences(x, y):
|
|
n = len(y)
|
|
a = np.zeros((n,n))
|
|
|
|
a[:,0] = y
|
|
|
|
for i in range(1,n):
|
|
for j in range(1,i+1):
|
|
a[i,j] = (a[i,j-1] - a[i-1,j-1]) / (x[i] - x[i-j]) # (links daneben - links darüber) / (x[zeile] - x[zeile - spalte])
|
|
|
|
return a
|
|
|
|
def newton_interpolation(a, data, x):
|
|
n = len(a)
|
|
|
|
p = np.zeros(len(x))
|
|
|
|
for i in range(1,n+1):
|
|
p = a[n-i] + (x - data[n-i]) * p # Horner Schema; '-' and '*' are overloaded for numpy arrays: so at least one argument 'data' or 'x' has to be numpy array
|
|
|
|
return p
|
|
|
|
####################################################################################################################################################################
|
|
n = 12
|
|
|
|
x = np.linspace(-1, 1, 200)
|
|
x_e = np.linspace(-1, 1, n) # equidistant grid points
|
|
|
|
# Chebyshev grid points
|
|
x_c = np.zeros(n)
|
|
for i in range(0,n):
|
|
x_c[i] = np.cos((2 * i + 1) * np.pi / (2 * n))
|
|
|
|
f = runge(x)
|
|
y_e = runge(x_e) # values for grid points for interpolation with equidistant grid points
|
|
y_c = runge(x_c) # values for grid points for interpolation with Chebyshev grid points
|
|
|
|
# Interpolation with equidistant grid points and evaluation of interpolated values at x
|
|
a_e = np.diag(divided_differences(x_e, y_e))
|
|
p_e = newton_interpolation(a_e, x_e, x)
|
|
|
|
# Interpolation with Chebyshev grid points and evaluation of interpolated values at x
|
|
a_c = np.diag(divided_differences(x_c, y_c))
|
|
p_c = newton_interpolation(a_c, x_c, x)
|
|
|
|
|
|
# Plotting of Runge and the two interpolated polynomials
|
|
plt.rcParams['text.usetex'] = True
|
|
|
|
plt.plot(x, f, label=r"$f(x) = \frac{1}{1 + 25 x^2}$")
|
|
plt.plot(x, p_e, label=r"$p_n(x)$ equidistant")
|
|
plt.plot(x, p_c, label=r"$p_n(x)$ Chebyshev")
|
|
plt.plot([], [], ' ', label="$n = {}$".format(n))
|
|
|
|
plt.title("A3)")
|
|
plt.xlabel("x")
|
|
plt.ylabel("y")
|
|
plt.legend()
|
|
|
|
plt.show()
|
|
|
|
# braendel@math.tu-freiberg.de
|