In Finite differences we used finite differences to turn a discrete collection of function values into an estimate of the derivative of the function at a point. Just as with differentiation in elementary calculus, we can generalize differences at a point into an operation that maps discretized functions to discretized derivatives.
10.3.1 Matrices for finite differences ¶ We first discretize the interval x ∈ [ a , b ] x\in[a,b] x ∈ [ a , b ] into equal pieces of length h = ( b − a ) / n h=(b-a)/n h = ( b − a ) / n , leading to the nodes
x i = a + i h , i = 0 , … , n . x_i=a+i h, \qquad i=0,\ldots,n. x i = a + ih , i = 0 , … , n . Note again that our node indexing scheme starts at zero. Our goal is to find a vector g \mathbf{g} g such that g i ≈ f ′ ( x i ) g_i \approx f'(x_i) g i ≈ f ′ ( x i ) for i = 0 , … , n i=0,\ldots,n i = 0 , … , n . Our first try is the forward difference formula (5.4.2) ,
g i = f i + 1 − f i h , t = 0 , … , n − 1. g_i = \frac{f_{i+1}-f_i}{h}, \qquad t=0,\ldots,n-1. g i = h f i + 1 − f i , t = 0 , … , n − 1. However, this leaves g n g_n g n undefined, because the formula would refer to the unavailable value f n + 1 f_{n+1} f n + 1 . For g n g_n g n we could resort to the backward difference
g n = f n − f n − 1 h . g_n = \frac{f_{n}-f_{n-1}}{h}. g n = h f n − f n − 1 . We can summarize the entire set of formulas by defining
f = [ f ( x 0 ) f ( x 1 ) ⋮ f ( x n − 1 ) f ( x n ) ] , \mathbf{f} =
\begin{bmatrix}
f(x_0) \\[1mm] f(x_1) \\[1mm] \vdots \\[1mm] f(x_{n-1}) \\[1mm] f(x_n)
\end{bmatrix}, \quad f = ⎣ ⎡ f ( x 0 ) f ( x 1 ) ⋮ f ( x n − 1 ) f ( x n ) ⎦ ⎤ , and then the vector equation
[ f ′ ( x 0 ) f ′ ( x 1 ) ⋮ f ′ ( x n − 1 ) f ′ ( x n ) ] ≈ D x f , D x = 1 h [ − 1 1 − 1 1 ⋱ ⋱ − 1 1 − 1 1 ] . \begin{bmatrix}
f'(x_0) \\[1mm] f'(x_1) \\[1mm] \vdots \\[1mm] f'(x_{n-1}) \\[1mm] f'(x_n)
\end{bmatrix}
\approx
\mathbf{D}_x \mathbf{f}, \qquad
\mathbf{D}_x
= \frac{1}{h}
\begin{bmatrix}
-1 & 1 & & & \\[1mm]
& -1 & 1 & & \\[1mm]
& & \ddots & \ddots & \\[1mm]
& & & -1 & 1 \\[1mm]
& & & -1 & 1
\end{bmatrix}. ⎣ ⎡ f ′ ( x 0 ) f ′ ( x 1 ) ⋮ f ′ ( x n − 1 ) f ′ ( x n ) ⎦ ⎤ ≈ D x f , D x = h 1 ⎣ ⎡ − 1 1 − 1 1 ⋱ ⋱ − 1 − 1 1 1 ⎦ ⎤ . Here as elsewhere, elements of D x \mathbf{D}_x D x that are not shown are zero. We call D x \mathbf{D}_x D x a differentiation matrix . Each row of D x \mathbf{D}_x D x gives the weights of the finite-difference formula being used at one of the nodes.
The differentiation matrix D x \mathbf{D}_x D x in (10.3.5) is not a unique choice. We are free to use whatever finite-difference formulas we like in each row. However, it makes sense to choose rows that are as similar as possible. Using second-order centered differences where possible and second-order one-sided formulas (see Table 5.4.2 ) at the boundary points leads to
D x = 1 h [ − 3 2 2 − 1 2 − 1 2 0 1 2 − 1 2 0 1 2 ⋱ ⋱ ⋱ − 1 2 0 1 2 1 2 − 2 3 2 ] . \mathbf{D}_x = \frac{1}{h}
\begin{bmatrix}
-\frac{3}{2} & 2 & -\frac{1}{2} & & & \\[1mm]
-\frac{1}{2} & 0 & \frac{1}{2} & & & \\[1mm]
& -\frac{1}{2} & 0 & \frac{1}{2} & & \\
& & \ddots & \ddots & \ddots & \\
& & & -\frac{1}{2} & 0 & \frac{1}{2} \\[1mm]
& & & \frac{1}{2} & -2 & \frac{3}{2}
\end{bmatrix}. D x = h 1 ⎣ ⎡ − 2 3 − 2 1 2 0 − 2 1 − 2 1 2 1 0 ⋱ 2 1 ⋱ − 2 1 2 1 ⋱ 0 − 2 2 1 2 3 ⎦ ⎤ . The differentiation matrices so far are banded matrices, i.e., all the nonzero values are along diagonals close to the main diagonal.[1]
10.3.2 Second derivative ¶ Similarly, we can define differentiation matrices for second derivatives. For example,
[ f ′ ′ ( x 0 ) f ′ ′ ( x 1 ) f ′ ′ ( x 2 ) ⋮ f ′ ′ ( x n − 1 ) f ′ ′ ( x n ) ] ≈ 1 h 2 [ 2 − 5 4 − 1 1 − 2 1 1 − 2 1 ⋱ ⋱ ⋱ 1 − 2 1 − 1 4 − 5 2 ] [ f ( x 0 ) f ( x 1 ) f ( x 2 ) ⋮ f ( x n − 1 ) f ( x n ) ] = D x x f . \begin{bmatrix}
f''(x_0) \\[1mm] f''(x_1) \\[1mm] f''(x_2) \\[1mm] \vdots \\[1mm] f''(x_{n-1}) \\[1mm] f''(x_n)
\end{bmatrix}
\approx
\frac{1}{h^2}
\begin{bmatrix}
2 & -5 & 4 & -1 & & \\[1mm]
1 & -2 & 1 & & & \\[1mm]
& 1 & -2 & 1 & & \\[1mm]
& & \ddots & \ddots & \ddots & \\[1mm]
& & & 1 & -2 & 1 \\[1mm]
& & -1 & 4 & -5 & 2
\end{bmatrix}
\begin{bmatrix}
f(x_0) \\[1mm] f(x_1) \\[1mm] f(x_2) \\[1mm] \vdots \\[1mm] f(x_{n-1}) \\[1mm] f(x_n)
\end{bmatrix} = \mathbf{D}_{xx} \mathbf{f}. ⎣ ⎡ f ′′ ( x 0 ) f ′′ ( x 1 ) f ′′ ( x 2 ) ⋮ f ′′ ( x n − 1 ) f ′′ ( x n ) ⎦ ⎤ ≈ h 2 1 ⎣ ⎡ 2 1 − 5 − 2 1 4 1 − 2 ⋱ − 1 − 1 1 ⋱ 1 4 ⋱ − 2 − 5 1 2 ⎦ ⎤ ⎣ ⎡ f ( x 0 ) f ( x 1 ) f ( x 2 ) ⋮ f ( x n − 1 ) f ( x n ) ⎦ ⎤ = D xx f . We have multiple choices again for D x x \mathbf{D}_{xx} D xx , and it need not be the square of any particular D x \mathbf{D}_x D x . As pointed out in Finite differences , squaring the first derivative is a valid approach but would place entries in D x x \mathbf{D}_{xx} D xx farther from the diagonal than is necessary.
Together the matrices (10.3.6) and (10.3.7) give second-order approximations of the first and second derivatives at all nodes. These matrices, as well as the nodes x 0 , … , x n x_0,\ldots,x_n x 0 , … , x n , are returned by Function 10.3.1 .
Example 10.3.1 (Differentiation matrices)
Example 10.3.1 We test first-order and second-order differentiation matrices for the function x + exp ( sin 4 x ) x + \exp(\sin 4x) x + exp ( sin 4 x ) over [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] .
f = x -> x + exp(sin(4 * x));
For reference, here are the exact first and second derivatives.
df_dx = x -> 1 + 4 * exp(sin(4x)) * cos(4x);
d2f_dx2 = x -> 4 * exp(sin(4x)) * (4 * cos(4x)^2 - 4 * sin(4x));
We discretize on equally spaced nodes and evaluate f f f at the nodes.
t, Dₓ, Dₓₓ = FNC.diffmat2(18, [-1, 1])
y = f.(t);
Then the first two derivatives of f f f each require one matrix-vector multiplication.
yₓ = Dₓ * y
yₓₓ = Dₓₓ * y;
The results show poor accuracy for this small value of n n n .
using Plots
plot(df_dx, -1, 1, layout = 2, xaxis = (L"x"), yaxis = (L"f'(x)"))
scatter!(t, yₓ, subplot = 1)
plot!(d2f_dx2, -1, 1, subplot = 2, xaxis = (L"x"), yaxis = (L"f''(x)"))
scatter!(t, yₓₓ, subplot = 2)
A convergence experiment confirms the order of accuracy. Because we expect an algebraic convergence rate, we use a log-log plot of the errors.
n = @. round(Int, 2^(4:0.5:11))
err = zeros(length(n), 2)
for (k, n) in enumerate(n)
t, Dₓ, Dₓₓ = FNC.diffmat2(n, [-1, 1])
y = f.(t)
err[k, 1] = norm(df_dx.(t) - Dₓ * y, Inf)
err[k, 2] = norm(d2f_dx2.(t) - Dₓₓ * y, Inf)
end
plot(n, err, m = :o, label = [L"f'" L"f''"])
plot!(n, 10 * 10 * n .^ (-2);
l = (:dash, :black),
label = "2nd order",
xaxis = (:log10, "n"),
yaxis = (:log10, "max error"),
title = "Convergence of finite differences")
Example 10.3.1 We test first-order and second-order differentiation matrices for the function x + exp ( sin 4 x ) x + \exp(\sin 4x) x + exp ( sin 4 x ) over [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] .
f = @(x) x + exp(sin(4*x));
For reference, here are the exact first and second derivatives.
df_dx = @(x) 1 + 4 * exp(sin(4*x)) .* cos(4*x);
d2f_dx2 = @(x) 4 * exp(sin(4*x)) .* (4*cos(4*x).^2 - 4*sin(4*x));
We discretize on equally spaced nodes and evaluate f f f at the nodes.
[t, Dx, Dxx] = diffmat2(12, [-1 1]);
y = f(t);
Not enough input arguments.
Error in LiveEditorEvaluationHelperEc482176a-104e-4b4a-bc53-b10e8b2e2b15>@(t,u,p)Dxx*u (line 2)
f = @(t, u, p) Dxx * u;
Then the first two derivatives of f f f each require one matrix-vector multiplication.
yx = Dx * y;
yxx = Dxx * y;
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.
The results show poor accuracy for this small value of n n n .
clf, subplot(2, 1, 1)
fplot(df_dx, [-1, 1]), hold on
plot(t, yx, 'ko')
xlabel('x'), ylabel('f''(x)')
subplot(2, 1, 2)
fplot(d2f_dx2, [-1, 1]), hold on
plot(t, yxx, 'ko')
xlabel('x'), ylabel('f''''(x)')
Unrecognized function or variable 'yx'.
An convergence experiment confirms the order of accuracy. Because we expect an algebraic convergence rate, we use a log-log plot of the errors.
n = round( 2 .^ (4:.5:11)' );
err = zeros(length(n), 2);
for k = 1:length(n)
[t, Dx, Dxx] = diffmat2(n(k), [-1, 1]);
y = f(t);
err(k, 1) = norm(df_dx(t) - Dx * y, Inf);
err(k, 2) = norm(d2f_dx2(t) - Dxx * y, Inf);
end
clf
loglog(n, err, 'o-'), hold on
loglog(n, 100 * n.^(-2), 'k--')
legend("f'", "f''", '2nd order')
xlabel('n'), ylabel('max error')
title('Convergence of finite differences')
Example 10.3.1 We test first-order and second-order differentiation matrices for the function x + exp ( sin 4 x ) x + \exp(\sin 4x) x + exp ( sin 4 x ) over [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] .
f = lambda x: x + exp( sin(4 * x) )
For reference, here are the exact first and second derivatives.
df_dx = lambda x: 1 + 4 * exp(sin(4 * x)) * cos(4 * x)
d2f_dx2 = lambda x: 4 * exp(sin(4 * x)) * (4 * cos(4 * x)**2 - 4 * sin(4 * x))
We discretize on equally spaced nodes and evaluate f f f at the nodes.
t, Dx, Dxx = FNC.diffmat2(12, [-1, 1])
y = f(t)
Then the first two derivatives of f f f each require one matrix-vector multiplication.
yx = Dx @ y
yxx = Dxx @ y
The results show poor accuracy for this small value of n n n .
x = linspace(-1, 1, 500)
subplot(2, 1, 1)
plot(x, df_dx(x))
plot(t, yx, "ko")
xlabel("$x$"), ylabel("$f'(x)$")
subplot(2, 1, 2)
plot(x, d2f_dx2(x))
plot(t, yxx, "ko")
xlabel("$x$"), ylabel("$f''(x)$");
A convergence experiment confirms the order of accuracy. Because we expect an algebraic convergence rate, we use a log-log plot of the errors.
N = array([int(2**k) for k in arange(4, 11.5, 0.5)])
err1 = zeros(len(N))
err2 = zeros(len(N))
for k, n in enumerate(N):
t, Dx, Dxx = FNC.diffmat2(n, [-1, 1])
y = f(t)
err1[k] = norm(df_dx(t) - Dx @ y, inf)
err2[k] = norm(d2f_dx2(t) - Dxx @ y, inf)
loglog(N, err1, "-o", label="$f'$")
loglog(N, err2, "-o", label="$f''$")
plot(N, 10 * 10 / N**2, "k--", label="2nd order")
xlabel("$n$"), ylabel("max error")
legend(loc="lower left")
title("Convergence of finite differences");
10.3.3 Spectral differentiation ¶ Recall that finite-difference formulas are derived in three steps:
Choose a node index set S S S near node i i i . Interpolate with a polynomial using the nodes in S S S . Differentiate the interpolant and evaluate at node i i i . We can modify this process by using a global interpolant, either polynomial or trigonometric, as in Chapter 9 . Rather than choosing a different index set for each node, we use all of the nodes each time.
In a nonperiodic setting we use Chebyshev second-kind points for stability:
x k = − cos ( k π n ) , k = 0 , … , n ; x_k = -\cos\left(\frac{k \pi}{n}\right), \qquad k=0,\ldots,n; x k = − cos ( n kπ ) , k = 0 , … , n ; then the resulting Chebyshev differentiation matrix has entries
D 00 = 2 n 2 + 1 6 , D n n = − 2 n 2 + 1 6 , D i j = { − x i 2 ( 1 − x i 2 ) , i = j , c i c j ( − 1 ) i + j x i − x j , i ≠ j , \begin{gathered}
D_{00} = \dfrac{2n^2+1}{6}, \qquad D_{n n} = -\dfrac{2n^2+1}{6}, \\
D_{ij} =
\begin{cases}
-\dfrac{x_i}{2(1-x_i^2)}, & i=j, \\[4mm]
\dfrac{c_i}{c_j}\, \dfrac{(-1)^{i+j}}{x_i-x_j}, & i\neq j,
\end{cases}
\end{gathered} D 00 = 6 2 n 2 + 1 , D nn = − 6 2 n 2 + 1 , D ij = ⎩ ⎨ ⎧ − 2 ( 1 − x i 2 ) x i , c j c i x i − x j ( − 1 ) i + j , i = j , i = j , where c 0 = c n = 2 c_0=c_n=2 c 0 = c n = 2 and c i = 1 c_i=1 c i = 1 for i = 1 , … , n − 1 i=1,\ldots,n-1 i = 1 , … , n − 1 . Note that this matrix is dense. The simplest way to compute a second derivative is by squaring D x \mathbf{D}_x D x , as there is no longer any concern about the bandwidth of the result.
Function 10.3.2 returns these two matrices. The function uses a change of variable to transplant the standard [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] for Chebyshev nodes to any [ a , b ] [a,b] [ a , b ] . It also takes a different approach to computing the diagonal elements of D x \mathbf{D}_x D x than the formulas in (10.3.9) (see Exercise 5 ).
Chebyshev differentiation matrices1
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
"""
diffcheb(n, xspan)
Compute Chebyshev differentiation matrices on `n`+1 points in the
interval `xspan`. Returns a vector of nodes and the matrices for the
first and second derivatives.
"""
function diffcheb(n, xspan)
x = [-cos(k * π / n) for k in 0:n] # nodes in [-1,1]
# Off-diagonal entries.
c = [2; ones(n - 1); 2] # endpoint factors
dij = (i, j) -> (-1)^(i + j) * c[i+1] / (c[j+1] * (x[i+1] - x[j+1]))
Dₓ = [dij(i, j) for i in 0:n, j in 0:n]
# Diagonal entries.
Dₓ[isinf.(Dₓ)] .= 0 # fix divisions by zero on diagonal
s = sum(Dₓ, dims = 2)
Dₓ -= diagm(s[:, 1]) # "negative sum trick"
# Transplant to [a,b].
a, b = xspan
x = @. a + (b - a) * (x + 1) / 2
Dₓ = 2 * Dₓ / (b - a) # chain rule
# Second derivative.
Dₓₓ = Dₓ^2
return x, Dₓ, Dₓₓ
end
Chebyshev differentiation matrices1
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
function [x,Dx,Dxx] = diffcheb(n,xspan)
%DIFFCHEB Chebyshev differentiation matrices.
% Input:
% n number of subintervals (integer)
% xspan interval endpoints (vector)
% Output:
% x Chebyshev nodes in domain (length n+1)
% Dx matrix for first derivative (n+1 by n+1)
% Dxx matrix for second derivative (n+1 by n+1)
x = -cos( (0:n)'*pi/n ); % nodes in [-1,1]
Dx = zeros(n+1);
c = [2; ones(n-1,1); 2]; % endpoint factors
i = (0:n)'; % row indices
% Off-diagonal entries
for j = 0:n
num = c(i+1).*(-1).^(i+j);
den = c(j+1)*(x-x(j+1));
Dx(:,j+1) = num./den;
end
% Diagonal entries
Dx(isinf(Dx)) = 0; % fix divisions by zero on diagonal
Dx = Dx - diag(sum(Dx,2)); % "negative sum trick"
% Transplant to [a,b]
a = xspan(1); b = xspan(2);
x = a + (b-a)*(x+1)/2;
Dx = 2*Dx/(b-a);
% Second derivative
Dxx = Dx^2;
Chebyshev differentiation matrices1
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
def diffcheb(n, xspan):
"""
diffcheb(n, xspan)
Compute Chebyshev differentiation matrices on n+1 points in the
interval xspan. Return a vector of nodes, and the matrices for the first
and second derivatives.
"""
x = -np.cos(np.arange(n + 1) * np.pi / n) # nodes in [-1,1]
Dx = np.zeros([n + 1, n + 1])
c = np.hstack([2.0, np.ones(n - 1), 2.0]) # endpoint factors
# Off-diagonal entries
Dx = np.zeros([n + 1, n + 1])
for i in range(n + 1):
for j in range(n + 1):
if i != j:
Dx[i, j] = (-1) ** (i + j) * c[i] / (c[j] * (x[i] - x[j]))
# Diagonal entries by the "negative sum trick"
for i in range(n + 1):
Dx[i, i] = -np.sum([Dx[i, j] for j in range(n + 1) if j != i])
# Transplant to [a,b]
a, b = xspan
x = a + (b - a) * (x + 1) / 2
Dx = 2 * Dx / (b - a)
# Second derivative
Dxx = Dx @ Dx
return x, Dx, Dxx
Example 10.3.2 (Chebyshev differentiation matrices)
Example 10.3.2 Here is a 4 × 4 4\times 4 4 × 4 Chebyshev differentiation matrix.
t, Dₓ = FNC.diffcheb(3, [-1, 1])
Dₓ
4×4 Matrix{Float64}:
-3.16667 4.0 -1.33333 0.5
-1.0 0.333333 1.0 -0.333333
0.333333 -1.0 -0.333333 1.0
-0.5 1.33333 -4.0 3.16667
We again test the convergence rate.
f = x -> x + exp(sin(4 * x));
df_dx = x -> 1 + 4 * exp(sin(4 * x)) * cos(4 * x);
d2f_dx2 = x -> 4 * exp(sin(4 * x)) * (4 * cos(4 * x)^2 - 4 * sin(4 * x));
n = 5:5:70
err1 = zeros(size(n))
err2 = zeros(size(n))
for (k, n) in enumerate(n)
t, Dₓ, Dₓₓ = FNC.diffcheb(n, [-1, 1])
y = f.(t)
err1[k] = norm(df_dx.(t) - Dₓ * y, Inf)
err2[k] = norm(d2f_dx2.(t) - Dₓₓ * y, Inf)
end
Since we expect a spectral convergence rate, we use a semi-log plot for the error.
plot(n, [err1 err2]; m = :o,
label=[L"f'" L"f''"],
xaxis=(L"n"), yaxis = (:log10, "max error"),
title="Convergence of Chebyshev derivatives")
Example 10.3.2 Here is a 4 × 4 4\times 4 4 × 4 Chebyshev differentiation matrix.
[t, Dx] = diffcheb(3, [-1, 1]);
format rat
Dx
We again test the convergence rate.
f = @(x) x + exp(sin(4*x));
df_dx = @(x) 1 + 4 * exp(sin(4*x)) .* cos(4*x);
d2f_dx2 = @(x) 4 * exp(sin(4*x)) .* (4*cos(4*x).^2 - 4*sin(4*x));
n = 5:5:70;
err = zeros(length(n), 2);
for k = 1:length(n)
[t, Dx, Dxx] = diffcheb(n(k), [-1, 1]);
y = f(t);
err(k, 1) = norm(df_dx(t) - Dx * y, Inf);
err(k, 2) = norm(d2f_dx2(t) - Dxx * y, Inf);
end
Since we expect a spectral convergence rate, we use a semi-log plot for the error.
clf, format
semilogy(n, err, 'o-'), hold on
legend("f'", "f''")
xlabel('n'), ylabel('max error')
title('Convergence of finite differences')
Example 10.3.2 Here is a 4 × 4 4\times 4 4 × 4 Chebyshev differentiation matrix.
t, Dx, Dxx = FNC.diffcheb(3, [-1, 1])
print(Dx)
[[-3.16666667 4. -1.33333333 0.5 ]
[-1. 0.33333333 1. -0.33333333]
[ 0.33333333 -1. -0.33333333 1. ]
[-0.5 1.33333333 -4. 3.16666667]]
We again test the convergence rate.
f = lambda x: x + exp(sin(4 * x))
df_dx = lambda x: 1 + 4 * exp(sin(4 * x)) * cos(4 * x)
d2f_dx2 = lambda x: 4 * exp(sin(4 * x)) * (4 * cos(4 * x) ** 2 - 4 * sin(4 * x))
N = range(5, 75, 5)
err1 = zeros(len(N))
err2 = zeros(len(N))
err = zeros((len(N), 2))
for k, n in enumerate(N):
t, Dx, Dxx = FNC.diffcheb(n, [-1, 1])
y = f(t)
err[k, 0] = norm(df_dx(t) - Dx @ y, inf)
err[k, 1] = norm(d2f_dx2(t) - Dxx @ y, inf)
Since we expect a spectral convergence rate, we use a semi-log plot for the error.
semilogy(N, err, "-o")
xlabel("$n$"), ylabel("max error")
legend(["$f'$", "$f''$"], loc="lower left")
title("Convergence of Chebyshev derivatives");
According to Theorem 9.3.1 , the convergence of polynomial interpolation to f f f using Chebyshev nodes is spectral if f f f is analytic (at least having infinitely many derivatives) on the interval. The derivatives of f f f are also approximated with spectral accuracy.
10.3.4 Exercises ¶ (a) ✍ Calculate D x 2 \mathbf{D}_x^2 D x 2 using (10.3.5) to define D x \mathbf{D}_x D x .
(b) ⌨ Repeat the experiment of Demo 10.3.1 , but using this version of D x 2 \mathbf{D}_x^2 D x 2 to estimate f ′ ′ f'' f ′′ . What is the apparent order of accuracy in f ′ ′ f'' f ′′ ?
(a) ✍ Find the derivative of f ( x ) = sign ( x ) x 2 f(x) =\operatorname{sign}(x)x^2 f ( x ) = sign ( x ) x 2 on the interval [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] . (If this gives you trouble, use an equivalent piecewise definition of f f f .) What is special about this function at x = 0 x=0 x = 0 ?
(b) ⌨ Adapt Demo 10.3.1 to operate on the function from part (a), computing only the first derivative. What is the observed order of accuracy?
(c) ✍ Show that for even values of n n n , there is only one node at which the error for computing f ′ f' f ′ in part (b) is nonzero.
⌨ To get a fourth-order accurate version of D x \mathbf{D}_x D x , five points per row are needed, including two special rows at each boundary. For a fourth-order D x x \mathbf{D}_{xx} D xx , five symmetric points per row are needed for interior rows and six points are needed for the rows near a boundary.
(a) Modify Function 10.3.1 to a function diffmat4
, which outputs fourth-order accurate differentiation matrices. You may want to use Function 5.4.1 .
(b) Repeat the experiment of Demo 10.3.1 using diffmat4
in place of Function 10.3.1 , and compare observed errors to fourth-order accuracy.
✍ Explain in detail how lines 23-24 in Function 10.3.2 correctly change the interval from [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] to [ a , b ] [a,b] [ a , b ] .
(a) ✍ What is the derivative of a constant function?
(b) ✍ Explain why for any reasonable differentiation matrix D \mathbf{D} D , we should find ∑ j = 0 n D i j = 0 \displaystyle \sum_{j=0}^nD_{ij}=0 j = 0 ∑ n D ij = 0 for all i i i .
(c) ✍ What does this have to do with Function 10.3.2 ? Refer to specific line(s) in the function for your answer.
Define the ( n + 1 ) × ( n + 1 ) (n+1)\times (n+1) ( n + 1 ) × ( n + 1 ) matrix T = [ 1 1 1 ⋮ ⋱ 1 1 ⋯ 1 ] \mathbf{T} = \displaystyle \begin{bmatrix}
1 & & & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & 1 & \cdots & 1
\end{bmatrix} T = ⎣ ⎡ 1 1 ⋮ 1 1 1 ⋱ ⋯ 1 ⎦ ⎤ .
(a) ✍ Write out T u \mathbf{T}\mathbf{u} Tu for a generic vector u \mathbf{u} u (start with a zero index). How is this like integration?
(b) ✍ Find the inverse of T \mathbf{T} T for any n n n . (You can use Julia to find the pattern, but show that the result is correct in general.) What does this have to do with the inverse of integration?