next up previous contents index
Next: 4.7 Interval arithmetic Up: 4. Nonlinear Polynomial Solvers Previous: 4.5 Auxiliary variable method   Contents   Index

4.6 Robustness issues

Current state-of-the-art CAD systems used to create and interrogate curved objects are based on geometric solid modeling methods that typically operate in floating point arithmetic (FPA). Arithmetic operations, especially division, in FPA lead to significant numerical errors. Division operation can be avoided by four-dimensional homogeneous processing proposed by Yamaguchi [285,456]. Furthermore CAD systems will frequently fail as a result of the limited precision that is inherent to the internal representation of floating point numbers [167,168]. One has to keep in mind that any sequence of operations on a digital computer is essentially equivalent to a finite sequence of manipulations on a discrete grid of points. For example, a floating point (FP) number in general form is given by [126]
$\displaystyle (\pm).b_1b_2\cdots b_p\cdot2^{E},$     (4.34)

where $ b_1\cdots b_{p}$ is the mantissa made up of binary digits 0 or 1, $ b_i=0$ or 1, with $ b_1 \neq 0$ , and $ p$ is the number of significant digits, and $ E$ is an integer exponent. If $ p=2$ and $ -2\leq E\leq3$ then a list of positive numbers in this system is

\begin{displaymath} \begin{array}{ll} .10*2^{-2}= \frac{1}{8},\quad\quad & .11*2... ...\;, \\ .10*2^{3}= 4,\quad\quad &.11*2^{3}= 6\;,\\ \end{array}\end{displaymath}


and are plotted in Fig. 4.7. Obviously, the resulting set of FP numbers is a finite subset of the rational numbers $ \frac{m}{n}$ (where $ m$ , $ n$ are integers) in the interval $ [\frac{1}{8},6]$ and they are distributed non-uniformly in this interval.
Figure 4.7: Nonnegative floating point numbers on the interval [0,6] (adapted from [126])
\begin{figure}\centerline{ \psfig{figure=fig/float.eps,height=1.0in} }\end{figure}

Nonlinear polynomial solvers operating in rational arithmetic (RA), where the arithmetic is done with rational numbers without approximation [204], are robust, but are generally memory intensive and time consuming due to the growth of the number of digits needed to represent rational numbers that result from arithmetic operations on other rational numbers. On the other hand, nonlinear solvers operating in floating point arithmetic are faster, but generally not robust. Interval methods, which are described in Sect. 4.7, effectively solve these two problems, namely, nonlinear polynomial solvers operating in interval arithmetic (IA) are inexpensive compared to rational arithmetic, and they are robust in eliminating regions not containing roots.

Figure 4.8: Curves $ y=x^{4}$ and $ y=0$ contact tangentially at the origin (adapted from [179])
\begin{figure}\vspace{0.0in} \centerline{ \psfig{figure=fig/t1_t4.ps,height=2.4in,width=3.4in} } \vspace*{-0.8in}\end{figure}

Example 4.6.1. Suppose we have a degree four planar Bézier curve whose control points are given by

$\displaystyle (-0.5,0.0625),\;(-0.25, -0.0625),\;(0,0.0625),\;(0.25, -0.0625),\;(0.5, 0.0625)\;,$      

as shown in Fig. 4.8. This Bézier curve is equivalent to the explicit curve $ y=x^4$ ( $ -0.5\leq x \leq 0.5$ ). Apparently the curve intersects with $ x$ -axis tangentially at $ (x,y)=(0,0)$ . However, if the curve has been translated by +1 in the $ y$ direction and translated back to the original position by moving by $ -\frac{1}{3}$ three times during a geometric processing session, the curve will generally not be the same as the original curve if floating point arithmetic is used for the computation. For illustration, let us assume a decimal computer with a four-digit mantissa, and the computer rounds off intelligently rather than truncating. Then the rational number $ -\frac{1}{3}$ will be stored in the decimal computer as $ -0.3333\times10^0$ and after the processing the new control points will be
$\displaystyle (-0.5,0.0631),\; (-0.25, -0.0624),\; (0,0.0631),\;(0.25, -0.0624),\; (0.5, 0.0631)\;.$      

If we evaluate the curve at parameter value $ t=0.5$ , we obtain (0, 0.00035) instead of (0,0). Therefore there exists a numerical gap which could later lead to inconsistency between topological structures and geometric equations. For example, if these new control points are used for computing intersections with the $ x$ -axis, the algorithm will return no solutions when the tolerance for the function value is smaller than 0.00035.

The above problem illustrates the case when the error is created during the formulation of the governing equations by various algebraic transformations.

Example 4.6.2. This example finds the roots of a cubic polynomial equation $ (x-0.1)(x-0.6)(x-0.7)=0$ by the PP method over an interval $ 0 \leq x \leq 1$ operating in FPA. Conversion to the Bernstein form was performed in exact arithmetic. This particular example was run at a tolerance of $ 10^{-4}$ and the binary subdivision was conducted when the box size did not reduce more than $ 5\%$ from the previous step. The algorithm output is listed in Table 4.1. At iteration 9, PP algorithm loses the root 0.7 due to floating point rounding.


Table 4.1: $ (x-0.1)(x-0.6)(x-0.7)=0$ solved by PP method operating in FPA (adapted from [255])

Iter		 Bounding Box (FPA) Message

1		 [0,1]  

2		 [0.0763636363636364, 0.856]  

3		 [0.098187732239346, 0.770083868323999]  

4		 [0.0999880766853688, 0.72387404781026] Binary Subdivision

5		 [0.402239977003124, 0.704479954527487]  

6		 [0.550441290533288, 0.700214508664293]  

7		 [0.591018492648952, 0.700000534482207]  

8		 [0.599458794784619, 0.700000000003332] Binary Subdivision

9		 [0.649998841568898, 0.699999999999999] No Root in Box

10		 [0.599997683137796, 0.649998841568898] Root Found in Box

11		 [0.099999999478761, 0.402239977003124] Root Found in Box

This example illustrates another serious problem which arises from the usage of FPA in shape interrogation. To remedy such problems interval arithmetic research in geometric modeling has become quite active as we will see in subsequent Sects. 4.7 to 4.9.

next up previous contents index
Next: 4.7 Interval arithmetic Up: 4. Nonlinear Polynomial Solvers Previous: 4.5 Auxiliary variable method   Contents   Index
December 2009