## `numeric::cubicSpline` -- interpolation by cubic splines

#### Introduction

`numeric::cubicSpline``([x[0],y[0]], [x[1],y[1]], ..)` returns the cubic spline interpoland through a sequence of coordinate pairs [x[i],y[i]].

#### Call(s)

`numeric::cubicSpline([x[0],y[0]], .. , [x[n],y[n]] <, BoundaryCondition> <, Symbolic> )`

#### Parameters

`x[0],x[1],..,x[n]` | - | numerical real values in ascending order
`y[0],y[1],..,y[n]` | - | arbitrary expressions

#### Options

``BoundaryCondition`` | - | the type of the boundary condition: either NotAKnot, Natural, Periodic, or Complete=`[a,b]` with arbitrary arithmetical expressions `a,b`.
`Symbolic` | - | prevents conversion of the input data to floating point numbers. With this option symbolic abscissae x[i] are accepted, which are assumed to be ordered.

#### Returns

the spline interpoland: a MuPAD procedure.

#### Related Functions

`numeric::lagrange`

#### Details

- The call

  `S:=numeric::cubicSpline([x[0],y[0]],..,[x[n],y[n]] <, Option>)`

  yields the cubic spline function `S` interpolating the data [x[0],y[0]],..,[x[n],y[n]], i.e, S(x[i])=y[i] for i=0..n. The spline function is a piecewise polynomial of degree <=3 on the intervals (-infinity, x[1]], [x[1],x[2]], .., [x[n-1], infinity) . `S` and its first two derivatives `S',S''` are continuous at the points x[1],..,x[n-1]. Note that `S` extends the polynomial representation on [x[0],x[1]], [x[n-1],x[n]] to (-infinity,x[1]] and [x[n-1],infinity), respectively.

- By default, NotAKnot boundary conditions are assumed, i.e., the third derivative `S'''` is continuous at the points x[1] and x[n-1]. With this boundary condition `S` is a polynomial on the intervals (-infty,x[2]] and [x_n-2,-infty).
- By default, all input data are converted to floating point numbers. This conversion may be suppressed by the option Symbolic.
- Without the option Symbolic the abscissae x[i] must be numerical real values in ascending order. If these data are not ordered, then `numeric::cubicSpline` reorders the abscissae internally, issuing a warning.
- The function `S` returned by `numeric::cubicSpline` may be called with one or two arguments.

  The call `S(z)` returns an explicit expression or a number, if `z` is a real number. Otherwise the unevaluated call `S(z)` is returned.

  The call `S(z,i)` is meant for symbolic arguments `z`. The argument `i` must be an integer. Internally, `z` is assumed to satisfy x[i]<= z <= x[i+1] and `S(z,i)` returns a polynomial expression in `z` representing the spline function on this interval.

- If `S` is generated with symbolic abscissae x[i] (necessarily using the option Symbolic), then the call `S(z)` with numerical `z` leads to an error. The call `S(z,i)` must be used for symbolic abscissae!

#### Option: Symbolic

- With this option no conversion of the input data to floating point numbers occurs.
- Symbolic abscissae x[i] are accepted.
- The ordering x[0] < x[1] < .. < x[n] is assumed by `numeric::cubicSpline`. This ordering is not checked, even if the abscissae are numerical!

#### Option: `BoundaryCondition`

- With the default boundary condition NotAKnot the third derivative `S'''` of the spline function is continuous at the points x[1] and x[n-1]. With this boundary condition `S` is a polynomial on the intervals (-infty,x[2]] and [x_n-2,-infty).
- The boundary condition Natural produces a spline function `S` satisfying S''(x[0])=S''(x[n])=0.
- The boundary condition Periodic produces a spline function `S` satisfying S(x[0])=S(x[n]),S'(x[0])=S'(x[n]), S''(x[0])=S''(x[n]). With this option the input data y[0],y[n] must coincide, otherwise an error occurs.
- The boundary condition Complete=`[a,b]` produces a spline function `S` satisfying S'(x[0])=a, S'(x[n])=b. Symbolic data `a,b` are accepted.

#### Example 1

We demonstrate some calls with numerical input data:

    >> data := [i, sin(i*PI/20)] $ i= 0..40:

    >> S1 := numeric::cubicSpline(data):

    >> S2 := numeric::cubicSpline(data, Natural):

    >> S3 := numeric::cubicSpline(data, Periodic):

    >> S4 := numeric::cubicSpline(data, Complete = [3, PI]):

At the abscissae the corresponding input data are reproduced:

    >> float(op(data, 6)[2]), S1(5), S2(5), S3(5), S4(5)

          0.7071067812, 0.7071067812, 0.7071067812, 0.7071067812,

             0.7071067812

Interpolation between the abscissae depends on the boundary condition:

    >> S1(4.5), S2(4.5), S3(4.5), S4(4.5)

              0.6494470263, 0.6494470123, 0.6494469992, 0.6517696766

These are the cubic polynomials in `z` defining the spline on the interval x[0]=0 <= z <= x[1]=1:

    >> expand(S1(z, 0)), expand(S2(z, 0)), expand(S3(z, 0)),
       expand(S4(z, 0))

                                             2                    3
          0.1570962007 z - 0.00002961951081 z  - 0.0006321161139 z ,

                                               3
             0.1570790998 z - 0.0006446347923 z ,

                                               2                    3
             0.157063067 z + 0.00002776961744 z  - 0.0006563716136 z ,

                                  2                3
             3.0 z - 4.924083441 z  + 2.080517906 z

    >> delete data, S1, S2, S3, S4:

#### Example 2

We demonstrate some calls with symbolic data:

    >> S := numeric::cubicSpline([i, y.i] $ i=0..3):

    >> S(1/2)

                   0.3125 y0 + 0.9375 y1 - 0.3125 y2 + 0.0625 y3

This is the cubic polynomial in `z` defining the spline on the interval x[0]=0 <= z <= x[1]=1:

    >> S(z, 0)

          y0 + z (3.0 y1 - 1.833333333 y0 - 1.5 y2 + 0.3333333333 y3 + z

             (1.0 y0 - 2.5 y1 + 2.0 y2 - 0.5 y3 +

             z (0.5 y1 - 0.1666666667 y0 - 0.5 y2 + 0.1666666667 y3)))

With the option Symbolic exact arithmetic is used:

    >> S := numeric::cubicSpline([i, y.i] $ i=0..3, Symbolic):

    >> S(1/2)

                             5 y0   15 y1   5 y2   y3
                             ---- + ----- - ---- + --
                              16     16      16    16

Also symbolic boundary data are accepted:

    >> S := numeric::cubicSpline([i, exp(i)] $ i=0..10,
                                 Complete = [a, b]):

    >> S(0.1)

               0.08341154273 a + 0.00000005947817812 b + 1.020064753

    >> S := numeric::cubicSpline([0, y0], [1, y1], [2, y2],
                                 Symbolic, Complete=[a, 5]):

    >> collect(S(z, 0), z)

                      3 / 3 a   5 y0          3 y2       \
          y0 + a z + z  | --- + ---- - 2 y1 + ---- - 5/4 | +
                        \  4     4             4         /

              2 /        9 y0   7 a   3 y2       \
             z  | 3 y1 - ---- - --- - ---- + 5/4 |
                \         4      4     4         /

    >> delete S:

#### Example 3

We demonstrate the use of symbolic abscissae. Here the option Symbolic is mandatory.

    >> S := numeric::cubicSpline([x.i, y.i] $ i=0..2, Symbolic):

The spline function `S` can only be called with 2 arguments. This is the cubic polynomial in `z` defining the spline on the interval x[0] <= z <= x[1]:

    >> S(z, 0)

                        / (y1 - y0) (x1 - 2 x0 + x2)
          y0 + (z - x0) | -------------------------- -
                        \    (x1 - x0) (x2 - x0)

             (x1 - x0) (y2 - y1)
             ------------------- + (z - x0)
             (x2 - x0) (x2 - x1)

             /       y2 - y1               y1 - y0       \ \
             | ------------------- - ------------------- | |
             \ (x2 - x0) (x2 - x1)   (x1 - x0) (x2 - x0) / /

    >> delete S:

#### Example 4

Spline functions can be plotted:

    >> S := numeric::cubicSpline([i, 1/(1 + i^2/100)] $ i=0..100):

    >> plotfunc2d(S(x), x = 0..100, Ticks = [10, Steps = 0.2])

    >> delete S:

#### Example 5

We demonstrate how to generate a phase plot of the differential equation x''(t)+x(t)^3 = sin(t), with initial conditions x(0)=x'(0)=0. First, we use `numeric::odesolve` to compute a numerical mesh of solution points [x[i],y[i]]=[x(t[i]),x'(t[i])] with n+1 equidistant time nodes t[0],..,t[n] in the interval [0,20]:

    >> DIGITS := 4: n := 100:

    >> for i from 0 to n do t[i] := 20/n*i: end_for:

    >> f := (t, x) -> [x[2], sin(t) - x[1]^3]:

    >> x[0] := 0: y[0] := 0:

    >> for i from 1 to n do
         [x[i], y[i]] :=
         numeric::odesolve(t[i-1]..t[i], f, [x[i-1], y[i-1]]):
       end_for:

The mesh of the (x(t),x'(t)) phase plot consists of the following points:

    >> Plotpoints := [point(x[i], y[i]) $ i=0..n]:

We wish to connect these points by a spline curve. We define a spline interpoland Sx(t) approximating the solution x(t) by interpolating the data [t[0],x[0]], .. , [t[n], x[n]]. A spline interpoland Sy(t) approximating x'(t) is obtained by interpolating the data [t[0],y[0]], .. , [t[n], y[n]]:

    >> Sx := numeric::cubicSpline([t[i], x[i]] $ i=0..n):

    >> Sy := numeric::cubicSpline([t[i], y[i]] $ i=0..n):

Finally, we plot the mesh points together with the interpolating spline curve:

    >> plot2d([Mode = List, Plotpoints, PointWidth = 30],
              [Mode = Curve,[Sx(z), Sy(z)], z = [0, 20], Grid = [5*n]])

The function `plot::ode` serves for displaying numerical solutions of ODEs. In fact, it is implemented as indicated by the previous commands. The following call produces the same plot:

    >> plot(plot::ode(
        [t[i] $ i=0..n], f, [x[0], y[0]],
        [(t, x) -> [x[1], x[2]], Style = Points, Color = RGB::Red],
        [(t, x) -> [x[1], x[2]], Style = Splines, Color = RGB::Blue])):

    >> delete DIGITS, n, i, t, f, x, y, Plotpoints, Sx, Sy:

#### Changes

- used to be called `spline`
- The routine was completely redesigned: the functionality was extended, performance improved.

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000