Worksheet 5: Spline Interpolation Solutions

(1)

Physik auf dem Computer SS 2017

Worksheet 5: Spline Interpolation Solutions

June 13, 2017

Author: Johannes Zeman Task 5.1: Cubic and Quadratic Splines (6 points)

In this task, you will spline-interpolate the following functions on the specified domains:

Name Definition Domain

Sine Function f (x) = sin x [0, 2π]

Runge Function g(x) = _1+x ¹

2

[−5, 5]

Lennard-Jones Function h(x) = x ⁻¹² − x ⁻⁶ [1, 5]

The python script ws5.py provides a demo which uses the class scipy.interpolate.interp1d to compute interpolating splines.

• 5.1.1 (2 points) Write a class CubicSplineInterpolation which implements cubic spline interpolation. As in the previous worksheet, the class shall provide a constructor init(self, ...) to initialize the interpolation and a method call(self,x) to compute the value of the interpolating function at x . On the boundaries, set the second derivative of the spline to zero.

Use the class to create the same plots as the demo does. Note that the splines in this class will not be identical to the splines in the demo as SciPy uses different boundary conditions.

Hint The constructor has to compute the spline coefficients by first generating the defining linear equations (see lecture notes pp. 35-36, especially the equation after equation 2.8) and then solving them using scipy.linalg.solve .

• 5.1.2 (2 points) Derive the equations for the quadratic spline, where the spline polynomial is defined as

S _j (x) = a _j + b _j (x − x _j ) + c _j (x − x _j ) ² . (1) The conditions for the quadratic spline are:

– It matches the function values at the supporting points x _j and x _j+1 . – The first derivative of the spline at x _j+1 is the same for S _j and S _j+1 .

– In your implementation, the first derivative should vanish at the lower boundary.

• 5.1.3 (2 points) Write a class QuadraticSplineInterpolation which implements the quadratic

spline interpolation. Redo the same plots as in task 5.1.1 with quadratic splines.

(2)

Solution of Task 5.1.1 (Cubic Spline): Derivation of Defining Linear Equations

A cubic spline interpolation connecting n points (x j |y _j ), j ∈ [0, 1, ...n − 1] can be expressed as splines S _j (x) of the general form

S j (x) = a j + b j (x − x j ) + c j (x − x j ) ² + d j (x − x j ) ³ , j ∈ [0, 1, ..., n − 1] , (2) where each segment S j (x) is the interpolating third-order polynomial on the interval [x j , x j+1 ].

The corresponding first and second derivatives with respect to x are thus d

dx S _j (x) = b _j + 2c _j (x − x _j ) + 3d _j (x − x _j ) ² (3) d ²

dx ² S _j (x) = 2c _j + 6d _j (x − x _j ) . (4)

Each segment has to fulfill the following conditions:

i) S j (x j ) = y j , j ∈ [0, 1, ..., n − 1]

ii) S _j+1 (x _j+1 ) = S _j (x _j+1 ) , j ∈ [0, 1, ..., n − 2]

iii) _dx ^d S j+1 (x j+1 ) = _dx ^d S j (x j+1 ) , j ∈ [0, 1, ..., n − 2]

iv) _dx ^d

²₂

S _j+1 (x _j+1 ) = _dx ^d

²₂

S _j (x _j+1 ) , j ∈ [0, 1, ..., n − 2]

Additionally, it has to fulfill the boundary conditions v) _dx ^d

²₂

S ₀ (x ₀ ) = A

vi) _dx ^d

²2

S n−2 (x n−1 ) = B .

Here, we impose the so-called natural boundary conditions by requiring A = B = 0.

Using these conditions, the defining linear equations for the coefficients a _j , b _j , c _j , and d _j are derived in the following.

For the coefficients a _j , it follows from condition i) together with Eq. (2) that

a _j = y _j . (5)

Applying condition iv) to Eq. (4) yields

2c _j+1 + 6d _j+1 (x _j+1 − x _j+1 )

| {z }

=0

= 2c _j + 6d _j (x _j+1 − x _j )

⇒ d j = c j+1 − c j

3(x _j+1 − x _j ) . (6)

Furthermore, applying conditions i) and ii) to Eq. (2) and inserting Eq. (5) and (6) yields y j+1 = y j + b j (x j+1 − x j ) + c j (x j+1 − x j ) ² + d j (x j+1 − x j ) ³

⇒ b _j = y _j+1 − y _j x _j+1 − x _j − 1

3 (x _j+1 − x _j )(2c _j + c _j+1 ) . (7)

Analogously, by applying condition iii) on Eq. (3) and inserting Eq. (6), it follows that b _j+1 = b _j + (c _j + c _j+1 )(x _j+1 − x _j ) .

2

(3)

Plugging in Eq. (7) leads to y j+2 − y j+1

x j+2 − x j+1

= y j+1 − y j

x j+1 − x j

+ 1

3 (c j (x j+1 − x j ) + 2c j+1 (x j+2 − x j ) + c j+2 (x j+2 − x j+1 ))

Inserting Eq. (7), bringing all c-dependent terms to the left side, dividing by (x j+2 − x j ) and shifting the index j → (j − 1) yields the system of linear equations

α j c j−1 + 2c j + β j c j+1 = 3γ j (8) with

α _j := x _j − x j−1

x j+1 − x j−1

, β _j := x _j+1 − x _j x j+1 − x j−1

, and γ _j :=

y

j+1

−y

j

x

j+1

−x

_j

− _x ^y

^j

^−y

^j−1

j

−x

j−1

x j+1 − x j−1

. In order to solve the equation system, we write Eq. (8) in matrix notation:







α ₁ 2 β ₁ 0

α ₂ 2 β ₂ . .. . .. . ..

0 α n−2 2 β n−2











 c ₁ c ₂ .. . c n−2







=





 3γ ₁ 3γ ₂ .. . 3γ n−2







(9)

Yet, there are two equations missing in order to solve Eq. (9), which are the boundary conditions v) and vi). With natural boundary conditions A = B = 0, applying these conditions to Eq. (4) together with Eq. (6) yields

c 0 = c n−1 = 0 .

Adding these equations to Eq. (9) finally leads to the complete system of linear equations:







1 0 0

α ₁ 2 β ₁ α 2 2 β 2

. .. . .. . ..

α n−2 2 β n−2

0 0 1











 c ₀ c ₁ c 2

.. . c n−2

c n−1







=





 0 3γ ₁ 3γ 2

.. . 3γ n−2

0 





(10)

Once the system of linear equations is solved, the coefficients b _j and d _j are obtained from Eq. (7) and (6), respectively.

For the implementation of the class CubicSplineInterpolation , see the Python script ws5_solution.py .

(4)

Solution of Task 5.1.2 (Quadratic Spline): Derivation of Defining Linear Equations A quadratic spline interpolation connecting n points (x j |y _j ), j ∈ [0, 1, ...n − 1] can be expressed as splines S _j (x) of the general form

S _j (x) = a _j + b _j (x − x _j ) + c _j (x − x _j ) ² , j ∈ [0, 1, ..., n − 1] , (11) where each segment S j (x) is the interpolating second-order polynomial on the interval [x j , x j+1 ].

The corresponding first derivatives with respect to x are thus d

dx S j (x) = b j + 2c j (x − x j ) . (12)

Each segment has to fulfill the following conditions:

i) S j (x j ) = y j , j ∈ [0, 1, ..., n − 1]

ii) S _j+1 (x _j+1 ) = S _j (x _j+1 ) , j ∈ [0, 1, ..., n − 2]

iii) _dx ^d S j+1 (x j+1 ) = _dx ^d S j (x j+1 ) , j ∈ [0, 1, ..., n − 2]

Additionally, it has to fulfill the boundary condition iv) _dx ^d S 0 (x 0 ) = A .

We impose natural boundary conditions by setting A = 0.

Using these conditions, the defining linear equations for the coefficients a j , b j , and c j are derived in the following.

For the coefficients a _j , it follows from condition i) together with Eq. (11) that

a _j = y _j . (13)

Applying condition iii) to Eq. (12) yields

b j+1 + 2c j+1 (x j+1 − x j+1 )

| {z }

=0

= b j + 2c j (x j+1 − x j )

⇒ c j = b j+1 − b j

2(x j+1 − x j ) . (14)

Furthermore, applying conditions i) and ii) to Eq. (11) and inserting Eq. (13) and (14) yields y j+1 = y j + b j (x j+1 − x j ) + c j (x j+1 − x j ) ²

⇒ b _j + b _j+1 = 2 y _j+1 − y _j

x _j+1 − x _j . (15)

With an index shift j → (j − 1) we write this as

b j−1 + b _j = 2γ _j (16)

with

γ _j := y _j − y j−1

x _j − x j−1

.

4

(5)

In order to solve the equation system, we write Eq. (16) in matrix notation:







1 1 0

1 1

. .. ...

0 1 1











 b 1

b 2

.. . b n−1







=





 2γ 1

2γ 2

.. . 2γ n−1







(17)

Yet, there is one equation missing in order to solve Eq. (17), which is the boundary condition iv).

With the natural boundary condition A = 0, applying this condition to Eq. (12) yields b ₀ = 0 .

Adding this equation to Eq. (17) finally leads to the complete system of linear equations:

Worksheet 5: Spline Interpolation Solutions

Physik auf dem Computer SS 2017