<< . .

. 52
( : 95)



. . >>

4




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 ,

<< . .

. 52
( : 95)



. . >>