Fundamentals of Linear Algebra and Optimization Project 1

The purpose of this project is to solve a curve interpolation problem using cubic splines.
This type of problem arises frequently in computer graphics and in robotics (path planning).
Recall from Homework 1 that polynomial curves of degree ? m can be de?ned in terms of control points and the B?ezier polynomials. Cubic B?ezier curves are often used because they are cheap to implement and give more exibility than quadratic B?ezier curves.
A cubic B?ezier curve C(t) (in R2 or R3) is speci?ed by a list of four control points
(b0; b2; b2; b3) and is given parametrically by the equation
C(t) = (1 􀀀 t)3 b0 + 3(1 􀀀 t)2t b1 + 3(1 􀀀 t)t2 b2 + t3 b3:
Clearly, C(0) = b0, C(1) = b3, and for t 2 [0; 1], the point C(t) belongs to the convex hull of the control points b0; b1; b2; b3.
Interpolation problems require ?nding curves passing through some given data points and possibly satisfying some extra constraints.
It is known that Lagrange interpolation is not very satisfactory when N is greater than 5, since Lagrange interpolants tend to oscillate in an undesirable manner. Thus, we turn to B?ezier spline curves.
A B?ezier spline curve F is a curve which is made up of curve segments which are B?ezier curves, say C1; : : : ;Cm (m ? 2). We will assume that F de?ned on [0;m], so that for
i = 1; : : : ;m,
F(t) = Ci(t 􀀀 i + 1); i 􀀀 1 ? t ? i:
Typically, some smoothness is required between any two junction points, that is, between any two points Ci(1) and Ci+1(0), for i = 1; : : : ;m 􀀀 1. We require that Ci(1) = Ci+1(0) (C0-continuity), and typically that the derivatives of Ci at 1 and of Ci+1 at 0 agree up to second order derivatives. This is called C2-continuity, and it ensures that the tangents agree as well as the curvatures.
There are a number of interpolation problems, and we consider one of the most common problems which can be stated as follows:
1
x0 = d−1
x1
x2
x3
x4
x5
x6
x7 = d8
d0
d1
d2
d3
d4
d5
d6
d7
Figure 1: A C2 cubic interpolation spline curve passing through the points x0; x1, x2; x3,
x4; x5, x6; x7
Problem: Given N +1 data points x0; : : : ; xN, ?nd a C2 cubic spline curve F, such that F(i) = xi, for all i, 0 ? i ? N (N ? 2).
A way to solve this problem is to ?nd N +3 auxiliary points d􀀀1; : : : ; dN+1 called de Boor control points from which N B?ezier curves can be found. Actually,
d􀀀1 = x0 and dN+1 = xN
so we only need to ?nd N + 1 points d0; : : : ; dN.
It turns out that the C2-continuity constraints on the N B?ezier curves yield only N 􀀀 1
equations, so d0 and dN can be chosen arbitrarily. In practice, d0 and dN are chosen according to various \end conditions," such as prescribed velocities at x0 and xN. For the time being, we will assume that d0 and dN are given.
Figure 1 illustrates an interpolation problem involving N + 1 = 7 + 1 = 8 data points.
The control points d0 and d7 were chosen arbitrarily.
It can be shown that d1; : : : ; dN􀀀1 are given by the linear system
2
0
BBBBB@
7
2 1
1 4 1 0
. . . . . . . . .
0 1 4 1
1 7
2
1
CCCCCA
0
BBBBB@
d1
d2
...
dN􀀀2
dN􀀀1
1
CCCCCA
=
0
BBBBB@
6x1 􀀀 3
2d0
6x2
...
6xN􀀀2
6xN􀀀1 􀀀 3
2dN
1
CCCCCA
:
The derivation of the above system assumes that N ? 4. If N = 3, this system reduces
to
7
2 1
1 7
2
!
d1
d2
!
=

6x1 􀀀 3
2d0
6x2 􀀀 3
2d3
!
When N = 2, it can be shown that d1 is given by
d1 = 2x1 􀀀
1
2
d0 􀀀
1
2
d2:
Observe that when N ? 3, the above matrix is strictly diagonally dominant, so it is
invertible. Actually, the above system needs to be solved for the x-coordinates and for the y-coordinates of the dis (and also for the z-coordinates, if the points are in R3). Once the above system is solved, the B?ezier cubics C1; : : : ;CN are determined as follows (we assume N ? 2): For 2 ? i ? N 􀀀 1, the control points (bi
0; bi
1; bi
2; bi
3) of Ci are given by
bi
0 = xi􀀀1
bi
1 =
2
3
di􀀀1 +
1
3
di
bi
2 =
1
3
di􀀀1 +
2
3
di
bi
3 = xi:
The control points (b10
; b11
; b12
; b13
) of C1 are given by
b10
= x0
b11= d0
b12
=
1
2
d0 +
1
2
d1
b13
= x1;
and the control points (bN0; bN1
; bN2
; bN3
) of CN are given by
bN0
= xN􀀀1
bN1
=
1
2
dN􀀀1 +
1
2
dN
bN2
= dN
bN3
= xN:
3
Prove that the tangent vectors m0 at x0 and mN at xN are given by
m0 = 3(d0 􀀀 x0)
mN = 3(xN 􀀀 dN):
End Conditions.
One method to determine the points d0 and dN is the natural end condition, which
consists in setting the second derivatives at x0 and at xN to be zero, that is,
C00 1 (0) = 0; C00 N(1) = 0:
(1) Prove that the second derivative at b0 of a B?ezier cubic speci?ed by the control points (b0; b1; b2; b3) is
6(b0 􀀀 2b1 + b2); and the second derivative at b3 is
6(b1 􀀀 2b2 + b3):
Prove that our system becomes
0
BBBBB@
4 1
1 4 1 0
. . . . . . . . .
0 1 4 1
1 4
1
CCCCCA
0
BBBBB@
d1
d2
...
dN􀀀2
dN􀀀1
1
CCCCCA
=
0
BBBBB@
6x1 􀀀 x0
6x2
...
6xN􀀀2
6xN􀀀1 􀀀 xN
1
CCCCCA
;
where d0; dN are given by
d0 =
2
3
x0 +
1
3
d1
dN =
1
3
dN􀀀1 +
2
3
xN:
Note that d0 is on the line segment (x0; d1) (1=3 of the way from x0) and dN is on the line segment (dN􀀀1; xN) (1=3 of the way from xN).
In the above derivation we assumed that N ? 4. If N = 3, show that this system reduces
to ?
4 1
1 4
??
d1
d2
?
=
?
6x1 􀀀 x0
6x2 􀀀 x3
?
:
Another method to determine the points d0 and dN is the quadratic end condition, which consists in requiring that the second derivatives at x0 and at x1 agree, and similarly for the second derivatives at xN􀀀1 and xN; this means that
C00 1 (0) = C00 1 (1) and C00 N(0) = C00 N(1):
4
(2) Prove that our system becomes
0
BBBBB@
5 1
1 4 1 0
. . . . . . . . .
0 1 4 1
1 5
1
CCCCCA
0
BBBBB@
d1
d2
...
dN􀀀2
dN􀀀1
1
CCCCCA
=
0
BBBBB@
7x1 􀀀 x0
6x2
...
6xN􀀀2
7xN􀀀1 􀀀 xN
1
CCCCCA
;
where d0; dN are given by
d0 = d1 +
2
3
x0 􀀀
2
3
x1
dN = dN􀀀1 +
2
3
xN 􀀀
2
3
xN􀀀1:
Geometrically, the vector d0􀀀d1 is equal to 2
3 (x0􀀀x1) and similarly the vector dN 􀀀dN􀀀1
is equal to 2 3 (xN 􀀀 xN􀀀1); in particular, the line segments (d0; d1) and (x0; x1) are parallel, and so are (dN; dN􀀀1) and (xN; xN􀀀1).
In the above derivation we assumed that N ? 4. If N = 3, show that this system reduces
to ?
5 1
1 5
??
d1
d2
?
=
?
7x1 􀀀 x0
7x2 􀀀 x3
?
:
Yet another method known as Bessel end condition is to require the ?rst derivative C01 (0) at x0 to be equal to the ?rst derivative of the unique parabola passing through x0; x1; x2, for t = 0; 1; 2, and to require the ?rst derivative C0N (1) at xN to be equal to the ?rst derivative of the unique parabola passing through xN􀀀2; xN􀀀1; xN, for t = 0; 1; 2.
A parabola passing through x0 and x2 as above is given by
C(t) =
(2 􀀀 t)2
4
x0 +
(2 􀀀 t)t
2
b1 +
t2
4
x2;
so to require that C(1) = x1 means that
x1 =
1
4
x0 +
1
2
b1 +
1
4
x2;
which yields
b1 = 􀀀
1
2
x0 + 2x1 􀀀
1
2
x2:
(3) Show that
m0 = C0(0) = b1 􀀀 x0 = 􀀀
3
2
x0 + 2x1 􀀀
1
2
x2;
5
and using the fact that d0 = x0 + 1
3m0, show that
d0 =
1
2
x0 +
1
2
?
x1 +
1
3
(x1 􀀀 x2)
?
:
Geometrically, d0 is the midpoint of the line segment from x0 to a point obtained by extrap-
olation from x1 and x2 (x1 + 1
3 (x1 􀀀 x2)).
Similarly, for the parabola interpolating xN􀀀2; xN􀀀1 and xN, we get
bN = 􀀀
1
2
xN􀀀2 + 2xN􀀀1 􀀀
1
2
xN;
and show that
mN = C0(1) = xN 􀀀 bN =
1
2
xN􀀀2 􀀀 2xN􀀀1 +
3
2
xN;
and using the fact that dN = xN 􀀀 1
3mN, show that
dN =
1
2
?
xN􀀀1 +
1
3
(xN􀀀1 􀀀 xN􀀀2)
?
+
1
2
xN:
Geometrically, dN is the midpoint of the line segment from xN to a point obtained by extrapolation from xN􀀀1 and xN􀀀2 (xN􀀀1+ 1
3 (xN􀀀1􀀀xN􀀀2)). Note that the above derivation
is correct for N ? 2. Finally, there is the not a knot end condition, which consists in forcing the ?rst two B?ezier segments C1 and C2 to belong to the same cubic curve, and similarly for the last two B?ezier segments CN􀀀1 and CN. This amounts to require that C000 1 (1) = C000 2 (0) and C000 N􀀀1(1) = C000 N(0).
(4) Prove that the third derivative at b0 and at b3 of a B?ezier cubic speci?ed by the
control points (b0; b1; b2; b3) is
6(􀀀b0 + 3b1 􀀀 3b2 + b3):
Prove that if N = 3, then
d0 =
7
18
x0 +
8
9
x1 +
7
18
x2 􀀀
2
3
d2
d1 = 􀀀
1
6
x0 +
4
3
x1 􀀀
1
6
x2
d2 = 􀀀
1
6
x1 +
4
3
x2 􀀀
1
6
x3
d3 =
7
18
x1 +
8
9
x2 +
7
18
x3 􀀀
2
3
d1
are already computed in terms of x0; : : : ; x3, and there is no need to solve any linear system.
6
Prove that if N = 4, then
d1 = 􀀀
1
6
x0 +
4
3
x1 􀀀
1
6
x2
d3 = 􀀀
1
6
x2 +
4
3
x3 􀀀
1
6
x4
d2 =
3
2
x2 􀀀
1
4
d1 􀀀
1
4
d3
d0 =
7
18
x0 +
8
9
x1 +
7
18
x2 􀀀
2
3
d2
d4 =
7
18
x2 +
8
9
x3 +
7
18
x4 􀀀
2
3
d2:
Prove that if N ? 5, our linear system becomes the (N 􀀀 3) ? (N 􀀀 3) system
0
BBBBB@
4 1 0 0
1 4 1 0
. . . . . . . . .
0 1 4 1
0 0 1 4
1
CCCCCA
0
BBBBB@
d2
d3
...
dN􀀀3
dN􀀀2
1
CCCCCA
=
0
BBBBB@
6x2 + 1
6x0 􀀀 4
3x1 + 1
6x2
6x3
...
6xN􀀀3
6xN􀀀2 + 1
6xN􀀀2 􀀀 4
3xN􀀀1 + 1
6xN
1
CCCCCA;
and d0; d1; dN􀀀1; dN are given by
d0 =
7
18
x0 +
8
9
x1 +
7
18
x2 􀀀
2
3
d2
d1 = 􀀀
1
6
x0 +
4
3
x1 􀀀
1
6
x2
dN􀀀1 = 􀀀
1
6
xN􀀀2 +
4
3
xN􀀀1 􀀀
1
6
xN
dN =
7
18
xN􀀀2 +
8
9
xN􀀀1 +
7
18
xN 􀀀
2
3
dN􀀀2:
If N = 5, this system reduces to

4 1
1 4
!
d2
d3
!
=

6x2 + 1
6x0 􀀀 4
3x1 + 1
6x2
6x3 + 1
6x3 􀀀 4
3x4 + 1
6x5
!
:
(5) Implement the Gaussian elimination method with partial pivoting as well as the
method for solving a triangular system by back-substitution.
Use your program to solve several instances of the interpolation problem. Verify that no pivoting is needed.
(6) Implement the LU-factorization method for tridiagonal matrices and test it on the same interpolation problems as in (5).

Do you notice any improvement over Gaussian elimination (running time, numerical precision)?
(7) After computing d1; : : : ; dN􀀀1, compute the control points for the B?ezier curves C1; : : :CN and write a program to plot these B?ezier segments (for t 2 [0; 1]) to visualize the interpolating spline. Experiment with the choice of end conditions.
TOTAL: 300 points.