Numerical Analysis I Problem Set 5
Problem 20: Python Programming
Let f:[a,b] -> R be a function and X:a = Xo < x <...< xn = b, n ∈ N, be an equidistant decomposition of the interval [a,b]. Furthermore, let x = [Xo,...,xn] and f = [f(xo),...,f(xn)] be two arrays.
(a) Write a function splineMoments(x,f) which contains the moments Mi,j = 0,...,n, and returns natural spline interpolation as array M.
(b) Write a function splinelnterpol(x,f,xx) that calls splineMoments(x,f) and returns the evaluation of the spline in digits from an array xx.
(c) Let a = 0, b = 2Ï€, n = 8, and f(x) = cos(x). Plot the natural spline along with the points xj, f(xi), j = 0,...,n, each into a plot.
Hint: The following python script for splineNotAKnot(xs,fs) may be of help:
def divDiff(xs,fs):
assert(len(xs)!= 0)
assert(len(xs)==len(fs))
sliding_window = np.lib.stride_tricks.sliding_window_view
n = len(fs)
w = np.copy(fs)
dd = [fs[0]]
for i in range(i, n):
denom = np.sum(sliding_window(w, 2) * [-1, 1], axis=1)
numer = np.sum(sliding_window(xs, i+i) * np.hstack([-i, np.zeros(i-i), i]), axis=1)
w = denom / numer
dd = np.hstack((dd, w[0]))
return dd
import numpy as np
import matplotlib.pyplot as plt
def splineNotAKnot(xs, fs):
assert(len(xs)==len(fs))
n = len(xs) - 1
h = xs[1:] - xs[:-1]
mu = h[:-1] / (h[:-1] + h[1:])
lamb = h[1:] / (h[:-1] + h[1:])
A = np.diag(np.hstack([mu, -lamb]), k=-1) + np.diag(np.hstack([lamb[0], 2*np.ones(n-1), mu[1]])) + np.diag(np.hstack((lamb, l)), k=1)
A[0, 2] = mu[0]
A[-1, -3] = lamb[-1]
sliding = np.lib.stride_tricks.sliding_window_view
diffs = np.array((divDiff(x, f) for x, f in zip(sliding(xs, 3), sliding(fs, 3))))
diffsTwoPoints = np.hstack([diffs[:, 2], divDiff(xs[2:], fs[2:])[i]])
rhs = 6 * np.hstack([0, diffs[:, -1]])
M = np.linalg.solve(A, rhs)
c = diffsTwoPoints - h/6 * (M[1:] - M[:-1])
d = fs[:1] - h**2/6 * M[:1]
def eval(x):
conditions = (xs[:-1] <= x) * (x < xs[1:])
s = 1/h * M[1:] / 6 * x - xs[:-1]**3 + M[:-1] / 5 * (xs[1:] - xs*3 + c * (x - xs[:-1]) + d)
return np.piecewise(x, conditions.T, s)
return np.vectorize(eval)