$SPAD/src/input spline.input

Timothy Daly December 3, 2016

Abstract

This is a problem picked up from a youtube lecture broadcast from the University of Southern Florida called “Quadratic Spline Interpolation”. [?]

Contents

2

— * —

)set break resume
)sys rm -f spline.output
)spool spline.output
)set message test on
)set message auto off
)clear all

———-

The upward velocity of a rocket is given as a function of time. Using quadratic splines

1. find the velocity at t=16 seconds

2. find the acceleration at t=16 seconds

3. find the distance covered between t=11 and t=16 seconds

t  v(t)
s  m/s
0  0
10 227.04
15 362.78
20 517.35
22.5 602.97
30 901.67

We are using quadratic splines to approximate the trajectory and estimate the values. There are 5 equations:

v(t) = a_{1}t^{2} + b_{1}t + c_{1} 0 ≤ t ≤ 10
v(t) = a_{2}t^{2} + b_{2}t + c_{2} 10 ≤ t ≤ 15
v(t) = a_{3}t^{2} + b_{3}t + c_{3} 15 ≤ t ≤ 20
v(t) = a_{4}t^{2} + b_{4}t + c_{4} 20 ≤ t ≤ 22.5
v(t) = a_{5}t^{2} + b_{5}t + c_{5} 22.5 ≤ t ≤ 30

Thus, there are 15 unknowns:

a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, b_{1}, b_{2}, b_{3}, b_{4}, b_{5}, c_{1}, c_{2}, c_{3}, c_{4}, c_{5}

So we need to develop 15 equations. We can get the first 10 equations because the splines need to go thru the endpoints. From the first two data points [0, 0] and [10, 227.04] we know that the first spline goes

through both data points:

v(t) = a_{1}t^{2} + b_{1}t + c_{1} 0 ≤ t ≤ 10
v(0) = a_{1}(0)^{2} + b_{1}(0) + c_{1} = 0

v(0) = a_{1}(10)^{2} + b_{1}(10) + c_{1} = 0

The next set of equations come from the points [10, 227.04] and [15, 362.78] since the second spline goes thru these points:

v(t) = a_{2}t^{2} + b_{2}t + c_{2} 10 ≤ t ≤ 15
v(10) = a_{2}(10)^{2} + b_{2}(10) + c_{2} = 227.04
v(15) = a_{2}(15)^{2} + b_{2}(15) + c_{2} = 362.78

The next set of equations come from the points [15, 362.78] and [20, 517.35] since the third spline goes thru these points:

v(t) = a_{3}t^{2} + b_{3}t + c_{3} 15 ≤ t ≤ 20
v(15) = a_{3}(15)^{2} + b_{3}(15) + c_{3} = 362.78
v(20) = a_{3}(20)^{2} + b_{3}(20) + c_{3} = 517.35

The next set of equations come from the points [20, 517.35] and [22.5, 602.97] since the fourth spline goes thru these points:

v(t) = a_{4}t^{2} + b_{4}t + c_{4} 20 ≤ t ≤ 22.5
v(20) = a_{4}(20)^{2} + b_{4}(20) + c_{4} = 517.35
v(22.5) = a_{4}(22.5)^{2} + b_{4}(22.5) + c_{4} = 602.97

The next set of equations come from the points [22.5, 602.97] and [30, 901.67] since the fifth spline goes thru these points:

v(t) = a_{5}t^{2} + b_{5}t + c_{5} 22.5 ≤ t ≤ 30
v(22.5) = a_{5}(22.5)^{2} + b_{5}(22.5) + c_{5} = 602.97
v(30) = a_{5}(30)^{2} + b_{5}(30) + c_{5} = 901.67

The next four equations come from the constraint that the slope of each of these curves has to match at the point where they touch. This means that they must have the same derivative at that point.

We skip the endpoint 0 since this is unconstrained. The point x = 10 has two equations that are valid at that point:

v(t) = a_{1}t^{2} + b_{1}t + c_{1} 0 ≤ t ≤ 10
w(t) = a_{2}t^{2} + b_{2}t + c_{2} 10 ≤ t ≤ 15

We compute the derivative of each equation, set them equal, and evaluate the result at the common point x = 10:

\frac{d}{dt} (a_{1}t^{2} + b_{1}t + c_{1})|_{t=10} = \frac{d}{dt} (a_{2}t^{2} + b_{2}t + c_{2})|_{t=10}

(2a_{1}t + b_{1})|_{t=10} = (2a_{2}t + b_{2})|_{t=10}

2a_{1}(10) + b_{1} = 2a_{2}(10) + b_{2}

20a_{1} + b_{1} − 20a_{2} − b_{2} = 0

We do this with the second pair:

\frac{d}{dt} (a_{2}t^{2} + b_{2}t + c_{2})|_{t=15} = \frac{d}{dt} (a_{3}t^{2} + b_{3}t + c_{3})|_{t=15}

(2a_{2}t + b_{2})|_{t=15} = (2a_{3}t + b_{3})|_{t=15}
2a_{2}(15) + b_{2} = 2a_{3}(15) + b_{3}
30a_{2} + b_{2} − 30a_{3} − b_{3} = 0

We do this with the third pair:

\frac{d}{dt} (a_{3}t^{2} + b_{3}t + c_{3})|_{t=20} = \frac{d}{dt} (a_{4}t^{2} + b_{4}t + c_{4})|_{t=20}

(2a_{3}t + b_{2})|_{t=20} = (2a_{4}t + b_{3})|_{t=20}
2a_{3}(20) + b_{3} = 2a_{4}(20) + b_{4}
40a_{3} + b_{3} − 40a_{4} − b_{4} = 0

We do this with the fourth pair:

\frac{d}{dt} (a_{4}t^{2} + b_{4}t + c_{4})|_{t=22.5} = \frac{d}{dt} (a_{5}t^{2} + b_{5}t + c_{5})|_{t=22.5}

(2a_{4}t + b_{4})|_{t=22.5} = (2a_{5}t + b_{5})|_{t=22.5}
2a_{4}(22.5) + b_{4} = 2a_{5}(22.5) + b_{5}
45a_{4} + b_{4} − 45a_{5} − b_{5} = 0

We have 14 of the 15 needed equations. The t = 0 and t = 30 points are exterior and only have 1 spline. We create the 15th

equation by assuming a linear equation at t = 0:

a_{1} = 0

which implies that

v(t) = a_{1}t^{2} + b_{1} + c_{1} = 0

becomes

v(t) = (0)t^{2} + b_{1} + c_{1} = 0
v(t) = b_{1} + c_{1} = 0

These can be created into a matrix equation of 15 equations and 15 unknowns.

[ \begin{matrix}0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 100 & 10 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 100 & 10 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 225 & 15 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 225 & 15 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 400 & 20 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 400 & 20 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 506.25 & 22.5 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 506.25 & 22.5 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 900 & 30 & 1 \\ 20 & 1 & 0 & −20 & −1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 30 & 1 & 0 & −30 & −1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 40 & 1 & 0 & −40 & −1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 45 & 1 & 0 & −45 & −1 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix} ] [ \begin{matrix}a_{1} \\ b_{1} \\ c_{1} \\ a_{2} \\ b_{2} \\ c_{2} \\ a_{3} \\ b_{3} \\ c_{3} \\ a_{4} \\ b_{4} \\ c_{4} \\ a_{5} \\ b_{5} \\ c_{5} \end{matrix} ] = [ \begin{matrix}0 \\ 227.04 \\ 227.04 \\ 362.78 \\ 362.78 \\ 517.35 \\ 517.35 \\ 602.97 \\ 602.97 \\ 901.67 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix} ]

The coefficient matrix is:

— * —

--S 1 of 16
incoef:=[_
[0  ,  0, 1,   0,  0, 0,   0,  0, 0,   0,     0,   0,   0,     0,   0],_
[100, 10, 1,   0,  0, 0,   0,  0, 0,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0, 100, 10, 1,   0,  0, 0,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0, 225, 15, 1,   0,  0, 0,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0, 225, 15, 1,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0, 400, 20, 1,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0,   0,  0, 0, 400,    20,   1,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0,   0,  0, 0, 506.25, 22.5, 1,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0,   0,  0, 0,   0,     0,   0, 506.25, 22.5, 1],_
[  0,  0, 0,   0,  0, 0,   0,  0, 0,   0,     0,   0, 900,    30,   1],_
[ 20,  1, 0, -20, -1, 0,   0,  0, 0,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0,  30,  1, 0, -30, -1, 0,   0,     0,   0,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0,  40,  1, 0, -40,    -1,   0,   0,     0,   0],_
[  0,  0, 0,   0,  0, 0,   0,  0, 0,  45,     1,   0, -45,    -1,   0],_
[  1,  0, 0,   0,  0, 0,   0,  0, 0,   0,     0,   0,   0,     0,   0]]
--R
--R
--R   (1)
--R   [[0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [100.0,10.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,100.0,10.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,225.0,15.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,225.0,15.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,400.0,20.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,400.0,20.0,1.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,506.25,22.5,1.0,0.0,0.0,0.0],

--R    [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,506.25,22.5,1.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,900.0,30.0,1.0],
--R    [20.0,1.0,0.0,- 20.0,- 1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,30.0,1.0,0.0,- 30.0,- 1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,40.0,1.0,0.0,- 40.0,- 1.0,0.0,0.0,0.0,0.0],
--R    [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,45.0,1.0,0.0,- 45.0,- 1.0,0.0],
--R    [1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]]
--R                                                      Type: List(List(Float))
--E 1

———-

and the values are:

— * —

--S 2 of 16
vals:Vector(Float):=_
[0,227.04,227.04,362.78,362.78,517.35,517.35,602.97,602.97,901.67,0,0,0,0,0]
--R
--R
--R   (2)
--R   [0.0, 227.04, 227.04, 362.78, 362.78, 517.35, 517.35, 602.97, 602.97,
--R    901.67, 0.0, 0.0, 0.0, 0.0, 0.0]
--R                                                          Type: Vector(Float)
--E 2
--S 3 of 16
digits(5)
--R
--R
--R   (3)  20
--R                                                        Type: PositiveInteger
--E 3
--S 4 of 16
outcoef:=solve(incoef,vals)
--R
--R
--R   (4)
--R   [
--R     particular =
--R       [0.0, 22.704, 0.0, 0.88882, 4.926, 88.888, - 0.1356, 35.66, - 141.6,
--R        1.6045, - 33.94, 554.41, 0.20905, 28.851, - 152.01]
--R     ,
--R    basis= [[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]]]
--RType: Record(particular: Union(Vector(Float),"failed"),basis: List(Vector(Float)))
--E 4
--S 5 of 16
D:=outcoef.particular

--R
--R
--R   (5)
--R   [0.0, 22.704, 0.0, 0.88882, 4.926, 88.888, - 0.1356, 35.66, - 141.6, 1.6045,
--R    - 33.94, 554.41, 0.20905, 28.851, - 152.01]
--R                                               Type: Union(Vector(Float),...)
--E 5
--S 6 of 16
s1:=D.1*t^2+D.2*t+D.3
--R
--R
--R   (6)  22.704 t
--R                                                      Type: Polynomial(Float)
--E 6
--S 7 of 16
s2:=D.4*t^2+D.5*t+D.6
--R
--R
--R                 2
--R   (7)  0.88882 t  + 4.926 t + 88.888
--R                                                      Type: Polynomial(Float)
--E 7
--S 8 of 16
s3:=D.7*t^2+D.8*t+D.9
--R
--R
--R                  2
--R   (8)  - 0.1356 t  + 35.66 t - 141.6
--R                                                      Type: Polynomial(Float)
--E 8
--S 9 of 16
s4:=D.10*t^2+D.11*t+D.12
--R
--R
--R                2
--R   (9)  1.6045 t  - 33.94 t + 554.41
--R                                                      Type: Polynomial(Float)
--E 9
--S 10 of 16
s5:=D.13*t^2+D.14*t+D.15
--R
--R
--R                  2
--R   (10)  0.20905 t  + 28.851 t - 152.01
--R                                                      Type: Polynomial(Float)

--E 10

———-

What is the velocity at t = 16?

— * —

--S 11 of 16
s3(16)
--R
--R
--R   (11)  394.26
--R                                                                  Type: Float
--E 11

———-

What is the acceleration at t = 16 — * —

--S 12 of 16
s3d:=differentiate(s3,t)
--R
--R
--R   (12)  - 0.27112 t + 35.66
--R                                                      Type: Polynomial(Float)
--E 12
--S 13 of 16
s3d(16)
--R
--R
--R   (13)  31.323
--R                                                                  Type: Float
--E 13

———-

Find the distance covered by the rocket from t = 11 to t = 16 Here we would want to do

∫{16 \atop 11}

but this uses two spline approximations so we have to break the integration at the interior point and compute:

∫_{11}^{16} = ∫_{11}^{15} s2 + ∫_{15}^{16} s3

— * —

--S 14 of 16
s2i:=integrate(s2,t)
--R
--R
--R                  3          2
--R   (14)  0.29628 t  + 2.463 t  + 88.888 t
--R                                                      Type: Polynomial(Float)
--E 14
--S 15 of 16
s3i:=integrate(s3,t)
--R
--R
--R                     3          2
--R   (15)  - 0.045187 t  + 17.83 t  - 141.6 t
--R                                                      Type: Polynomial(Float)
--E 15

———-

Now we can evaluate the difference of both of these integrals at their endpoints and sum the results to get the answer.

— * —

--S 16 of 16
(s2i(15)-s2i(11)) + (s3i(16)-s3i(15))
--R
--R
--R   (16)  1595.9
--R                                                                  Type: Float
--E 16
)spool
)lisp (bye)

———-

References

[1] http://www.youtube.com/watch?v=ES5zkgyMeAM