##### Document information

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 differential 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

Y ′(x0) Y ′(x1)

... Y ′(xn−1) Y ′(xn)

= 3

− 1 h20

1 h20

− 1 h20

1 h20 − 1

h21

1 h21

. . . . . . . . .

− 1 h2n−2

1 h2n−2

− 1 h2n−1

1 h2n−1

− 1 h2n−1

1 h2n−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 simplifies by multiplying by h. Note that all of the

1

diagonal elements in the interior of the second coefficient matrix become 0. Multiplying by the resulting matrix of coefficients of Y (x),

2 1 1 4 1

. . . . . . . . .

1 4 1 1 2

Y ′(x0) Y ′(x1)

... Y ′(xn−1) Y ′(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 coefficient matrices

First, identify the values of the function,

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

Substituting these values into (1),2 11 4 1 1 2

Y ′(x0)Y ′(x1) Y ′(x2)

= 3 1

10 −1

Solving the system yields, Y ′(x0)Y ′(x1)

Y ′(x2)

= 320 −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 Y ′′(x0) = Y

′′(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 different cubic poly- nomials: one on (0, 1) and the other on (1, 2). The general form of a cubic polynomial with an inflection 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 y′(1) = a− 3b = 0

Solving for the coefficients,

y(x) = 3

2 x− 1

2 x3

Hence the first 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 Infinite spline

Define some point x0 as the origin of the coordinates. In this case, the coefficient matrix is still tri-diagonal. However, the matrix is now infinite. The calculation is essentially the same as that in the finite case neglecting the endpoints of the interval,

. . . . . . . . .

1 4 1 1 4 1

1 4 1 . . . . . . . . .

...

Y ′(x−1) Y ′(x0) Y ′(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 find a solution for the points xk, k ≥ 0, we can construct the solution for k < 0 by symmetry.

Using this symmetry,

Y ′(x0) = 0

Y ′(xk) = −Y ′(x−k)∀k 6= 0

For k = 1,

4Y ′(x1) + Y ′(x2) = −3 (2)

For k > 1,

Y ′(xk−1) + 4Y ′(xk) + Y

′(xk+1) = 0 (3)

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

2.3.2 Solving difference equations

Elementary solution Solutions to iterative difference equations such as (3) can be solved by assuming a solution of a special type. The elementary solution to a difference equation of this type is

Y ′(xk) = r k

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,

Y ′±(xk) = ( −2±

√ 3 )k

4

Particular solution Because the difference equation is linear, any linear combination of these two solutions is also a solution. The general solution is of the form,

Y ′(xk) = a ( −2 +

√ 3 )k

+ b ( −2−

√ 3 )k

The absolute value of the second term is greater than 0. It will oscillate with increasing magnitude as k → ∞. To avoid such undesired oscillatory behavior, we must choose b = 0.

Now use (2) to solve for the coefficient a,

a ( 4r + r2

) = −3

Taking r = −2 + √ 3,

r2 + 4r = −1 ⇒ a = 3

The final form of the solution is

Y ′(xk) = 3 ( −2 +

√ 3 )k

Note. • This solution alternates in sign, as expected.

• The decay in the derivative is fast, which is also desirable,∣∣∣∣gk+1gk ∣∣∣∣ = −2 +√3 ≈ −26795

5