Cubic splines Tobin A. Driscoll Richard J. Braun
A piecewise linear interpolant is continuous but has discontinuities in its derivative. We often desire a smoother interpolant, i.e., one that has some continuous derivatives. By far the most popular choice is piecewise cubic.
We use S ( x ) S(x) S ( x ) t 0 < t 1 < ⋯ < t n t_0 < t_1 < \cdots < t_n t 0 < t 1 < ⋯ < t n y 0 , … , y n y_0,\ldots,y_n y 0 , … , y n k = 1 , … , n k=1,\ldots,n k = 1 , … , n S ( x ) S(x) S ( x ) [ t k − 1 , t k ] [t_{k-1},t_k] [ t k − 1 , t k ] S k ( x ) S_k(x) S k ( x )
S k ( x ) = a k + b k ( x − t k − 1 ) + c k ( x − t k − 1 ) 2 + d k ( x − t k − 1 ) 3 , k = 1 , … , n , S_k(x) = a_k + b_k(x-t_{k-1}) + c_k(x-t_{k-1})^2 + d_k(x-t_{k-1})^3, \qquad k=1,\ldots,n, S k ( x ) = a k + b k ( x − t k − 1 ) + c k ( x − t k − 1 ) 2 + d k ( x − t k − 1 ) 3 , k = 1 , … , n , where a k , b k , c k , d k a_k,b_k,c_k,d_k a k , b k , c k , d k 4 n 4n 4 n
5.3.1 Smoothness conditions ¶ We are able to ensure that S S S
1. Interpolation by S k S_k S k
Algebraically we require S k ( t k − 1 ) = y k − 1 S_k(t_{k-1})=y_{k-1} S k ( t k − 1 ) = y k − 1 S k ( t k ) = y k S_k(t_k)=y_k S k ( t k ) = y k k = 1 , … , n k=1,\dots,n k = 1 , … , n (5.3.1)
a k = y k − 1 , a_k = y_{k-1}, a k = y k − 1 , a k + b k h k + c k h k 2 + d k h k 3 = y k , k = 1 , … , n , a_k + b_k h_k + c_k h_k^2 + d_k h_k^3 = y_{k}, \qquad k=1,\ldots,n, a k + b k h k + c k h k 2 + d k h k 3 = y k , k = 1 , … , n , where we have used the definition
h k = t k − t k − 1 , k = 1 , … , n . h_k = t_{k}-t_{k-1}, \qquad k=1,\ldots,n. h k = t k − t k − 1 , k = 1 , … , n . The values of h k h_k h k linearly in the constraint equations. So we will express the constraints using linear algebra. The left endpoint interpolation constraints (5.3.2)
[ I 0 0 0 ] [ a b c d ] = [ y 0 ⋮ y n − 1 ] , \begin{bmatrix}
\mathbf{I} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0}
\end{bmatrix}
\begin{bmatrix}
\mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d}
\end{bmatrix}
=
\begin{bmatrix}
y_0 \\ \vdots \\ y_{n-1}
\end{bmatrix}, [ I 0 0 0 ] ⎣ ⎡ a b c d ⎦ ⎤ = ⎣ ⎡ y 0 ⋮ y n − 1 ⎦ ⎤ , with I \mathbf{I} I n × n n\times n n × n (5.3.3)
[ I H H 2 H 3 ] [ a b c d ] = [ y 1 ⋮ y n ] , \begin{bmatrix}
\mathbf{I} & \mathbf{H} & \mathbf{H}^2 & \mathbf{H}^3
\end{bmatrix}
\begin{bmatrix}
\mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d}
\end{bmatrix}
=
\begin{bmatrix}
y_1 \\ \vdots \\ y_{n}
\end{bmatrix}, [ I H H 2 H 3 ] ⎣ ⎡ a b c d ⎦ ⎤ = ⎣ ⎡ y 1 ⋮ y n ⎦ ⎤ , where we have defined the diagonal matrix
H = [ h 1 h 2 ⋱ h n ] . \mathbf{H} =
\begin{bmatrix}
h_1 & & & \\ & h_2 & & \\ & & \ddots & \\ & & & h_n
\end{bmatrix}. H = ⎣ ⎡ h 1 h 2 ⋱ h n ⎦ ⎤ . Collectively, (5.3.5) (5.3.6) 2 n 2n 2 n
2. Continuity of S ′ ( x ) S'(x) S ′ ( x )
We do not know what the slope of the interpolant should be at the nodes, but we do want the same slope whether a node is approached from the left or the right. Thus, we obtain constraints at the nodes that sit between two neighboring piecewise definitions, so that S 1 ′ ( t 1 ) = S 2 ′ ( t 1 ) S_1'(t_1)=S_2'(t_1) S 1 ′ ( t 1 ) = S 2 ′ ( t 1 )
b k + 2 c k h k + 3 d k h k 2 = b k + 1 , k = 1 , … , n − 1. b_k + 2 c_k h_k + 3 d_k h_k^2 = b_{k+1}, \qquad k=1,\dots,n-1. b k + 2 c k h k + 3 d k h k 2 = b k + 1 , k = 1 , … , n − 1. Moving the unknowns to the left side, as a system these become
E [ 0 J 2 H 3 H 2 ] [ a b c d ] = 0 , \mathbf{E}
\begin{bmatrix}
\boldsymbol{0} & \mathbf{J} & 2\mathbf{H} & 3\mathbf{H}^2
\end{bmatrix}
\begin{bmatrix}
\mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d}
\end{bmatrix}
= \boldsymbol{0}, E [ 0 J 2 H 3 H 2 ] ⎣ ⎡ a b c d ⎦ ⎤ = 0 , where now we have defined
J = [ 1 − 1 1 − 1 ⋱ ⋱ 1 − 1 1 ] , \mathbf{J} =
\begin{bmatrix}
1 & -1 & & & \\ & 1 & -1 & & \\ & & \ddots & \ddots & \\ & & &1 & -1 \\ & & & & 1
\end{bmatrix}, J = ⎣ ⎡ 1 − 1 1 − 1 ⋱ ⋱ 1 − 1 1 ⎦ ⎤ , and E \mathbf{E} E ( n − 1 ) × n (n-1)\times n ( n − 1 ) × n
E = [ 1 0 1 0 ⋱ ⋱ 1 0 ] . \mathbf{E} =
\begin{bmatrix}
1 & 0 & & & \\ & 1 & 0 & & \\ & & \ddots & \ddots & \\ & & & 1& 0
\end{bmatrix}. E = ⎣ ⎡ 1 0 1 0 ⋱ ⋱ 1 0 ⎦ ⎤ . Left-multiplying by E \mathbf{E} E (5.3.9) n − 1 n-1 n − 1 n − 1 n-1 n − 1
3. Continuity of S ′ ′ ( x ) S''(x) S ′′ ( x )
These again apply only at the interior nodes t 1 , … , t n − 1 t_1,\dots,t_{n-1} t 1 , … , t n − 1 S 1 ′ ′ ( t 1 ) = S 2 ′ ′ ( t 1 ) S_1''(t_1)=S_2''(t_1) S 1 ′′ ( t 1 ) = S 2 ′′ ( t 1 ) (5.3.1)
2 c k + 6 d k h k = 2 c k + 1 , k = 1 , … , n − 1. 2 c_k + 6 d_k h_k = 2c_{k+1}, \qquad k=1,\dots,n-1. 2 c k + 6 d k h k = 2 c k + 1 , k = 1 , … , n − 1. In system form (after canceling a factor of 2 from each side) we get
E [ 0 0 J 3 H ] [ a b c d ] = 0 . \mathbf{E}
\begin{bmatrix}
\boldsymbol{0} & \boldsymbol{0} & \mathbf{J} & 3\mathbf{H}
\end{bmatrix}
\begin{bmatrix}
\mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d}
\end{bmatrix}
= \boldsymbol{0}. E [ 0 0 J 3 H ] ⎣ ⎡ a b c d ⎦ ⎤ = 0 . 5.3.2 End conditions ¶ So far the equations (5.3.5) (5.3.6) (5.3.9) (5.3.13) 2 n + ( n − 1 ) + ( n − 1 ) = 4 n − 2 2n+(n-1)+(n-1)=4n-2 2 n + ( n − 1 ) + ( n − 1 ) = 4 n − 2 4 n 4n 4 n (5.3.1) S ′ S' S ′ S ′ ′ S'' S ′′
Natural spline: S 1 ′ ′ ( t 0 ) = S n ′ ′ ( t n ) = 0 \quad S_1''(t_0)=S_n''(t_n)=0 S 1 ′′ ( t 0 ) = S n ′′ ( t n ) = 0 Not-a-knot spline: S 1 ′ ′ ′ ( t 1 ) = S 2 ′ ′ ′ ( t 1 ) , S n − 1 ′ ′ ′ ( t n − 1 ) = S n ′ ′ ′ ( t n − 1 ) \quad S_1'''(t_1)=S_2'''(t_1), \; S_{n-1}'''(t_{n-1})=S_n'''(t_{n-1}) S 1 ′′′ ( t 1 ) = S 2 ′′′ ( t 1 ) , S n − 1 ′′′ ( t n − 1 ) = S n ′′′ ( t n − 1 ) While natural splines have important theoretical properties, not-a-knot splines give better pointwise accuracy, and they are the only type we consider further.
In the not-a-knot spline, the values and first three derivatives of the cubic polynomials S 1 S_1 S 1 S 2 S_2 S 2 t 1 t_1 t 1 S n − 1 S_{n-1} S n − 1 S n S_n S n [1]
d 1 = d 2 , d n − 1 = d n . d_1=d_2, \quad d_{n-1}=d_n. d 1 = d 2 , d n − 1 = d n . Collectively, (5.3.5) (5.3.6) (5.3.9) (5.3.13) (5.3.14) 4 n 4n 4 n (5.3.1) cardinal basis, so the solution of a linear system cannot be avoided.
5.3.3 Implementation ¶ 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
"""
spinterp(t, y)
Construct a cubic not-a-knot spline interpolating function for data
values in `y` given at nodes in `t`.
"""
function spinterp(t, y)
n = length(t) - 1
h = [t[k+1] - t[k] for k in 1:n]
# Preliminary definitions.
Z = zeros(n, n)
In = I(n)
E = In[1:n-1, :]
J = diagm(0 => ones(n), 1 => -ones(n - 1))
H = diagm(0 => h)
# Left endpoint interpolation:
AL = [In Z Z Z]
vL = y[1:n]
# Right endpoint interpolation:
AR = [In H H^2 H^3]
vR = y[2:n+1]
# Continuity of first derivative:
A1 = E * [Z J 2 * H 3 * H^2]
v1 = zeros(n - 1)
# Continuity of second derivative:
A2 = E * [Z Z J 3 * H]
v2 = zeros(n - 1)
# Not-a-knot conditions:
nakL = [zeros(1, 3 * n) [1 -1 zeros(1, n - 2)]]
nakR = [zeros(1, 3 * n) [zeros(1, n - 2) 1 -1]]
# Assemble and solve the full system.
A = [AL; AR; A1; A2; nakL; nakR]
v = [vL; vR; v1; v2; 0; 0]
z = A \ v
# Break the coefficients into separate vectors.
rows = 1:n
a = z[rows]
b = z[n.+rows]
c = z[2*n.+rows]
d = z[3*n.+rows]
S = [Polynomial([a[k], b[k], c[k], d[k]]) for k in 1:n]
# This function evaluates the spline when called with a value
# for x.
return function (x)
if x < t[1] || x > t[n+1] # outside the interval
return NaN
elseif x == t[1]
return y[1]
else
k = findlast(x .> t) # last node to the left of x
return S[k](x - t[k])
end
end
end1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
function S = spinterp(t, y)
% SPINTERP Cubic not-a-knot spline interpolation.
% Input:
% t interpolation nodes (vector, length n+1)
% y interpolation values (vector, length n+1)
% Output:
% S not-a-knot cubic spline (function)
t = t(:); y = y(:); % ensure column vectors
n = length(t) - 1;
h = diff(t); % differences of all adjacent pairs
% Preliminary definitions.
Z = zeros(n);
I = eye(n); E = I(1:n-1, :);
J = I - diag(ones(n-1, 1), 1);
H = diag(h);
% Left endpoint interpolation:
AL = [ I, Z, Z, Z ];
vL = y(1:n);
% Right endpoint interpolation:
AR = [ I, H, H^2, H^3 ];
vR = y(2:n+1);
% Continuity of first derivative:
A1 = E*[ Z, J, 2*H, 3*H^2 ];
v1 = zeros(n-1, 1);
% Continuity of second derivative:
A2 = E*[ Z, Z, J, 3*H ];
v2 = zeros(n-1, 1);
% Not-a-knot conditions:
nakL = [ zeros(1, 3*n), [1, -1, zeros(1, n-2)] ];
nakR = [ zeros(1, 3*n), [zeros(1, n-2), 1, -1] ];
% Assemble and solve the full system.
A = [ AL; AR; A1; A2; nakL; nakR ];
v = [ vL; vR; v1; v2; 0; 0 ];
z = A \ v;
% Break the coefficients into separate vectors.
rows = 1:n;
a = z(rows );
b = z(n + rows);
c = z(2 * n + rows);
d = z(3 * n + rows);
S = @evaulate;
% This function evaluates the spline when called with a value for x.
function f = evaulate(x)
f = zeros(size(x));
for k = 1:n % iterate over the pieces
% Evalaute this piece's cubic at the points inside it.
index = (x >= t(k)) & (x <= t(k+1));
f(index) = polyval( [d(k), c(k), b(k), a(k)], x(index) - t(k) );
end
end
end1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def spinterp(t, y):
"""
spinterp(t, y)
Create a cubic not-a-knot spline interpolating function for data values in y given at nodes in t.
"""
n = len(t) - 1
h = [t[i + 1] - t[i] for i in range(n)]
# Preliminary definitions.
Z = np.zeros([n, n])
I = np.eye(n)
E = I[: n - 1, :]
J = np.eye(n) + np.diag(-np.ones(n - 1), 1)
H = np.diag(h)
# Left endpoint interpolation:
AL = np.hstack([I, Z, Z, Z])
vL = y[:-1]
# Right endpoint interpolation:
AR = np.hstack([I, H, H**2, H**3])
vR = y[1:]
# Continuity of first derivative:
A1 = E @ np.hstack([Z, J, 2 * H, 3 * H**2])
v1 = np.zeros(n - 1)
# Continuity of second derivative:
A2 = E @ np.hstack([Z, Z, J, 3 * H])
v2 = np.zeros(n - 1)
# Not-a-knot conditions:
nakL = np.hstack([np.zeros(3 * n), np.hstack([1, -1, np.zeros(n - 2)])])
nakR = np.hstack([np.zeros(3 * n), np.hstack([np.zeros(n - 2), 1, -1])])
# Assemble and solve the full system.
A = np.vstack([AL, AR, A1, A2, nakL, nakR])
v = np.hstack([vL, vR, v1, v2, 0, 0])
z = solve(A, v)
# Break the coefficients into separate vectors.
rows = np.arange(n)
a = z[rows]
b = z[n + rows]
c = z[2 * n + rows]
d = z[3 * n + rows]
S = [np.poly1d([d[k], c[k], b[k], a[k]]) for k in range(n)]
# This function evaluates the spline when called with a value for x.
def evaluate(x):
f = np.zeros(x.shape)
for k in range(n):
# Evaluate this piece's cubic at the points inside it.
index = (x >= t[k]) & (x <= t[k + 1])
f[index] = S[k](x[index] - t[k])
return f
return evaluateFunction 5.3.1 E \mathbf{E} E S S S
5.3.4 Conditioning and convergence ¶ Example 5.3.1 (Cubic splines)
Example 5.3.1 For illustration, here is a spline interpolant using just a few nodes.
using Plots
f = x -> exp(sin(7x))
plot(f, 0, 1, label="function", xlabel=L"x", ylabel=L"y")
t = [0, 0.075, 0.25, 0.55, 0.7, 1] # nodes
y = f.(t) # values at nodes
scatter!(t, y, label="values at nodes")S = FNC.spinterp(t, y)
plot!(S, 0, 1, label="spline")Now we look at the convergence rate as the number of nodes increases.
x = (0:10000) / 1e4 # sample the difference at many points
n = @. round(Int, 2^(3:0.5:7)) # numbers of nodes
err = zeros(length(n))
for (k, n) in enumerate(n)
t = (0:n) / n
S = FNC.spinterp(t, f.(t))
dif = @. f(x) - S(x)
err[k] = norm(dif, Inf)
end
@pt :header=["n", "max-norm error"] [n[1:2:end] err[1:2:end]]Since we expect convergence that is O ( h 4 ) = O ( n − 4 ) O(h^4)=O(n^{-4}) O ( h 4 ) = O ( n − 4 )
order4 = @. (n / n[1])^(-4)
plot(n, [err order4];
m=[:o :none], l=[:solid :dash],
label=["error" "4th order"],
xaxis=(:log10, "n"), yaxis=(:log10, L"|| f-S\,||_\infty"),
title="Convergence of spline interpolation")Example 5.3.1 For illustration, here is a spline interpolant using just a few nodes.
clf
f = @(x) exp(sin(7 * x));
fplot(f, [0, 1], displayname="function")
t = [0, 0.075, 0.25, 0.55, 0.7, 1]; % nodes
y = f(t); % values at nodes
hold on, scatter(t, y, displayname="values at nodes")
S = spinterp(t, y);
fplot(S, [0, 1], displayname="spline")
xlabel("x"); ylabel("y")
legend();Unrecognized function or variable 'spinterp'. Now we look at the convergence rate as the number of nodes increases.
x = (0:10000)' / 1e4; % sample the difference at many points
n = round(2 .^ (3:0.5:7))'; % numbers of nodes
maxerr = zeros(size(n));
for i = 1:length(n)
t = (0:n(i))' / n(i);
S = spinterp(t, f(t));
err = f(x) - S(x);
maxerr(i) = norm(err, Inf);
end
table(n(1:2:end), maxerr(1:2:end), variableNames=["n", "inf-norm error"])Unrecognized function or variable 'spinterp'. Since we expect convergence that is O ( h 4 ) = O ( n − 4 ) O(h^4)=O(n^{-4}) O ( h 4 ) = O ( n − 4 )
clf
loglog(n, maxerr, "-o", displayname="error")
order4 = 0.5 * maxerr(end) * (n / n(end)) .^ (-4);
hold on
loglog(n, order4, "k--", displayname="O(n^{-4})")
xlabel("n"); ylabel("|| f-S ||_{\infty}")
title(("Convergence of spline interpolation"));Example 5.3.1 For illustration, here is a spline interpolant using just a few nodes.
f = lambda x: exp(sin(7 * x))
x = linspace(0, 1, 500)
fig, ax = subplots()
ax.plot(x, f(x), label="function")
t = array([0, 0.075, 0.25, 0.55, 0.7, 1]) # nodes
y = f(t) # values at nodes
xlabel("$x$")
ylabel("$y$")
ax.scatter(t, y, label="nodes")S = FNC.spinterp(t, y)
ax.plot(x, S(x), label="spline")
ax.legend()
figNow we look at the convergence rate as the number of nodes increases.
N = floor(2 ** linspace(3, 8, 17)).astype(int)
err = zeros(N.size)
for i, n in enumerate(N):
t = linspace(0, 1, n + 1) # interpolation nodes
p = FNC.spinterp(t, f(t))
err[i] = max(abs(f(x) - p(x)))
print(err)[3.05633432e-02 2.39601586e-02 1.68054365e-02 7.64098319e-03
2.89472870e-03 1.34574135e-03 5.43142890e-04 2.28104055e-04
9.17629364e-05 3.71552636e-05 1.56015311e-05 6.34890672e-06
2.53866817e-06 9.98323636e-07 4.35498457e-07 1.75251504e-07
6.59321329e-08]
Since we expect convergence that is O ( h 4 ) = O ( n − 4 ) O(h^4)=O(n^{-4}) O ( h 4 ) = O ( n − 4 )
order4 = (N / N[0]) ** (-4)
loglog(N, err, "-o", label="observed error")
loglog(N, order4, "--", label="4th order")
xlabel("$n$")
ylabel("$\|f-S\|_\infty$")
legend();<>:5: SyntaxWarning: invalid escape sequence '\|'
<>:5: SyntaxWarning: invalid escape sequence '\|'
/var/folders/0m/fy_f4_rx5xv68sdkrpcm26r00000gq/T/ipykernel_23602/3396180238.py:5: SyntaxWarning: invalid escape sequence '\|'
ylabel("$\|f-S\|_\infty$")
Besides having more smoothness than a piecewise linear interpolant, the not-a-knot cubic spline improves the order of accuracy to 4.
Suppose that f ( x ) f(x) f ( x ) [ a , b ] [a,b] [ a , b ] f ∈ C 4 [ a , b ] f\in C^4[a,b] f ∈ C 4 [ a , b ] S n ( x ) S_n(x) S n ( x ) ( t i , f ( t i ) ) \bigl(t_i,f(t_i)\bigr) ( t i , f ( t i ) ) i = 0 , … , n i=0,\ldots,n i = 0 , … , n t i = a + i h t_i=a+i h t i = a + ih h = ( b − a ) / n h=(b-a)/n h = ( b − a ) / n h h h C > 0 C>0 C > 0
∥ f − S n ∥ ∞ ≤ C h 4 . \bigl\| f - S_n \bigr\|_\infty \le Ch^4. ∥ ∥ f − S n ∥ ∥ ∞ ≤ C h 4 . The conditioning of spline interpolation is much more complicated than for the piecewise linear case. First, the fact that the coefficients of all the cubics must be solved for simultaneously implies that each data value in y \mathbf{y} y S S S S S S y \mathbf{y} y Exercise 5.3.5
5.3.5 Exercises ¶ ✍ In each case, write out the entries of the matrix and right-hand side of the linear system that determines the coefficients for the cubic not-a-knot spline interpolant of the given function and node vector.
(a) cos ( π 2 x 2 ) , t = [ − 1 , 1 , 4 ] \cos (\pi^2 x^2 ), \: \mathbf{t} = [-1,1,4] cos ( π 2 x 2 ) , t = [ − 1 , 1 , 4 ]
(b) cos ( π 2 x 2 ) , t = [ 0 , 1 2 , 3 4 , 1 ] \cos (\pi^2 x^2), \: \mathbf{t} = [0,\tfrac{1}{2},\tfrac{3}{4},1] cos ( π 2 x 2 ) , t = [ 0 , 2 1 , 4 3 , 1 ]
(c) ln ( x ) , t = [ 1 , 2 , 3 ] \ln(x), \: \mathbf{t} = [1,2,3] ln ( x ) , t = [ 1 , 2 , 3 ]
(d) sin ( x 2 ) , t = [ − 1 , 0 , 1 ] \sin(x^2),\: \mathbf{t} = [-1,0,1] sin ( x 2 ) , t = [ − 1 , 0 , 1 ]
⌨ For each given function, interval, and value of n n n n + 1 n+1 n + 1 Function 5.3.1
(a) cos ( π x 2 ) \cos(\pi x^2) cos ( π x 2 ) [ 0 , 4 ] [0,4] [ 0 , 4 ] n = 18 n=18 n = 18
(b) ln ( x ) \ln(x) ln ( x ) [ 1 , 20 ] [1,20] [ 1 , 20 ] n = 4 n=4 n = 4
(c) sin ( 1 x ) \sin\left(\frac{1}{x}\right) sin ( x 1 ) [ 1 2 , 7 ] \left[\frac{1}{2},7\right] [ 2 1 , 7 ] n = 9 n=9 n = 9
⌨ For each given function and interval, perform piecewise linear interpolation using Function 5.3.1 n + 1 n+1 n + 1 n = 10 , 20 , 40 , 80 , 160 , 320 n=10,20,40,80,160,320 n = 10 , 20 , 40 , 80 , 160 , 320 n n n
E ( n ) = ∥ f − p ∥ ∞ = max x ∣ f ( x ) − p ( x ) ∣ E(n) = \| f-p \|_\infty = \max_x | f(x) - p(x) | E ( n ) = ∥ f − p ∥ ∞ = x max ∣ f ( x ) − p ( x ) ∣ by evaluating the function and interpolant at 1600 points in the interval. Make a log-log plot of E E E n n n E = C n − 4 E=Cn^{-4} E = C n − 4 C C C
(a) cos ( π x 2 ) \cos(\pi x^2) cos ( π x 2 ) [ 0 , 4 ] [0,4] [ 0 , 4 ]
(b) ln ( x ) \ln(x) ln ( x ) [ 1 , 20 ] [1,20] [ 1 , 20 ]
(c) sin ( 1 x ) \sin\left(\frac{1}{x}\right) sin ( x 1 ) [ 1 2 , 7 ] \left[\frac{1}{2},7\right] [ 2 1 , 7 ]
⌨ Although the cardinal cubic splines are intractable in closed form, they can be found numerically. Each cardinal spline interpolates the data from one column of an identity matrix. Define the nodes t = [ 0 , 0.075 , 0.25 , 0.55 , 1 ] \mathbf{t} = \bigl[0,\, 0.075,\, 0.25,\, 0.55,\, 1] t = [ 0 , 0.075 , 0.25 , 0.55 , 1 ] [ 0 , 1 ] [0,1] [ 0 , 1 ] [ 0 , 1 ] [0,1] [ 0 , 1 ]
(a) ✍ If y 0 = y n y_0=y_n y 0 = y n S ( x ) S(x) S ( x ) S ′ S' S ′ S ′ ′ S'' S ′′
(b) ⌨ Modify Function 5.3.1 f ( x ) = exp ( sin ( 3 x ) ) f(x) =\exp(\sin(3x)) f ( x ) = exp ( sin ( 3 x )) [ 0 , 2 π / 3 ] [0,2\pi/3] [ 0 , 2 π /3 ] n = 8 n=8 n = 8
This explains the name of the not-a-knot spline—for splines, “knots” are the points at which different piecewise definitions meet.