Interpolation as a linear system

The polynomial interpolation problem in the definition above seeks a polynomial of the form \[p(t) = c_1 + c_2 t + c_3 t^2 + \cdots + c_n t^{n-1}\] such that \(y_i = p(t_i)\) for all \(i\). We can rewrite this as

These equations can be written succinctly in our usual linear system form

\[\begin{bmatrix} 1 & t_1 & \cdots & t_1^{n-2} & t_1^{n-1} \\ 1 & t_2 & \cdots & t_2^{n-2} & t_2^{n-1} \\ 1 & t_3 & \cdots & t_3^{n-2} & t_3^{n-1} \\ \vdots & \vdots & & \vdots & \vdots \\ 1 & t_n & \cdots & t_n^{n-2} & t_n^{n-1} \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_n \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n \end{bmatrix},\]

which we denote \(\mathbf{V} \mathbf{c} = \mathbf{y}\).

We begin with an example of data about the population of China.

year = [1982, 2000, 2010, 2015] #years
pop = [1008.18, 1262.64, 1337.82, 1374.62] #population in millions

4-element Vector{Float64}:
 1008.18
 1262.64
 1337.82
 1374.62

t = year .- 1980.0 #indexing years after 1980
y = pop;

V = [t[i]^j for i=1:4, j=0:3]

4×4 Matrix{Float64}:
 1.0   2.0     4.0      8.0
 1.0  20.0   400.0   8000.0
 1.0  30.0   900.0  27000.0
 1.0  35.0  1225.0  42875.0

Now we can solve the system \[\mathbf{V}\mathbf{c} = \mathbf{y}.\]

c = V \ y

4-element Vector{Float64}:
 962.2387878787875
  24.127754689754774
  -0.5922620490620537
   0.00684386724386731

In the next couple of weeks, we’ll study the algorithms used under the hood when we use the Julia backslash operator. Even the algorithms solving a simple problem like a linear system are mathematically rich!

y - V*c

4-element Vector{Float64}:
 0.0
 0.0
 2.2737367544323206e-13
 0.0

Remember that this solution contains the coefficients of a polynomial that interpolates these four points!

using Polynomials
p = Polynomial(c)
p(2005-1980)

1302.2043001443

To evaluate the polynomial over this interval, we use the broadcasting operation . to evaluate it at every element of an array of points in the interval.

tt = range(0,35,length=500) #500 times in the interval [0,35]
yy = p.(tt)                 #evaluate the polynomial at all vector elements
foreach(println,yy[1:4])

962.2387878787875
963.9282039963299
965.6118068288089
967.2896105457506

Block matrix expressions

It is often helpful to break a given matrix into smaller, named submatrices. For instance, we might write that matrix broken into a \(2 \times 3\) array of six blocks as \[\begin{bmatrix} \mathbf{A}_{11} & \mathbf{A}_{12} & \mathbf{A}_{13} \\ \mathbf{A}_{21} & \mathbf{A}_{22} & \mathbf{A}_{23} \end{bmatrix}.\]

Vectors and matrices in Julia

Once you get the hang of vector and matrix manipulation in Julia, it tends to be fairly straightforward and intuitive. We’ll get started with some practice!

Get on the server at https://149-165-153-51.js2proxy.cacao.run/ and download our in class playground notebook using wget "https://raw.githubusercontent.com/jamiehadd/164-ScientificComputing/refs/heads/main/Lecture3/Lecture3_playground.ipynb". Run the included code blocks, understand the output, and note any questions you have!

The Big Lesson

You probably remember the process of calculating inverse matrices from Linear Algebra – challenging! It turns out calculating the inverse via numerical computation is also challenging and almost never the best tool at your disposal.

A = randn(10000,10000)
x_true = ones(10000,);
b = A*x_true;

@time x = A\b;

  3.749634 seconds (9 allocations: 763.092 MiB, 0.43% gc time)

@time x = inv(A)*b;

 13.757109 seconds (13 allocations: 767.975 MiB, 0.06% gc time)

One way to check a computed answer of a linear system is to compute the residual. The entries of the residual should be ideally close to machine precision (relative to the size of the entries in the solution).

A = [1 0 -1; 2 2 1; -1 -3 0]
b = [1,2,3]
x = A\b

3-element Vector{Float64}:
  2.1428571428571432
 -1.7142857142857144
  1.1428571428571428

residual = b - A*x

3-element Vector{Float64}:
 -4.440892098500626e-16
 -4.440892098500626e-16
  0.0

Triangular Systems

It is especially easy to solve a triangular system of linear equations – one defined with a triangular matrix.

This process is given by the formulas

\[\begin{align*} x_1 &= \frac{b_1}{L_{11}} \\x_2 &= \frac{b_2 - L_{21}x_1}{L_{22}} \\x_3 &= \frac{b_3 - L_{31}x_1 - L_{32}x_2}{L_{33}} \\x_4 &= \frac{b_4 - L_{41}x_1 - L_{42}x_2 - L_{43}x_3}{L_{44}}. \end{align*}\]

Let’s implement the forward and backward substitution processes.

"""
    forwardsub(L,b)

Solve the lower triangular lienar system with matrix 'L' and
right-hand side vector 'b'.
"""
function forwardsub(L,b)
    n = size(L,1)
    x = zeros(n)
    x[1] = b[1]/L[1,1]
    for i in 2:n
        s = sum( L[i,j]*x[j] for j in 1:i-1 )
        x[i] = ( b[i] - s ) / L[i,i]
    end
    return x
end

forwardsub

"""
    backsub(U,b)

Solve the upper triangular lienar system with matrix 'U' and
right-hand side vector 'b'.
"""
function backsub(U,b)
    n = size(U,1)
    x = zeros(n)
    x[n] = b[n]/U[n,n]
    for i in n-1:-1:1
        s = sum( U[i,j]*x[j] for j in i+1:n )
        x[i] = ( b[i] - s ) / U[i,i]
    end
    return x
end

backsub

Since we don’t know the solution to this system, to see how accurate the approximation is, we use the residual.

b - L*x

5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0

Now, we’ll consider another example where we know the solution due to how we built the system.

α = 0.3;
β = 1e12;
U = diagm( 0=>ones(5), 1=>[-1,-1,-1,-1] )
U[1,[4,5]] = [ α-β, β ]
U

5×5 Matrix{Float64}:
 1.0  -1.0   0.0  -1.0e12   1.0e12
 0.0   1.0  -1.0   0.0      0.0
 0.0   0.0   1.0  -1.0      0.0
 0.0   0.0   0.0   1.0     -1.0
 0.0   0.0   0.0   0.0      1.0

x_exact = ones(5)
b = [α,0,0,0,1]

5-element Vector{Float64}:
 0.3
 0.0
 0.0
 0.0
 1.0

x = backsub(U,b)
err_x = x - x_exact

5-element Vector{Float64}:
 -4.882812499995559e-5
  0.0
  0.0
  0.0
  0.0

We got only four digits of accuracy in the first entry of the solution, which is not great since we started with 16 digits of accuracy!