2

0

’2

’4

’2 0 2 4 6

FIGURE 8.13. Parametric splines for a spiral-like node distribution. The spline

of cumulative length is drawn in solid line

dimensional case). Composite parametric splines can be generated as well

by enforcing suitable continuity conditions (see [SL89]).

Program 68 - par spline : Parametric splines

function [xi,yi] = par spline (x, y)

t (1) = 0;

for i = 1:length (x)-1

t (i+1) = t (i) + sqrt ( (x(i+1)-x(i))ˆ2 + (y(i+1)-y(i))ˆ2 );

end

z = [t(1):(t(length(t))-t(1))/100:t(length(t))];

xi = spline (t,x,z);

yi = spline (t,y,z);

8.7.1 B´zier Curves and Parametric B-splines

e

The B´zier curves and parametric B-splines are widely employed in graph-

e

ical applications, where the nodes™ locations might be a¬ected by some

uncertainty.

Let P0 , P1 , . . . , Pn be n + 1 points ordered in the plane. The oriented

polygon formed by them is called the characteristic polygon or B´zier poly-

e

gon. Let us introduce the Bernstein polynomials over the interval [0, 1]

de¬ned as

n!

n

tk (1 ’ t)n’k = tk (1 ’ t)n’k ,

bn,k (t) =

k!(n ’ k)!

k

360 8. Polynomial Interpolation

for n = 0, 1, . . . and k = 0, . . . , n. They can be obtained by the following

recursive formula

bn,0 (t) = (1 ’ t)n

bn,k (t) = (1 ’ t)bn’1,k (t) + tbn’1,k’1 (t), k = 1, . . . , n, t ∈ [0, 1].

It is easily seen that bn,k ∈ Pn , for k = 0, . . . , n. Also, {bn,k , k = 0, . . . , n}

provides a basis for Pn . The B´zier curve is de¬ned as follows

e

n

0 ¤ t ¤ 1.

Bn (P0 , P1 , . . . , Pn , t) = Pk bn,k (t), (8.58)

k=0

This expression can be regarded as a weighted average of the points Pk ,

with weights bn,k (t).

The B´zier curves can also be obtained by a pure geometric approach

e

starting from the characteristic polygon. Indeed, for any ¬xed t ∈ [0, 1],

we de¬ne Pi,1 (t) = (1 ’ t)Pi + tPi+1 for i = 0, . . . , n ’ 1 and, for t

¬xed, the piecewise line that joins the new nodes Pi,1 (t) forms a polygon

of n ’ 1 edges. We can now repeat the procedure by generating the new

vertices Pi,2 (t) (i = 0, . . . , n ’ 2), and terminating as soon as the polygon

comprises only the vertices P0,n’1 (t) and P1,n’1 (t). It can be shown that

P0,n (t) = (1 ’ t)P0,n’1 (t) + tP1,n’1 (t) = Bn (P0 , P1 , . . . , Pn , t),

that is, P0,n (t) is equal to the value of the B´zier curve Bn at the points

e

corresponding to the ¬xed value of t. Repeating the process for several val-

ues of the parameter t yields the construction of the curve in the considered

region of the plane.

8

6

4

2

0

’2

’4

2 4 6 8 10 12 14 16

FIGURE 8.14. Computation of the value of B3 relative to the points (0,0), (4,7),

(14,7), (17,0) for t = 0.5, using the graphical method described in the text

Notice that, for a given node con¬guration, several curves can be con-

structed according to the ordering of points Pi . Moreover, the B´zier curve

e

8.7 Splines in Parametric Form 361

Bn (P0 , P1 , . . . , Pn , t) coincides with Bn (Pn , Pn’1 , . . . , P0 , t), apart from

the orientation.

Program 69 computes bn,k at the point x for x ∈ [0, 1].

Program 69 - bernstein : Bernstein polynomials

function [bnk]=bernstein (n,k,x)

if k == 0,

C = 1;

else,

C = prod ([1:n])/( prod([1:k])*prod([1:n-k]));

end

bnk = C * xˆk * (1-x)ˆ(n-k);

Program 70 plots the B´zier curve relative to the set of points (x, y).

e

Program 70 - bezier : B´zier curves

e

function [bezx,bezy] = bezier (x, y, n)

i = 0; k = 0;

for t = 0:0.01:1,

i = i + 1; bnk = bernstein (n,k,t); ber(i) = bnk;

end

bezx = ber * x (1); bezy = ber * y (1);

for k = 1:n

i = 0;

for t = 0:0.01:1

i = i + 1; bnk = bernstein (n,k,t); ber(i) = bnk;

end

bezx = bezx + ber * x (k+1); bezy = bezy + ber * y (k+1);

end

plot(bezx,bezy)

In practice, the B´zier curves are rarely used since they do not provide a

e

su¬ciently accurate approximation to the characteristic polygon. For this

reason, in the 70™s the parametric B-splines were introduced, and they are

used in (8.58) instead of the Bernstein polynomials. Parametric B-splines

are widely employed in packages for computer graphics since they enjoy

the following properties:

1. perturbing a single vertex of the characteristic polygon yields a local

perturbation of the curve only around the vertex itself;

2. the parametric B-spline better approximates the control polygon than

the corresponding B´zier curve does, and it is always contained within

e

the convex hull of the polygon.

362 8. Polynomial Interpolation

In Figure 8.15 a comparison is made between B´zier curves and para-

e

metric B-splines for the approximation of a given characteristic polygon.

FIGURE 8.15. Comparison of a B´zier curve (left) and a parametric B-spline

e

(right). The vertices of the characteristic polygon are denoted by —

We conclude this section by noticing that parametric cubic B-splines

allow for obtaining locally straight lines by aligning four consecutive ver-

tices (see Figure 8.16) and that a parametric B-spline can be constrained

at a speci¬c point of the characteristic polygon by simply making three

consecutive points of the polygon coincide with the desired point.

FIGURE 8.16. Some parametric B-splines as functions of the number and posi-

tions of the vertices of the characteristic polygon. Notice in the last ¬gure (right)

the localization e¬ects due to moving a single vertex

8.8 Applications

In this section we consider two problems arising from the solution of fourth-

order di¬erential equations and from the reconstruction of images in axial

tomographies.

8.8 Applications 363

8.8.1 Finite Element Analysis of a Clamped Beam

Let us employ piecewise Hermite polynomials (see Section 8.4) for the nu-

merical approximation of the transversal bending of a clamped beam. This

problem was already considered in Section 4.7.2 where centered ¬nite dif-

ferences were used.

The mathematical model is the fourth-order boundary value problem

(4.74), here presented in the following general formulation

(±(x)u (x)) = f (x), 0 < x < L

(8.59)

u(0) = u(L) = 0, u (0) = u (L) = 0.

In the particular case of (4.74) we have ± = EJ and f = P ; we assume

henceforth that ± is a positive and bounded function over (0, L) and that

f ∈ L2 (0, L).

We multiply (8.59) by a su¬ciently smooth arbitrary function v, then,

we integrate by parts twice, to obtain

L L

L L

±u v dx ’ [±u v]0 + [±u v ]0 = f vdx.

0 0

Problem (8.59) is then replaced by the following problem in integral form

L L

¬nd u ∈ V such that ∀v ∈ V, (8.60)

±u v dx = f vdx,

0 0

where

V = v : v (k) ∈ L2 (0, L), k = 0, 1, 2, v (k) (0) = v (k) (L) = 0, k = 0, 1 .

Problem (8.60) admits a unique solution, which represents the deformed

con¬guration that minimizes the total potential energy of the beam over

the space V (see, for instance, [Red86], p. 156)

L

1

±(u )2 ’ f u dx.

J(u) =

2

0

In view of the numerical solution of problem (8.60), we introduce a partition

Th of [0, L] into K subintervals Tk = [xk’1 , xk ], (k = 1, . . . , K) of uniform

length h = L/K, with xk = kh, and the ¬nite dimensional space

vh ∈ C 1 ([0, L]), vh |T ∈ P3 (T )

Vh =

(8.61)

(k) (k)

∀T ∈ Th , vh (0) = vh (L) = 0, k = 0, 1 .

364 8. Polynomial Interpolation

Let us equip Vh with a basis. For this purpose, we associate with each

internal node xi (i = 1, . . . , K ’ 1) a support σi = Ti ∪ Ti+1 and two

functions •i , ψi de¬ned as follows: for any k, •i |Tk ∈ P3 (Tk ), ψi |Tk ∈ P3 (Tk )

and for any j = 0, . . . , K,

±

•i (xj ) = δij , •i (xj ) = 0,

(8.62)

ψi (xj ) = 0, ψi (xj ) = δij .

Notice that the above functions belong to Vh and de¬ne a basis

Bh = {•i , ψi , i = 1, . . . , K ’ 1}. (8.63)

ˆ

These basis functions can be brought back to the reference interval T =

ˆ

[0, 1] for 0 ¤ x ¤ 1, by the a¬ne maps x = hˆ + xk’1 between T and Tk ,

ˆ x

for k = 1, . . . , K.

(0)

ˆ

Therefore, let us introduce on the interval T the basis functions •0ˆ

(1) (0) (1)

and •0 , associated with the node x = 0, and •1 and •1 , associated

ˆ ˆ ˆ ˆ

with node x = 1. Each of these is of the form • = a0 + a1 x + a2 x2 +

ˆ ˆ ˆ ˆ

3

a3 x ; in particular, the functions with superscript “0” must satisfy the

ˆ

¬rst two conditions in (8.62), while those with superscript “1” must ful¬ll

the remaining two conditions. Solving the (4—4) associated system, we get

(0) (1)

•0 (ˆ) = 1 ’ 3ˆ2 + 2ˆ3 , •0 (ˆ) = x ’ 2ˆ2 + x3 ,

ˆx x x ˆx ˆ x ˆ

(8.64)

(0) (1)

= 3ˆ ’ 2ˆ , = ’ˆ + x .

2 3 2 3

•1 (ˆ)

ˆx x x •1 (ˆ)

ˆx x ˆ

The graphs of the functions (8.64) are drawn in Figure 8.17 (left), where

(0) (0) (1) (1)

(0), (1), (2) and (3) denote •0 , •1 , •0 and •1 , respectively.

ˆ ˆ ˆ ˆ

The function uh ∈ Vh can be written as

K’1 K’1

(1)

uh (x) = ui •i (x) + ui ψi (x). (8.65)

i=1 i=1

The coe¬cients and the degrees of freedom of uh have the following mean-

(1)

ing: ui = uh (xi ), ui (xi ) = uh (xi ) for i = 1, . . . , K ’ 1. Notice that (8.65)

is a special instance of (8.32), having set mi = 1.

The discretization of problem (8.60) reads

L L

¬nd uh ∈ Vh such that ∀vh ∈ Bh . (8.66)

±uh vh dx = f vh dx,

0 0

This is called the Galerkin ¬nite element approximation of the di¬erential

problem (8.59). We refer to Chapter 12, Sections 12.4 and 12.4.5, for a

more comprehensive discussion and analysis of the method.

8.8 Applications 365

5

1 10

(0) (1)

0.8 0

10

0.6

’5

10

0.4

’10

10

0.2

(2)

Prec. No Prec.

’15

10

0

(3)

’20

’0.2 10

0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 120 140 160

FIGURE 8.17. Canonical Hermite basis on the reference interval 0 ¤ x ¤ 1 (left);

ˆ

convergence histories for the conjugate gradient method in the solution of system

(8.69) (right). On the x-axis the number of iterations k is shown, while the y-axis

represents the quantity r(k) 2 / b1 2 , where r is the residual of system (8.69)

Using the representation (8.65) we end up with the following system in

(1) (1) (1)

the 2K ’ 2 unknowns u1 , u2 , . . . , uK’1 , u1 , u2 , . . . uK’1

±

±

K’1 L L L

(1)

uj ±•j •i dx + uj ±ψj •i dx = f •i dx,

j=1

± 0L 0L

0 (8.67)

L

K’1

(1)

uj ±•j ψi dx + uj ±ψj ψi dx = f ψi dx,

j=1 0 0 0

for i = 1, . . . , K ’1. Assuming, for the sake of simplicity, that the beam has

unit length L, that ± and f are two constants and computing the integrals

in (8.67), the ¬nal system reads in matrix form

Au + Bp = b1

(8.68)

BT u + Cp = 0,

(1)

where the vectors u, p ∈ RK’1 contain the nodal unknowns ui and ui ,