numeric::cubicSpline
--
interpolation by cubic splinesnumeric::cubicSpline
([x[0],y[0]], [x[1],y[1]],
..)
returns the cubic spline interpoland through a sequence of
coordinate pairs [x[i],y[i]].
numeric::cubicSpline([x[0],y[0]], .. , [x[n],y[n]]
<, BoundaryCondition> <, Symbolic> )
x[0],x[1],..,x[n] |
- | numerical real values in ascending order |
y[0],y[1],..,y[n] |
- | arbitrary expressions |
|
- | 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. |
the spline interpoland: a MuPAD procedure.
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.
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).numeric::cubicSpline
reorders the abscissae internally, issuing a warning.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.
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!numeric::cubicSpline
. This ordering is not checked,
even if the abscissae are numerical!BoundaryCondition
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).S
satisfying
S''(x[0])=S''(x[n])=0.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.[a,b]
produces a spline function
S
satisfying S'(x[0])=a, S'(x[n])=b.
Symbolic data a,b
are accepted.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:
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:
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:
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:
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:
spline