Matrix Formulation of Variational Method

Interpolation

Adrian Down

March 03, 2006

1 Matrix formulation

1.1 Unequally-spaced mesh points

Recall the matrix formulation of the diﬀerential equation considered last

lecture,

2

h0

1

h0

1

h0

2

h0+2

h1

1

h1

.........

1

hn−2

2

hn−2+2

hn−1

1

hn−1

1

hn−1

2

hn−1

Y0(x0)

Y0(x1)

.

.

.

Y0(xn−1)

Y0(xn)

= 3

−1

h2

0

1

h2

0

−1

h2

0

1

h2

0−1

h2

1

1

h2

1

.........

−1

h2

n−2

1

h2

n−2−1

h2

n−1

1

h2

n−1

−1

h2

n−1

1

h2

n−1

Y(x0)

Y(x1)

.

.

.

Y(xn−1)

Y(xn)

where hj=xj+1 −xj, for j∈ { 0, . . . , n −1}. This is general formulation

for non-uniform node spacing.

1.2 Equally-spaced mesh points

Suppose that the mesh points are uniformly spaced, so that hi=hj=

h, ∀i, j. The matrix simpliﬁes by multiplying by h. Note that all of the

1

diagonal elements in the interior of the second coeﬃcient matrix become 0.

Multiplying by the resulting matrix of coeﬃcients of Y(x),

2 1

1 4 1

.........

1 4 1

1 2

Y0(x0)

Y0(x1)

.

.

.

Y0(xn−1)

Y0(xn)

=3

h

−Y(x0) + Y(x1)

−Y(x0) + Y(x2)

.

.

.

−Y(xn−2) + Y(xn)

−Y(xn−1) + Y(xn)

(1)

2 Debugging examples

2.1 Three points, symmetric about the center

We take the three nodes to be

(0,0) (1,1) (2,0)

In this case, we expect the interpolating polynomials to take an arch shape.

2.1.1 Construct the coeﬃcient matrices

First, identify the values of the function,

Y(x0) = 0 Y(x1) = 1 Y(x2) = 0

Substituting these values into (1),

2 1

141

1 2

Y0(x0)

Y0(x1)

Y0(x2)

=3

1

1

0

−1

Solving the system yields,

Y0(x0)

Y0(x1)

Y0(x2)

=

3

2

0

−3

2

This result is consistent with what we would expect by symmetry.

2

2.1.2 Restrictions on the second derivative

The second derivative is clearly continuous at the point (1,1) by symmetry.

We must still enforce that Y00(x0) = Y00 (xn) = 0. We could write the explicit

form of the interpolating polynomial and compute the derivative directly.

However, we can analyze the second derivative in a more clever method.

Recall that our interpolation method relies on two diﬀerent cubic poly-

nomials: one on (0,1) and the other on (1,2). The general form of a cubic

polynomial with an inﬂection point at the origin is

y=ax −bx3

To be consistent with the other initial conditions, we must enforce that,

y(1) = a−b= 1

y0(1) = a−3b= 0

Solving for the coeﬃcients,

y(x) = 3

2x−1

2x3

Hence the ﬁrst derivative of the proposed solution Y(x) is consistent with

the general form of a polynomial having second derivative equal to 0 at the

origin.

Note. A similar argument could be conducted on the interpolating polyno-

mial on (1,2).

2.2 Inﬁnite spline

Deﬁne some point x0as the origin of the coordinates. In this case, the

coeﬃcient matrix is still tri-diagonal. However, the matrix is now inﬁnite.

The calculation is essentially the same as that in the ﬁnite case neglecting

the endpoints of the interval,

.........

1 4 1

141

1 4 1

.........

.

.

.

Y0(x−1)

Y0(x0)

Y0(x1)

.

.

.

=3

h

.

.

.

−Y(x−2) + Y(x0)

−Y(x−1) + Y(x1)

−Y(x0) + Y(x2)

.

.

.

3

2.3 Y(x0) = 1, all other f(xk) = 0

2.3.1 Setup

We expect a curve that looks like a symmetric sinusoid that decays away

from the origin. Thus if we can ﬁnd a solution for the points xk,k≥0, we

can construct the solution for k < 0 by symmetry.

Using this symmetry,

Y0(x0) = 0

Y0(xk) = −Y0(x−k)∀k6= 0

For k= 1,

4Y0(x1) + Y0(x2) = −3 (2)

For k > 1,

Y0(xk−1)+4Y0(xk) + Y0(xk+1) = 0 (3)

We also enforce the boundary condition that Y0(xk)→0 as k→ ∞.

2.3.2 Solving diﬀerence equations

Elementary solution Solutions to iterative diﬀerence equations such as

(3) can be solved by assuming a solution of a special type. The elementary

solution to a diﬀerence equation of this type is

Y0(xk) = rk

Substituting this assumed form of the solution into (3),

rk−1+ 4rk+rk+1 = 0

Dividing by rk−1,

1+4r+r2= 0

Using the quadratic formula,

r=−2±√3

Hence we can construct two linearly independent solutions,

Y0

±(xk) = −2±√3k

4