### 2.24.27 Natural splines: `spline`

#### Definition

Let σn be a subdivision of a real interval [a,b] :

a=x0, x1, ..., xn=b

s is a spline function of degree l, if s is an application from [a,b] in ℝ such that :

- s has continuous derivatives up to the order l−1,
- on each interval of the subdivision, s is a polynomial of degree less or equal than l.

#### Theorem

The set of spline functions of degree l on σn is a ℝ-vectorial subspace of dimension n+l.

Proof
On [a,x1], s is a polynomial A of degree less or equal to l, hence on [a,x1], s=A(x)=a0+a1x+...alxl and A is a linear combinaison of 1,x,...xl.
On [x1,x2], s is a polynomial B of degree less or equal to l, hence on [x1,x2], s=B(x)=b0+b1x+...blxl.
s has continuous derivatives up to order l−1, hence :

∀ 0 ≤ j ≤ l−1, B(j)(x1)−A(j)(x1)=0

therefore B(x)−A(x)=α1(x−x1)l or B(x)=A(x)+α1(x−x1)l.
Define the function :

q1(x) =

⎧
⎨
⎩

0 | sur | [a,x1]
(x−x1)l | sur | [x1,b]

Hence :

s|[a,x2]=a0+a1x+...alxl+α1q1(x)

On [x2,x3], s is a polynomial C of degree less or equal than l, hence on [x2,x3], s=C(x)=c0+c1x+...clxl.
s has continuous derivatives until l−1, hence :

∀ 0 ≤ j ≤ l−1, C(j)(x2)−B(j)(x2)=0

therefore C(x)−B(x)=α2(x−x2)l or C(x)=B(x)+α2(x−x2)l.
Define the function :

q2(x) =

⎧
⎨
⎩

0 | on | [a,x2]
(x−x2)l | on | [x2,b]

Hence : s|[a,x3]=a0+a1x+...alxl+α1q1(x)+α2q2(x)
And so on, the functions are defined by :

∀ 1 ≤ j ≤ n−1, qj(x) =

⎧
⎨
⎩

0 | on | [a,xj]
(x−xj)l | on | [xj,b]

hence,

s|[a,b]=a0+a1x+...alxl+α1q1(x)+....+αn−1qn−1(x)

and s is a linear combination of n+l independant functions 1,x,..xl,q1,..qn−1.

#### Interpolation with spline functions

If we want to interpolate a function f on σn by a spline function s of degree l, then s must verify s(xk)=yk=f(xk) for all 0≥ k≥ n. Hence there are n+1 conditions, and l−1 degrees of liberty. We can therefore add l−1 conditions, these conditions are on the derivatives of s at a and b.

Hermite interpolation, natural interpolation and periodic interpolation are three kinds of interpolation obtained by specifying three kinds of constraints. The unicity of the solution of the interpolation problem can be proved for each kind of constraints.

If l is odd (l=2m−1), there are 2m−2 degrees of liberty. The constraints are defined by :

- Hermite interpolation

  ∀ 1≤ j≤ m−1, s(j)(a)=f(j)(a), s(j)(b)=f(j)(b)

- Natural interpolation

  ∀ m ≤ j ≤ 2m−2, s(j)(a)=s(j)(b)=0

- periodic interpolation

  ∀ 1≤ j≤ 2m−2, s(j)(a)=s(j)(b)

If l is even (l=2m), there are 2m−1 degrees of liberty. The constraints are defined by :

- Hermite interpolation

  ∀ 1≤ j≤ m−1, s(j)(a)=f(j)(a), s(j)(b)=f(j)(b)

  and

  s(m)(a)=f(m)(a)

- Natural interpolation

  ∀ m ≤ j ≤ 2m−2, s(j)(a)=s(j)(b)=0

  and

  s(2m−1)(a)=0

- Periodic interpolation

  ∀ 1≤ j≤ 2m−1, s(j)(a)=s(j)(b)

A natural spline is a spline function which verifies the natural interpolation constraints.

`spline` takes as arguments a list of abscissa (by increasing order), a list of ordinates, a variable name, and a degree.
`spline` returns the natural spline function (with the specified degree and crossing points) as a list of polynomials, each polynomial being valid on an interval.

Examples:

1. a natural spline of degree 3, crossing through the points x0=0,y0=1, x1=1,y1=3 and x2=2, y2=0, input :
   `spline([0,1,2],[1,3,0],x,3)`
   Output is a list of two polynomial expressions of x :

   [ −5*x3/4+13*x/4+1, 5*(x−1)3/4−15*(x−1)2/4+(x−1)/−2+3 ]

   defined respectivly on the intervals [0,1] and [1,2].
2. a natural spline of degree 4, crossing through the points x0=0,y0=1, x1=1,y1=3, x2=2, y2=0 and x3=3, y3=−1, input :
   `spline([0,1,2,3],[1,3,0,-1],x,4)`
   Output is a list of three polynomial functions of x :

   [(−62*x4+304*x)/121+1,

   (201*(x−1)4−248*(x−1)3−372*(x−1)2+56*(x−1))/121+3,

   (−139*(x−2)4+556*(x−2)3+90*(x−2)2+−628*(x−2))/121]

   defined respectivly on the intervals [0,1], [1,2] and [2,3].
3. The natural spline interpolation of cos on [0,π/2,3π/2], input :
   `spline([0,pi/2,3*pi/2],cos([0,pi/2,3*pi/2]),x,3)`
   Output :

   [((3*π3+(−7*π2)*x+4*x3)*1/3)/(π3),

   ((15*π3+(−46*π2)*x+36*x2−8*x3)*1/12)/(π3)]