Before you turn this problem in, make sure everything runs as expected. First, restart the kernel (in the menubar, select Kernel → Restart) and then run all cells (in the menubar, select Cell → Run All).
Make sure you fill in any place that says BEGIN SOLUTION and END SOLUTION, as well as your name and collaborators below:
[ ]:
NAME = ""
COLLABORATORS = ""
2025-10-31 Piecewise Interpolation¶
We covered global polynomial interpolation, piecwise constant interpolation, and piecewise spline interpolation. In this assignment, we will explore piecewise polynomial interpolation with each piece being a polynomial of degree three or higher. We will then use these piecewise polynomials for interpolation.
To start, we will use our matrix for polynomial interpolation with Chebyshev polynomials from class. The Chebyshev polynomials can be generated with the recursion
and our interpolating polynomial is given by \(p \left( x \right) = \sum_{i = 0}^{k-1} c_i T_i \left( x \right)\), where \(k\) is the number of interpolation points and \(k - 1\) is the highest order (degree) of our Chebyshev polynomials.
Recall that with this Vandermonde matrix for the Chebyshev polynomials, we create a matrix equation \(A c = p\), where \(A\) is the Vandermonde matrix with columns given by the Chebyshev polynomials \(T_i \left( x \right)\) and the rows given by our interpolation points \(x_j\), \(c\) is the vector of coefficients for each individual Chebyshev polynomial, and \(p\) is the values of our interpolating polynomial \(p \left( x \right)\) at our interpolation points \(x_j\).
Fitting on one interval¶
[ ]:
using LinearAlgebra
using Plots
using Test
# Our Vandermonde matrix for Chebyshev polynomials
function vander_chebyshev(x, k=nothing)
"""Return a Vandermonde matrix of the first n Chebyshev polynomials evaluated at the points x.
The returned matrix will be of size n x k, where n is the length of the vector x.
If no k is provided, the returned matrix will be square.
"""
if isnothing(k)
k = length(x)
end
n = length(x)
A = ones(n, k)
A[:, 2] = x
for i in 3:k
A[:, i] = 2 * x .* A[:, i-1] - A[:, i-2]
end
A
end
# Let's try it
k = 4
A = vander_chebyshev(Vector(LinRange(-1, 1, k)))
println("A =")
display(A);
We’ll solve a small interpolation problem to make sure this function works as expected.
[ ]:
# The target function
k = 5
x = LinRange(-1, 1, k)
y = sin.(x)
# The Vandermonde matrix
A = vander_chebyshev(x)
# Solve for the coefficents
c = inv(A) * y
# And use the coefficients
scatter(x, y, label="data")
x = LinRange(-1, 1, 100)
p = vander_chebyshev(x, k) * c
plot!(x, p, label="fit")
Fitting on multiple subintervals¶
We can now use this Vandermonde matrix to fix polynomials piecewise.
Consider the interval from \(\left[ -1, 3 \right]\). We can divide this into two sub-intervals, \(\left[ -1, 1 \right]\) and \(\left[ 1, 3 \right]\) and fit a polynomial of degree \(k\) to each subinterval via our Vandermonde matrix. We will have separate sets of coefficients \(c_{0, i}\) and \(c_{1, i}\) such that
and
For simplicity, we’ll use a small ‘trick’. Instead of creating a new Vandermonde matrix for each subinterval, we’ll just reuse the same Vandermonde matrix on the interval \(\left[ -1, 1 \right]\) and fit two separate sets of coefficients. When evaluating these polynomials, we’ll shift the input \(x\) back into the range \(\left[ -1, 1 \right]\) for the current piecewise piece of our polynomial fit.
[ ]:
# Our first range
k = 5
x1 = LinRange(-1, 1, k)
y1 = sin.(x1)
scatter(x1, y1, label="data, subinterval 1")
# The first fit
A = vander_chebyshev(x1)
A_inv = inv(A)
@show c1 = A_inv * y1
x = LinRange(-1, 1, 100)
plot!(x, vander_chebyshev(x, k) * c1, label="fit, subinterval 1")
# Our second range
x2 = LinRange(1, 3, k)
y2 = sin.(x2)
scatter!(x2, y2, label="data, subinterval 2")
# The second fit
@show c2 = A_inv * y2
x = LinRange(-1, 1, 100)
plot!(x .+ 2, vander_chebyshev(x, k) * c2, label="fit, subinterval 2")
We could create a single large matrix equation representing this computation, but that matrix would have multiple copies of our Vandermonde matrix, which would be inefficient. Let’s instead make a single function that solves these systems of equations for each piece of our piecewise polynomial.
[ ]:
# Fit piecewise Chebyshev polynomials
function fit_piecewise_chebyshev(x, y, n, k)
"""Return the polynomial coefficents for a piecewise fit of Chebyshev polynomials
at the points x with the values f(x_i) = y_i. on subintervals of length w.
The parameter n is the number of subintervals.
The parameter k is the number of points in each subinterval.
The subinterval length is given by (x[end] - x[1]) / n.
This function returns (a, C), where a is the start of the first subinterval, and C
is a list of sets of Chebyshev coefficients, where C_j corresponds to the interval
[a + w*(j - 1), a + w*j].
"""
# Get dimensions
a = x[1]
w = (x[end] - x[1]) / n # length of a subinterval
# Build our Vandermonde matrix
A = vander_chebyshev((x[1:k] .- a) .* (2 / w) .- 1) # Why this ".- a) .* (2 / w) .- 1" part?
A_inv = inv(A)
# And a place to hold coefficients
C = []
for j in 1:n
# Now solve for the current set of coeffients for the current subset of y
yj = zeros(k) # What is the correct indices from y?
### BEGIN SOLUTION
### END SOLUTION
push!(C, c)
end
a, C
end
# Let's make the same data combined
k = 5
n = 2
x = LinRange(-1, -1 + 2*n, n*(k-1) + 1)
y = sin.(x)
p = scatter(x, y, label="data")
display(p)
# And let's try to fit our piecewise polynomial
a, C = fit_piecewise_chebyshev(x, y, n, k)
@show C;
[ ]:
k = 5
n = 5
x = LinRange(1, 1 + n*2, n*(k-1) + 1)
y = cos.(x) .* log.(x)
a, C = fit_piecewise_chebyshev(x, y, n, k)
@test a == 1
@test length(C) == 5
@test C[1] ≈ [-0.4164026546438594, -0.6167932864672474, -0.12767897738782175, 0.07298432524001167, 0.00027267080444547744]
@test C[2] ≈ [-0.6019173600399995, 0.8000865534138122, 0.29530092118109597, -0.028009254382361723, -0.008924184564117534]
@test C[3] ≈ [1.332621217129435, 0.5434549235122322, -0.379306420655918, -0.038210237255288715, 0.00846656539185615]
@test C[4] ≈ [-0.2857915440827067, -1.8185885681934457, 0.017546468197354503, 0.08409664184838822, 0.000779197662611586]
@test C[5] ≈ [-1.4529213285142264, 1.0464432985128442, 0.4681804363327252, -0.040158216086755155, -0.01093182996020823]
Applying our fit¶
Naturally, having the sets of polynomial coefficients is not very useful if we cannot use them. Now we need a function to apply our coefficients to produce values at a given set of points, \(x\).
In this function, we will assume that the subintervals overlap. That means that if we have 3 subintervals with 5 points each, then
subinterval 1 contains points 1, 2, 3, 4, and 5
subinterval 2 contains points 5, 6, 7, 9, and 9
subinterval 3 contains points 9, 10, 11, 12, and 13
These subintervals overlap, so we will have an ‘extra’ point on each subinterval.
We will call the number of points we are providing interpolated values for \(k_\text{interp}\) while we call the number of original points we fit against \(k_\text{fit}\). Thus, we will use a Vandermonde matrix of size \(k_\text{interp} \times k_\text{fit}\) to compute the interpolated values \(y_i\) on each subinterval.
We will therefore compute
for subintervals \(j \in \left\lbrace 1, 2, \dots \right\rbrace\) and \(w\) is the width of the subintervals.
Notice that like above, we shift and rescale the values of \(x\) for each subinterval \(j\) so that \(\left[ a + w \left( j - 1 \right), a + w j \right]\) becomes \(\left[ -1, 1 \right]\).
[ ]:
# Interpolate piecewise Chebyshev polynomials
function interpolate_piecewise_chebyshev(x, a, C, n, k_interp)
"""Return the interpolation from a piecewise fit of Chebyshev polynomials on the interval x.
n is a parameter sating how many subintervals there are.
k_interp is a parameter stating how many points are on each subinterval.
The width of the subintervals is therefore (x[end] - x[1]) / n.
This function returns y, where a is the start of the first sub interval, and C
is a list of sets of Chebyshev coefficients, where C_i corresponds to the interval
[a + 2*(i - 1), a + 2*i].
"""
# Get dimensions
k_fit = length(C[1])
@assert length(C) == n # these need to match
w = (x[end] - x[1]) / n
# Build our Vandermonde matrix
A = vander_chebyshev((x[1:k_interp] .- a) .* (2 / w) .- 1, k_fit)
# And a place to hold output
y = zeros(length(x))
for j in 1:n
# Now solve for the current set of coeffients for the current subset of y
yj = zeros(k_interp)
### BEGIN SOLUTION
### END SOLUTION
y[1 + (j - 1)*(k_interp - 1) : 1 + j*(k_interp - 1)] = yj
end
y
end
# Let's make the same data combined
k = 5
n = 2
x = LinRange(-1, -1 + 2*n, n*(k-1) + 1)
y = sin.(x)
scatter(x, y, label="data")
# And let's try to fit our piecewise polynomial
a, C = fit_piecewise_chebyshev(x, y, n, k)
# Finally, let's use that fit
k = 50
x = LinRange(-1, 3, 2*(k-1) + 1)
y = interpolate_piecewise_chebyshev(x, a, C, n, k)
p = plot!(x, y, label="fit")
[ ]:
k = 4
n = 2
x = LinRange(-1, -1+2*n, n*(k-1) + 1)
y = x.^3
a, C = fit_piecewise_chebyshev(x, y, n, k)
k = 50
x = LinRange(-1, -1+2*n, n*(k-1) + 1)
y = interpolate_piecewise_chebyshev(x, a, C, n, k)
@test y ≈ x.^3
Runge returns¶
Now let’s compare a fit with a piecewise polynomial to a single global polynomial.
You should notice that the Runge phenomenon is less pronounced but still present for each piece of the piecewise polynomial.
[ ]:
# Let's try a longer range
k = 7
n = 3
num_fit = n*(k-1) + 1
x = LinRange(-1, -1+2*n, num_fit)
runge(x) = 1 / (1 + 10*x^2)
y = runge.(x)
scatter(x, y, label="data")
# And let's try to fit our piecewise polynomial
a, C_piecewise = fit_piecewise_chebyshev(x, y, n, k)
# And a global fit
c_global = inv(vander_chebyshev(x)) * y
# Finally, let's use that fit
k = 50
x = LinRange(-1, -1+2*n, n*(k-1) + 1)
## Here's the true function
y_true = runge.(x)
plot!(x, y_true, label="true")
## The piecewise fit
y_piecewise = interpolate_piecewise_chebyshev(x, a, C_piecewise, n, k)
@show norm(y_true - y_piecewise)
plot!(x, y_piecewise, label="fit, piecewise")
## And the global fit
y_global = vander_chebyshev(x, num_fit) * c_global
@show norm(y_true - y_global);
plot!(x, y_global, label="fit, global")
Previously, we repaired this Runge phenomenon by selecting better interpolation points.
We will try this again here by using a cosine range on each subinterval. Fill in the function below so that our points \(x\) use the CosRange distribution on each subinterval.
[ ]:
# Here's our CosRange function
CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))
# And a piecewise range utility
function build_range(a, b, n, k, range_function)
"""Return a range of points for a piecewise interpolation on the interval [a, b].
There will be n subintervals with k points on each subinterval.
`range_function` is the function used to distribute points on one of the subintervals.
"""
w = (b - a) / n # Size of the subinterval
x1 = range_function(a, a + w, k)
x = zeros(n*(k - 1) + 1)
for j = 1:n
### BEGIN SOLUTION
### END SOLUTION
end
x
end
# And our true data to fit
k = 7
n = 3
num_fit = n*(k-1) + 1
x = build_range(-1, -1 + 2*n, n, k, CosRange)
y = runge.(x)
scatter(x, y, label="data")
# And let's try to fit our piecewise polynomial
a, C_piecewise = fit_piecewise_chebyshev(x, y, n, k)
# And a global fit
c_global = inv(vander_chebyshev(x)) * y
# Finally, let's use that fit
k = 50
x = LinRange(-1, -1+2*n, n*(k-1) + 1)
## Here's the true function
y_true = runge.(x)
plot!(x, y_true, label="true")
## The piecewise fit
y_piecewise = interpolate_piecewise_chebyshev(x, a, C_piecewise, n, k)
@show norm(y_true - y_piecewise)
plot!(x, y_piecewise, label="fit, piecewise")
## And the global fit
y_global = vander_chebyshev(x, num_fit) * c_global
@show norm(y_true - y_global);
plot!(x, y_global, label="fit, global")
Discuss the tradeoffs between using global and piecewise polynomial interpolation. Include any additional charts with other functions, numbers of subintervals (\(n\)), and number of points per subinterval (\(k\)) to support your discussion.
BEGIN SOLUTION¶
END SOLUTION¶
Integration¶
One important use of interpolation is integration. We have the polynomial
for the Chebyshev polynomials
Suppose we wanted to compute the integral
Using the powers of alegbraic simplification, we can instead compute
This feels like something we can get traction with, especially with the special properties of the Chebyshev polynomials.
where
We can use this recursion to turn our interpolation coefficients into integration coefficients.
[ ]:
# Function to get Chebyshev polynomial integrals
function chebyshev_integrals(k)
"""Return the values of the integrals of the Chebyshev polynomials on the interval [-1, 1].
"""
T_a = zeros(k+1)
T_b = zeros(k+1)
# Initalize
a = -1
b = 1
T_a[1] = 1 # Note that T_a[1] is T_0
T_a[2] = a
T_b[1] = 1
T_b[2] = b
# Follow the Chebyshev polynomial recursion to compute T_n
for i = 3:k+1
### BEGIN SOLUTION
### END SOLUTION
end
# And now compute the integrals
interp_T = zeros(k)
interp_T[1] = b - a
interp_T[2] = (1/2) * b^2 - (1/2) * a^2
for i = 3:k
# NOTE: Be careful with indexing here!
# Julia uses 1 based indexing but the math equations use 0 based indexing!!
### BEGIN SOLUTION
### END SOLUTION
end
interp_T'
end
# Let's try it
a = -1
b = 1
k = 5
x = LinRange(a, b, k)
y = x.^4 # We should be able to exactly integrate this!
p = scatter(x, y, label="true")
A = vander_chebyshev(x)
# Fit the Chebyshev polynomials
c = inv(A) * y
x_plt = LinRange(a, b, 100)
p = plot!(x_plt, vander_chebyshev(x_plt, k) * c, label="fit")
display(p)
# Integrate the fit polynomial
interp_T = chebyshev_integrals(k)
@show interp_T * c ≈ (1/5) * b^5 - (1/5) * a^5;
Piecewise polynomial integration¶
And now we can combine these two ideas. Complete the function below to compute an integral of a piecewise polynomial fit.
[ ]:
# Integrate using piecewise Chebyshev polynomials
function integrate_piecewise_chebyshev(x, y, n, k)
"""Return the value of integrating a piecewise fit of Chebyshev polynomials on the interval x for the values y.
n is a parameter stating how many subintervals there are.
k is a parameter stating how many points are on each subinterval.
This function returns the integral of the piecewise Chebyshev polynomial interpolation.
"""
# Get fit
a, C = fit_piecewise_chebyshev(x, y, n, k)
# Build our integral values
w = (x[end] - x[1]) / n
interp_T = chebyshev_integrals(k)
interp_T *= (w / 2) # Why do I need this line?
# And a place to hold output
interp = 0
for j in 1:n
# Sum up the integrals on each interval
### BEGIN SOLUTION
### END SOLUTION
end
interp
end
Let’s use our new function to compute the integral
[ ]:
# Let's integrate cosine with increasingly many subintervals and points
plot()
a = 0
b = 5
max_n = 5
max_k = 5
for n in 1:max_n # number of subintervals
println("---- $n subintervals ----")
for k in 2:max_k # number of points per subinterval
x = build_range(a, b, n, k, LinRange)
y = cos.(x)
# Interpolate the fit
a, C = fit_piecewise_chebyshev(x, y, n, k)
# Integrate the fit
integral = integrate_piecewise_chebyshev(x, y, n, k)
# Plot the fix
x = build_range(a, b, n, 100, LinRange)
y = interpolate_piecewise_chebyshev(x, a, C, n, 100)
if k == max_k
plot!(x, y, label="n = $n, k = $k")
end
# And print out info
println("$k points each, integral ≈ $integral")
end
println()
end
# And the true fit
println("---- true value ----")
@show sin(b) - sin(a)
x = LinRange(a, b, 300)
y = cos.(x)
p = plot!(x, y, color=:black, label="true");
[ ]:
display(p)
Exploration¶
Run some experiments using different values of \(n\) and \(k\) and using
LinRangeandCosRangefor our integration example above. Generate a plot summarizing your results and discuss your plot.
[ ]:
### BEGIN SOLUTION
### END SOLUTION
This technique can exactly integrate polynomials up to degree \(k - 1\). If we correctly select the location of the integration points, then we can see accurate integration for polynomials of up to degree \(2k - 1\). This is known as Gaussian quadrature. Discuss the pros and cons of this approach vs the one in this assignment.
[ ]:
### BEGIN SOLUTION
### END SOLUTION
Select a Julia package that implements numerical integration. Use the package for a few integrals. Discuss the ease of use and quality of documentation for the package you picked.
[ ]:
### BEGIN SOLUTION
### END SOLUTION