The application of piecewise polynomials to problems of curve and surface approximation

RIJKSW A TERSTAAT COMMUNICATIONS

THE APPLICATION

OF

PIECEWISE

POLYNOMIALS

TO PROBLEMS

OF CURVE

AND

SURFACE

APPROXIMATION

Editor:

Dr. Kurt Kubik

Contributors:

Ir. EricR. Bosman, Dienst Informatieverwerking Rijkswaterstaat Prof. David Eckhart, Adviser of Rijkswaterstaat

Dr. Kurt Kubik, Dienst Informatieverwerking Rijkswaterstaat

Any correspondence should be addressed to

DIRECTIE WATERHUISHOUDING EN WATERBEWEGING THE HAGUE - NETHERLANDS

Contents

page

4 I'reface

5 Introduction

7 _{I. The interpolation of smootb curves}

23

_n.

Tbe computation of elastic beams witb large deftection

34 ill. Approximation of measured data by piecewise bicubic polynomial functions

52 IV. The numerical solution of variational problems, using piecewise polynomial functions and a least squares approximation to the boundary conditions

55 V. The use of piecewise polynomial functions for the approximation of tbe relationship between any two geodetic coordinate systems 63 VI. Fast and precise transformation of Decca coordinates to VTM

coordinates or to other coordinate systems

78 Vl], Methods for the transformation of points from sidelooking radar and infrared photograpbs to the map system

Preface

This publication gives the approximate solution to a variety of technical problems, the majority of which were encountered at the Meetkundige Dienst of Rijkswater-staat. The use of piecewise polynomial functions is a common factor to all the solutions. This publication consists of seven Papers, preceded by a general introduction. The first Paper deals with the approximation of a function of one variabie and the third with approximation of functions of two variables. The other Papers deal primarily with applications to technical problems.

The authors would like to thank their co-workers, who have contributed to the general development, for their assistance in the system analysis, program ming and data processing of the various applications of the theory. In particular the work of Mr. E. Clerici, Mr. J. W. J. Koelet, Mr. A.Kranendonk, Mr. J. Scheringa and Miss B. J.Vreugdenhil must be mentioned. Drs. D. Eckhart, professor at the ITC, con-tributed to this publication in his capacity as scientific adviser of Rijkswaterstaat, Dienst Informatieverwerking.

K. Kubik

Delft, the Netherlands May 1971

Introduetion

What are piecewise polynomials?

Itseems appropriate to begin this publication by defining the concept of piecewise polynomials and also to indicate the idea leading to this definition.

For many years, long, thin strips of wood or some other material have been used by draftsmen to fair in a smooth curve between specified points. These strips or mechanical splines are kept in place by attaching lead weights called "ducks " at points along the strip.

In order to represent sirnilar smooth curves mathematically, the concept of piece-wise polynomial functions was introduced. To this purpose, the individual curve sections of the mechanical spline between two adjacent ducks are replaced by different polynomials, which join smoothly at the locations of the ducks. The piecewise polynomial is then a smooth function consisting of different polynomials between each pair of adjacentjunction points (ducks).

Various conditions may now be imposed upon the shape ofthe piecewise polynomial function. Possible conditions are, that either the slope or the curvature of the curve is

in the average as small as possible. .

The idea of piecewise polynornial functions can also be extended to functions of two or more variables. For two variables the piecewise polynomial function, which may be visualized as a surface, consists of parts of polynomial surfaces, which join continuously with the adjacent surface sections.

Piecewise polynomial functions are found to have highly desirabIe characteristics as approximating, interpolating and curve fitting functions.

I.

The interpolation of smooth

curves

1. Introduction

This report documents a number of methods for the interpolation of smooth curves. The problem consists of computing a curve which passes through a given set of points (x., yJ where some of the derivatives are also given. A result is required which is precise to 7-9 places.This problem originated in road design work, but is also common to many other scientific and engineering tasks. The method to be described is a rather elaborate one. lts use is justified by the high level of precision which it affords.

To make the interpolating curve independent of the orientation of the coordinate system chosen for the computation, the curve is represented in its natural equation

H = H(s)

with H inclination angle of the curve and s length of are of the curve.

Piecewise polynomial functions of various degree (degrees 1, 2 and 3) are used to describe the relationship between Hand s. The mesh locations are either fixed or variabie. Variabie mesh locations are used to achieve a smooth function with a limited number of mesh locations. The detailed computation methods, along with numerical examples, are presented in the following sections.

2. Statement of tbeproblem

The task involves computation of a smooth curve y = y(x) in E2' satisfying the following constraints :

The curve passes through a number of points (x., Yj), i =1... N;

The angular inclination of the curve and/or the curvature and/or (he change of curvature have predetermined values at some of the points (x., yj) or within certain given intervals of the curve.

Using the natural representation of the curve, these constraints may be formally writ

or in more detail

f

Si+1

1

(Xi+1

-

x

J

-

cosH ds = 0 Si

Si+1

(Yi+1-Yi)-f sinHds =O;i = 1,... N-1 Si

(2a)

Si length of are of data points

(2b) for some values j

out of (1,2... N) or for all values s in

some given mesh (sa, Sa+1)

ofthe curve

(2c)

(2d)

The smoothness of the curve is ensured by imposing the following condition L

I =

f

0(QH' H2 + QK' K2 + QA .A2) ds -+min (3)

withQH,QKand QA properly chosen parameters,

L total length of are of the curve.

3. Solution of the problem

A piecewise polynornial function is chosen for the function H'(s), For this purpose we assume a partition of the interval (O,L) into meshes according to the mesh locations 0

=

So <Sl <S2 ... <sa

=

L.

A piecewise polynomial function equals in each mesh a polynomial of degree n, and has continuous m derivatives in the interval (0,L), m ~ n- 1.In this piecewise polynomial function H (s) the magnitudes of a number of parameters bare unknown. Those magnitudes are determined according to the constraints (I) together with con-dition (3) or by the equivalent formulation

J(b) = leb)

+

LLAp' Cp(b) -+ min

p

(4)

The total number ofunknowns b should alwaysbe larger than the total number of constraints imposed on the curve. The number of mesh points must be selected according to thisrule. The more unknowns there are in the problem, the greater will be the effect of the smoothness condition (3)on the computation.

Piecewise polynomials of degrees 1,2and 3 are used for the function H(s). The corresponding three methods of computation are described in the following three sections.

4. Piecewiselinear functions 4.1. Definition and parametrization

Let A:

°

= so < SI' .. <Sa= L be aset of mesh points. The function H(s) is a linear function of s within each of the meshes (S"-I' s,,),a =1... a, and is

con-tinuous over the interval (O,L).

As parameters of the curvewe choose:

The inclination angles H"at the mesh points, a =

°

.

..

a. The length of the are S",a =

°

.

..

a at the mesh points.

For the mesh (s",Sdi) the function H(s) may then be written as

(5)

When the choice of parameters is made in this way, the function Hts) will be continuous over the interval (O,L).

4.2. The unknowns andthe conditions

The mesh points do not in general coincide with the values s,at the data points. A typical interpolation curve into N data points is shown in Figure 1.The magnitudes of H,,(a

=

0, 1... a), s,,(a

=

1, ... (a- 1)) at the mesh locations are still un-known, as are the values of s,(i = 2 N) at the data points. The value SI at the first data point must be given to assure a solution.

These unknowns are determined by finding the minimum value of J,

J = I

+

L LAp' C, --+ min (4)

p

with constraints of type (2a),(2b) and (2c),while the parameters QA>QHin formula (3) are set equal to 0.The condition which specifies that the total number ofunknowns

Figure 1. The parameters of the piecewise polynomial of degree one

o

MESH POINT Je DATA POINT

a

KNOWN PARAMETER

a UNKNOWN PARAMETER

must be larger than the total number of constraints results in the following equality for the number a of mesh points in the curve

a

~

HN

+

q-1)

with q total number of constraints of type (2b - d).

4.3. The solution

4.3.1. The problem is solved by equating the fint derivates of J to the unknowns to

zero

8J

8b = 0; bT (6)

This results in a nonlinear system of equations which is solved iteratively using the

Newton process. The initial approximations for sa must be given; to select them

properly is not an easy task, particularly if the number of unknowns equals the

number of constraints. Initial approximations for Ha and Sjare computed using local

parabolic approximation. The initial approximation 0 is substituted for the

un-knowns LAp. The algorithm of the Newton process, which is used to improve the approximate values, is as follows

( 8b8b8JT )bp-dbp

+

( 8J )-Sb bp =0',

bp+1 = bp

+

db,

~ =

0,1,2, ...

This computational cycle is repeated until the largest residual- maximum value of the moduli of the elements in (

!~

)-

is smaller than a pre-determined tolerance. 4.3.2. Practical computations indicated that the results of an iteration step do not depend on the approximate values assigned to LAp. Consequently, the value

°

is used as an approximation to LAp within each iteration step, thus simplifying the computations of the matrix coefficients (

oJ )

to a considerable extent.

obob

T

4.3.3. After each iteration step the point sequence in the curve must be re-examined. A change in the sequence of the mesh locations (sa) relative to the data points (s.) will be tolerated. A change in the sequence oftwo or more mesh locations (sa+1 <sa) however, will not be tolerated. The computation wiJlnot continue should the sequence be altered in this way.

No further difficulties are encountered in the numerical computation ofthis problem.

For further details on the computation, reference is made to Section 5.3.

5. Piecewisequadratic functions 5.1. Definition and parametrization

LetLl :

°

= So <SI < ... <Sa = L be a set of mesh points. The function H(s) is

aquadratic function in s within each of the meshes (sa _ 1, SJ, a

=

1 ... a; its first derivative is continuous OVerthe interval (O,L).

The following magnitudes are selected as curve parameters: The inclination angles Ha at the mesh points, a = 0, ... a;

02H

The magnitudes Aa. = --2- for each mesh (Sa.-l' sa.), a = 1, ... a;

os

The length of are Sa.at the mesh points; a = 0, ... a. For the mesh (Sa.-l' sa) the function may then be written as

This function is continuous over the interval (O,L).lts first derivative is represented by

At the mesh location Sa this derivative is equal to the expressions

for the mesh (Sa+l' sj

oH

At the mesh location Sa'_- will be continuous, when the following condition is

os

enforced

----. Ha+1 sa+1 - Sa

(9)

5.2. The unknownsandthe conditions

The mesh locations do not in general coincide with the values s, at the data points. A typical interpolation curve into Ndata points is shown in Figure 2. The magnitudes Aa(a = 1 ... a), Ha(a = O ... a), sa(u = 1 ... (a- 1)) at the mesh locations and the values Sj(i = 2 ... N) at the data points are still unknown. The value S1at the first data point must be given to assure a solution.

The unknowns are determined from the extremum problem (4) with constraints of type (2a), (2b),(2c) and (2d) and the additional constraints (9).

The parameter QHin I, (3) is set equal to O.

The condition that the total number ofunknowns is larger than the total number of constraints results in the following relationship for the number of mesh points in the curve

a

~

HN

+

q-2)

with q total number of constraints of type (2b- d).

This relation holds for variabIe mesh locations. For the case where the mesh locations Sa are regarded as fixed (thus making the problem more stabIe), the total number of unknowns is decreased. Consequently a larger number of mesh points are required to compute the interpolative curve

Figure 2. The parameters of the piecewise polynomial ofdegreetwo

A

1 5

2

51

X

2

S1 H1

Y2

X1 SO

Y

1 Ho

o

MESH POINT x DATA POINT

a

KNOWN PARAMETER a UNKNOWN PARAMETER

5.3. The solution

5.3.1. The problem (4) can be reduced to the following non-linear equation system for the unknowns

for variable mesh locations for fixed mesh locations The Newton process will be employed for the so!ution of this equation system. If the specific va!ues are not known, initia! approximations for Samust be given. lt is recommended that these initial approximations be chosen ha!fway between the known data points. If locally more data (H], Kj, Aj) are given,a locally dense sub-division of the curve into meshes is required. Initia! approximations for Aa, Ha and s, are computed from the data points by loca! parabo!ic interpo!ation, and

°

is chosen asinitia! approximation to LAp'

5.3.2. In the computation of the coefficient matrix ( öJ ) (7)within the Newton öböbT

process, the most cumbersome operation is the computation of the integrals

c

,

=

f

12SkcosH(s) ds ;Sk

=

f

12SksinH(s) ds ;

II II

k

=

0,1,2,3,4, sa~ t1, t2 ~ Sa+1

The use of Simpsou's formula for integration proved to be too time consuming. Consequently an analytica! eva!uation of the integrals was sought.

At present they are evaluated as follows (b2 - ac)

[

J

21t

J}

S = t2

+

sin a /' S --;-(as

+

b) s = ti ( b2 - ac)

[

J

21t

J}

S = t2 -sin a ·C --;-(as

+

b) s

=

ti with ( 1t

)2n

00 (-lt

-f

z (1t 2) ~ 2 4n+1 C(z)

=

0cos

2

t dt

=

L..

i

2n)! (4n

+

1) z n;O ( 1t

)

2

n

+

l

00 (-lt

-f

z

(1t)

I

2 _ . _ 2 _ z4n+3 S(z) - 0 sm 2 t dt - (2n

+

1)!(4n

+

3) n;O

and a, b, c denoting the coefficients of the polynomial H(s) = as2

+

bs

+

c; s,~ s~ S,,+ 1

(10)

For the evaluation of the integrals (10) the series expansions for Fresnel Integrals are used. The other integrals Ck, Sk are then computed from Co, So by partial integra.. tion.

A number of supplementary integration routines are provided for those cases,

where either the above expressions are not sufficiently accurate or do not hold (for a close to 0).

5.3.3. Theconvergenee of the magnitudes S"is frequently of oscillatory nature, with

oscillation of considerable magnitude. In order to dampen this oscillation, diagonal terms are added to the equation matrix (7) in the columns of ds;Their magnitude are chosen such that the largest corrections ds, are in the order of 0.1 times the mesh

length. This measure does of course delay the average rate of convergency of the process; however, it prevents divergencein individual computations.

For the computation with fixed mesh locations sa' the diagonal terms 10600 are added into the equation matrix inthecolumns ofds;

5.3.4. Within each iteration step of the Newton process the linear equation system is solved by direct elirnination of the unknowns. Advantage isthereby taken of the small percentage of non-zero coefficients in the equation matrix. In practical com-putations 3-5iteration steps are necessary to obtain the solution ofthe problem with the required precision of7-9 places.

5.3.5. It proved bypractical computations, that the results ofoneiteration step do not depend on the approximate value for LAp. The value 0 is therefore used as approximation for LAp within every iteration step,thus simplifying the computations of the matrix coefficients to a considerable extent.

The computed curve may showawavelike course at either end. This behaviour is suppressed by enforcing the constraint (

°o~)

= 0 on the first and last mesh of the

a

curve. Another possibility isto add the term Qs· 1: (Sa-1 - sa)2to the function J (4),

a=1

whereQsis a properly chosen constant. This term also contributes to the stabilization ofthe mesh locations.

The curve resulting from the computation with fixed mesh locations, may have

oH

unnecessarily large local values of K = --. A carefulselection of the mesh locations

O

s

is required in order to avoid this effect. This condition will not be observed when computing with varia bIe mesh positions - their displacement then compensates for those large valuesof K.

5.3.6. The computer programme mayalso be used for computation with degenerate mesh points - points where the magnitude K isdiscontinuous. For those mesh points the continuity relation (9) is disregarded during the computation.

6. Piecewisecubicfunctions

6.1. Definition andparametrization

Let 11: 0 = So< S1< S2... <sa= L be a set of mesh points. The function H(s) is a cubic function in swithin each ofthe meshes (Sa-1' sa) n = 1... a, and is

continuous in its first two derivatives on (O,L).As parameters of the curve are chosen: 02H

The magnitudes Aa

=

--

(Sa),Ha and sa at the mesh points, a

=

0... a.

OS2

For the mesh (sa-I' sa) the function may then be written to

(11)

where ha = Sa- Sa-I'

This function Hls) is continuous and has continuous second derivatives as may be readily verified.

lts first derivative is represented by

in (O,L)

oH os

At the mesh location Sa' this derivative is equal to the expressions:

oH ha+l ha+l Ha+l - Ha

--(,)----A ---A

+

----os

c

_

-

3 a 6 a+I ha+1

At the mesh location Sa' oH will therefore be continuous, when the following con-os

dition is satisfied

Figure 3. The parameters of the piecewise polynornial ofdegree three ~Sa s 8-1 Ha Ha-1 Aa A.lI_1 Xa_1 Ya-1

o

MESH POINT II DATA POINT

a

KNOWN PARAMETER

'3 UNKNOWN PARAM ETER

6.2. Theunknowns andthe conditions

The values s, of the given data points coincide with some of the mesh locations SO:'

A typical interpolative curve into N data points is shownin Figure 3. The magnitudes

Aa, Ha(a =

o

.

..

2)inall the mesh points and the values s,for the given data points are still unknown. The values,for the mesh locations in between the data points are given.

The unknowns are determined from the extremum problem (4) with constraints of type (2a), (2b), (Ze) and (2d) and the additional constraints (12). The parameter QHin I (3) is set equal to O.

The condition, which specifies that the total number of unknowns must be larger than the total number of constraints, resultsin the following relationship for the total number of mesh points in the curve

a~N+q-4

with q the total number of constraints of type (2b-d).

6.3. The solution

The problem (4) can be reduced to the following non-linear equation system for the

unknowns

(13)

Ha and s, are computed by local para bol ie interpolation. The value 0 is used as initial approximation to LAp' Also in this method, the most cumbersome operation is the evaluation ofthe integrals

f

t

2

_ft

₂

Ck = Skcos H(s)ds; Sk= Sksin H(s)ds;

ti ti

No analytical expression is known to the author for the evaluation of these integrals. Consequently Simpson's rule is used for the computation, although its application is time consuming.

For further detail of the numerical computation, reference is made to Section 5, Paragraphs 5.3.5 and 5.3.6.Usually three iteration steps are required for the solution.

7. Applications of the computer programmes

The main application area of the approaches outlined in the preceding sections is road alignment (for another possible application cf. Paper 11). The principle aim in the development of the program me-package was to unburden the road-designer. With the methods existing so far a considerable number of curve parameters have had to be known a priori, in order to compute the curve. This condition is eliminated when the approaches, here documented, are used. Moreover, the designer may impose constraints of different types when computing the curve.

In actual computation several runs, usually 2 to 3, are necessary in order to obtain satisfactory results. After every run the designer re-examines the resulting curve and modifies,if necessary, his constraints. The need of this repetitive computation sterns from the fact that the designer's criteria in evaluating the curve couJd so far not be properly modeled by a mathematical expression.

After every run of the computation the following results are provided: A list of the parameters of the curve;

a list of intermediate points along the curve at regular intervals, and of the data points;

a graphical representation of the curve.

Two examples of the computations are given in the following. They iIlustrate the use of a piecewise quadratic and of a piecewise cubic function for a typical problem of road alignment, namely: Compute a smooth curve, which joins two given curve sections and passes through a set of given points (cf. Figures 4 and 5 and Table 1).

Figure 4. Computation with piecewise quadratic polynomials -48000 X >-o o o 'Cl N + >-o o o

..,

N

..

>-o o o N N +

x

0008'-- 49400.000 + 20000.000 Given data:

Straight line defined by two points PI -49000.00 P2 -48500.00 + N ()I o o o -<

Data points along the curve P3 -47555.00 P4 -47655.00 P5 -47995.00 P6 -48710.00

Circle defined by three points P7 -49212.00 P8 -49260.00 P9 -48780.00

..

N

,..

o o o -<

..

N N o o o -< +19860.00 +20410.00 +22035.00 +23000.00 +23910.00 +26280.00 +27480.00 +28750.00 +30000.00 Explanation to Figure 4:

+

given data points

*

mesh locations after computation

The figure shows the curve portion between the data points P2 and P7.

Table 1. Computation with pieeewise quadratic polynomials. List of eomputed intermediate points along the curve

Length of are X Y 0.000000 --48764.566 +20118.978 0.200000 --48630.031 +20266.966 0.400000 --48495.496 +20414.954 0.500000 ro --48428.229 +20488.948 0.600000 --48361.048 +20563.020 0.800000 --48228.763 +20713.014 1.000000 --48102.866 +20868.383 1.100000 00 --48043.770 +20949.046 1.200000 --47987.955 +21032.012 1.400000 --47885.934 +21203.984 1.600000 --47795.455 +21382.314 1.800000 --47714.766 +21565.293 1.900000 00 --47677.489 +21658.084 2.000000 --47642.093 +21751.608 2.200000 --47579.631 +21941.561 2.400000 --47533.838 +22136.157 2.600000 --47511.651 +22334.766 2.75000Q 00 --47514.701 +22484.638 2.8ooOQO -47519.981 +22534.355 3.000OÖÓ -47560.121 +22730.113 3.200000 --47623.165 +22919.837 3.400000 -47699.516 +23104.669 3.600000 --47780.152 +23287.692 3.750000 00 -47838.165 +23426.013 3.800000 --47856.321 +23472.599 4.000000 --47922.481 +23661.315 4.200000 --47979.596 +23852.971 4.400000 -48029.488 +24046.639 4.600000 --48074.050 +24241.607 4.800000 --48115.220 +24437.322 5.000000 --48154.957 +24633.335 5.200000 --48195.235 +24829.237 5.400000 --48238.018 +25024.604 5.600000 --48285.259 +25218.938 5.800000 --48338.873 +25411.606 6.000000 --48400.718 +25601.783 6.150000 00 --48453.568 +25742.152 6.200000 --48472.544 +25788.411 6.400000 --48554.342 +25970.899 6.600000 --48642.937 +26150.198 6.800000 --48734.941 +26327.779 6.950000 00 --48804.198 +26460.833 7.000000 --48827.062 +26505.300 7.200000 --48916.368 +26684.248 7.400000 --49000.570 +26865.649 7.600000 --49077.373 +27050.295 7.750000 oe --49128.681 +27191.234 7.800000 --49144.385 +27238.704 8.000000 --49199.841 +27430.829 8.200000 --49243.372 +27626.002 8.250000 --49252.371 +27675.184 ro, 00,oe mesh loeations

Figure 5. Computation with piecewise cubic polynomials

.50 Y ·40Y +30 Y

Given data:

16 _{Straight line defined by two points}

x

0 ₁₅

+ +P1 47.00 24.31 á,

"i 0

P 2 44.53 27.32 x

Data points along the curve

P 3 42.58 30.00 P 4 40.78 33.20 P 5 39.53 37.01 x _{P 6 39.23 40.00}

.

0 _1. + +P7 39.80 43.91 __, r; 0 P 8 41.41 47.62 P 9 43.86 50.70 P10 47.09 53.06 P11 50.81 54.52 P12 54.74 55.15 P13 60.00 55.24 () + + P14 70.00 54.86

.

13 o-'; 0

Straight line defined by two points

P15 80.00 54.45 P16 81.90 54.40 12 x _{11 +} 0 + + +

_'"

"'

.

0 x 10 9 2 3 x B + 0

+

•

+

..

.., 0 0 + 0 6 5 x 0 .; M Áos- Áos- ). 01;+

8. Final remarks

At present all three methods of solution to the problem are in practical use. The preferenee of the dient governs the selection of one of the methods for the solution of a particular problem.

Preferably, the number ofunknowns should always be weIl beyond the total number of constraints in the problem. Obviously it holds, that the larger the total number of mesh locations, the better the theoretical solution to the problem (4) will be approxi-mated (Ritz Method; cf. Courant and Hilbert, 1968).

In practical applications however, the engineer often requires to solve the problem with a limited number of mesh locations, particularly if he is interested in the para-meters of the curve sections. Then, of course, only a fair approximation to the theoretical solution of the problem (4) win be obtained. Also the location of the mesh points is then a more critical task and this motivated the use of variable mesh locations for the solutions given in this Paper. For a given number of meshes the variable mesh locations enable to better approximating the theoretical solution.

All three methods are programmed for a slightly extended formulation of the problem discussed in this Paper: it is also possible to compute a network of lines. Options exist to join the lines with the same angle of inclination and with the same curvature.

Studies are in progress to further refine and extend the existing computer pro-grammes. The results of these studies will be reported in due time.

References

J.H. Ahlberg, E. N. Nielsen, J.L. Walsh (1967).

The theory of splines and their application, Academie Press.

R. Courant and D. Hilbert (1968).

IT. Tbe computation of elastic beams with large deftection

This Paper documents a method ofcomputation for minimizing the energy of an idealized elastic beam constrained to pass through given data points. According to this method, the deflection of the beam is represented by its natural equation H =

=

H(s) (H being angular inclination and s length ofthe are) and approximated by a second degree mathematical spline. The numerical method, documented in Paper I, was employed in the solution of the problem.

1. Formulation of the problem

The strain energy of an elastic beam with the natural equation H = Hts), So::(s ::(sais givenby

J

~s50. ('ÖH(sY)2 dsÖS

(1)

Data points (x., yJ, i = I... n along with the angular inclinations H, at someof these points are given. We want to find that curve H'(s) of class Cl (so, Sa), which passes through the n data points, possesses the predetermined values H, at these points, and is of sucha nature the (1)is minimized. This curve will be of specialvalue when construction problems such as arising in the shipbuilding andaircraft industries must be dealt with.

2. Themethod ofsalution

The function H(s) is approximated by a second degree piecewise polynomial function, with mesh points So::(Sl ::(... ::(Sa'The mesh locations remain fixed throughout the computation. The unknown parameters of the spline are determined

from the extremum problem (I) with the constraints

f

Si+ 1

-(Xi+l -xJ

+

cos H(s)ds = 0

Si

f

Si+1

-(Yi+l -Yi)

+

sin H(s)ds

=

0; i

=

1, ... (n-l)

H(sj)

=

H, for some values j of(1 ... n)

For a detailed account of the numerical computation reference is made to Paper I.

3. Numerical example

Two typical applications are quoted below:

1. The computation of a beam passing through the point (x.; Yl) under inclination Hl =0 and passing through (x2, Y2)' The computation is carried out with 1,3, 6 and 12 meshes in order to demonstrate the convergence of the approximation to the theoretical solution (Ritz Method, cf. Courant and Hilbert, 1968).The results of the computations are shown in Figure 1 and Tables 1,2, 3 and 4.Looking at Figure 1, no difference is to be observed between the spline approximations made using 6 and 12 meshes. The largest difference between the two curves was found to be 0.009 (orthogonal distance) for s

=

135 units.

2. Computation of a beam passing through three points with the inclination at point 1 set equal to Hl

=

O.The computation is carried out with 1, 3, 6 and 12 meshes.

The results of the computations are presented in Figure 2 and Tables 5, 6 and 7. Again observing Figure 2, no noticeable difference remains between the two final runs of the computation, indicating that a fair approximation to the theoretical solution has been obtained. The largest difference between those two curves amounts to 0.30 (orthogonal distance) for s

=

130 units.

Reference

R. Courant and D. Hilbert (1968).

Figure 1. Deflection of an elastic beam, forced through two given points - 0 X ~tOO X +200 X )0-e I >-e e )0-e e ('.I I

1

+

+GiVen data:

-t-Point1:X~0, y~0, H·~0; Point 2:X~180, Y~ - 180. --j Given points - --1mesh - 3 meshes. 6 meshes. 12 meshes

+

+

x

0 -

x

OO~+ 20.000 250.000 X OOl+ e -< e o -<

,

I\) o o -<

Table I. Results of the computation for 1 mesh List of curve points

s

X Y

+

0.000 0.000 0.000 mesh location sn 0

+

10.000

₊

9.990 0.393

+

20.000

₊

19.920 1.553

+

30.000

₊

29.736 3.449

+

40.000

₊

39.390 6.048

+

50.000

₊

48.840 9.313

+

60.000

₊

58.049 13.206

+

70.000

₊

66.988 17.686

+

80.000

₊

75.629 22.716

+

90.000

₊

83.954 28.254

+

100.000

₊

91.947 34.261

+

110.000

₊

99.597 40.699

+

120.000

₊

106.898 47.531

+

130.000

₊

113.847 54.720

+

140.000

₊

120.446 62.233

+

150.000

₊

126.697 70.037

+

160.000

₊

132.608 78.103

+

170.000

₊

138.188 86.400

+

180.000

₊

143.448 94.904

+

190.000

₊

148.403 103.590

+

200.000

₊

153.068 112.434

+

210.000

₊

157.459 121.419

+

220.000

₊

161.594 130.523

+

230.000

₊

165.493 139.731

+

240.000

₊

169.175 149.029

+

250.000

₊

172.662 158Aol

+

260.000

₊

175.975 167.836

+

270.000

₊

179.135 177.323

+

280.000

₊

182.166 186.853

+

290.000

₊

185.090 196.416

+

300.000

₊

187.930 206.004

+

310.000

₊

190.709 215.610

+

320.000

₊

193.450 225.227

+

330.000

₊

196.178 234.848 mesh location SI =330.000 -I

Table 2. Results of the computation for 3 meshes List ofcurve points

s

X Y

+

0.000 0.000 0.000 mesh location so 0

+

10.000

₊

9.989 0.402

+

20.000

₊

19.916 1.589

+

30.000

₊

29.723 3.532

+

40.000

₊

39.359 6.195

+

50.000

₊

48.780 9.543

+

60.000

₊

57.946 13.536

+

70.000

₊

66.824 18.134

+

80.000

₊

75.387 23.297

+

90.000

₊

83.611 28.983

+

100.000

₊

91.480 35.151

+

109.999

₊

98.980 41.762 mesh location SI =109.999

+

110.000

₊

98.982 41.761

+

_, 120.000

₊

106.108 48.775 1 '](\ ()fV\

+

112.863 56.148 ,- ~.JV.vvv

+

140.000

₊

119.256 63.836

+

150.000

₊

125.303 71.799

+

160.000

₊

131.022 80.002

+

170.000

₊

136.433 88.411

+

180.000

₊

141.561 96.995

+

190.000

₊

146.432 105.728

+

200.000

₊

151.075 114.585

+

210.000

₊

155.517 123.544

+

219.999

₊

159.791 132.584 mesh location S2 =219.999

+

220.000

₊

159.788 132.586

+

230.000

₊

163.923 141.691

+

240.000

₊

167.937 150.850

+

250.000

₊

171.847 160.054

+

260.000

₊

175.670 169.294

+

270.000

₊

179.423 178.563

+

280.000

₊

183.123 187.853

+

290.000

₊

186.788 197.157

+

300.000

₊

190.436 206.469

+

310.000

₊

194.083 215.780

+

320.000

₊

197.748 225.084

+

329.999

₊

201.447 234.374 mesh location Sa =329.999 -I

Table 3. Results of the computation for 6 meshes List of curve points

s

X Y

+

0.000 0.000 0.000 mesh locationso 0

+

10.000

₊

9.989 0.396

+

20.000

₊

19.918 1.567

+

30.000

₊

29.729 3.489

+

40.000

₊

39.372 6.132

+

50.000

₊

48.798 9.463

+

54.999

₊

53.417 11.375 mesh locationS1 54.999

+

60.000

₊

57.968 13.447

+

70.000

₊

66.846 18.044

+

80.000

₊

75.407 23.210

+

90.000

₊

83.630 28.898

+

100.000

₊

91.501 35.064

+

109.999

₊

99.009 41.665 mesh locationS2 =109.999

+

110.000

₊

99.010 41.666

+

120.000

₊

106.154 48.662

+

130.000

₊

112.935 56.010

+

140.000

₊

119.361 63.671

+

150.000

₊

125.443 71.608

+

160.000

₊

131.195 79.787

+

164.999

₊

133.952 83.957 mesh locationS3 =164.999

+

170.000

₊

136.635 88.177

+

180.000

₊

141.784 96.749

+

190.000

₊

146.666 105.476

+

200.000

₊

151.306 114.334

+

210.000

₊

155.730 123.302

+

219.999

₊

159.967 132.359 mesh locationS4 =219.999

+

220.000

₊

159.967 132.360

+

230.000

₊

164.047 141.490

+

240.000

₊

167.998 150.676

+

250.000

₊

171.849 159.905

+

260.000

₊

175.631 169.162

+

270.000

₊

179.372 178.436

+

274.999

₊

181.236 183.074 mesh locationS5 =274.999

+

280.000

₊

183.101 187.714

+

290.000

₊

186.835 196.991

+

300.000

₊

190.572 206.267

+

310.000

₊

194.311 215.541

+

320.000

₊

198.050 224.816

+

329.999

₊

201.787 234.090 mesh locationS6 =329.999 -1

Table 4. Results of thecomputation for 12 meshes List of curve points S X Y

+

0.000 0.000 0.000 mesh location sn 0

+

10.000

+

9.990 0.394

+

20.000

+

19.919 1.562

+

27.499

₊

27.291 2.932 mesh location SI 27.499

+

30.000

₊

29.730 3.481

+

40.000

+

39.372 6.125

+

50.000

₊

48.798 9.458

+

54.999

₊

53.416 11.370 mesh location S2 54.999

+

60.000

+

57.968 13.442

+

70.000

₊

66.848 18.037

+

80.000

₊

75.410 23.200

+

82.498

₊

77.497 24.574 mesh location Sa 82.498

+

90.000

₊

83.633 28.887

+

100.000

₊

.91.503 35.054

+

109.998

₊

99.011 41.656 mesh location S4 =109.998

+

110.000

₊

'99.012 41.657

+

120.000

₊

106.156 48.653

+

130.000 _,-I _{112.938 56.000}

+

137.498

₊

117.789 61.717 mesh location S5 = 137.498

+

140.000

₊

119.364 63.661

+

150.000

₊

125.445 71.599

+

160.000

₊

131.197 79.778

+

164.998

₊

133.954 83.947 mesh location S6 = 164.998

+

170.000

₊

136.637 88.168

+

180.000

₊

141.787 96.740

+

190.000

₊

146.668 105.467

+

192.498

₊

147.849 107.668 mesh location S7 =192.498

+

200.000

₊

151.308 114.325

+

210.000

₊

155.731 123.294

+

219.998

₊

159.965 132.350 mesh location Ss =219.998

+

220.000

₊

159.967 132.352

+

230.000

₊

164.043 141.484

+

240.000

₊

167.990 150.671

+

247.498

₊

170.884 157.588 mesh location S9 =247.498

+

250.000

₊

171.840 159.900

+

260.000

₊

175.624 169.157

+

270.000

₊

179.369 178.429

+

274.998

+

181.235 183.065 mesh location SlO=274.998

+

280.000

₊

183.103 187.706

+

290.000

₊

186.840 196.982

+

300.000

₊

190.578 206.256

+

302.498

₊

191.512 208.573 mesh location S11=302.498

+

310.000

₊

194.316 215.532

+

320.000

+

198.054 224.807

+

329.998

₊

201.790 234.080 mesh location S12=329.998 -1

Figure 2. Computation of the deflection of an elastic beam foreed through three given points - 0 X +100 X .200 X >-o >-o o >-o o N I

1

---+

+

+GiVen data: Point1.X~0, Point 2:X~170, Point3:X~180, -1 Given points - -- 3 meshes - 6 meshes. 12 meshes y~0, y~ - 20; y~ - 180. H~0;

+

o -< o o -< I I\) o o -<

+

+

x

0 -20.000 250.000 X 001+

x

ooz-Table 5. Results ofthe computation for 3meshes List of curve points

s

X Y

+

0.000 0.000

₊

0.000 mesh location SO 0

+

10.000

₊

9.995

₊

0.280

+

20.000

₊

19.965

₊

1.035

+

30.000

₊

29.904

₊

2.138

+

40.000

₊

39.816

₊

3.462

+

50.000

₊

49.715

₊

4.882

+

60.000

₊

59.618

₊

6.271

+

70.000

₊

69.541

₊

7.504

+

80.000

₊

79.495

₊

8.455

+

90.000

₊

89.480

₊

8.997

+

100.000

₊

99.478

₊

9.003

+

110.000

₊

109.454

₊

8.346

+

120.000

₊

119.347

₊

6.903

+

130.000

₊

129.065

₊

4.562

+

132.999

₊

131.927

₊

3.669 mesh location SI =132.999

+

140.000

₊

138.488

₊

1.229

+

150.000

₊

147.503 3.088

+

160.000

₊

156.017 8.324

+

170.000

₊

163.950 14.406

+

180.000

₊

171.230 21.255

+

190.000

₊

177.799 28.789

+

200.000

₊

183.610 36.923

+

210.000

₊

188.625 45.571

+

220.000

₊

192.816 54.646

+

230.000

₊

196.166 64.065

+

240.000

₊

198.668 73.743

+

250.000

₊

200.321 83.603

+

260.000

₊

201.131 93.567

+

265.996

₊

201.219 99.562 mesh location S2 =265.996

+

270.000

₊

201.116 103.564

+

280.000

₊

200.318 113.530

+

290.000

₊

198.821 123.415

+

300.000

₊

196.710 133.188

+

310.000

₊

194.072 142.833

+

320.000

₊

190.992 152.346

+

330.000

₊

187.552 161.735

+

340.000

₊

183.831 171.017

+

350.000

₊

179.908 180.215

+

360.000

₊

175.857 189.358

+

370.000

₊

171.753 198.477

+

380.000

₊

167.670 207.605

+

390.000

₊

163.682 216.776

+

399.996

₊

159.865 226.014 mesh location S3 =399.996 -1

Table 6. Results of the computation for 6 meshes List of curve points

S X Y

+

0.000 0.000

₊

0.000 mesh location se 0

+

10.000

₊

9.997

₊

0.211

+

20.000

₊

19.981

₊

0.777

+

30.000

₊

29.946

₊

1.600

+

40.000

₊

39.898

₊

2.581

+

50.000

₊

49.844

₊

3.623

+

60.000

₊

59.793

₊

4.627

+

65.999

₊

65.768

₊

5.171 mesh location SI 65.999

+

70.000

₊

69.755

₊

5.496

+

80.000

₊

79.735

₊

6.126

+

90.000

₊

89.730

₊

6.408

+

100.000

₊

99.728

₊

6.231

+

110.000

₊

109.698

₊

5.486

+

120.000

₊

119.595

₊

4.067

+

130.000

₊

129.348

₊

1.872

+

132.999

₊

132.231

₊

1.04& mesh location S2 = 132.999

+

140.000

₊

138.865 1.186

+

150.000

₊

148.044 5.143

+

160.000

₊

156.781 10.000

+

170.000

₊

164.962 15.741

+

180.000

₊

172.472 22.337

+

190.000

₊

179.189 29.738

+

199.997

₊

184.993 37.870 mesh location S3 = 199.997

+

200.000

₊

185.012 37.860

+

210.000

₊

189.808 46.628

+

220.000

₊

193.583 55.883

+

230.000

₊

196.369 65.482

+

240.000

₊

198.217 75.306

+

250.000

₊

199.190 85.256

+

260.000

₊

199.365 95.252

+

265.996

₊

199.120 101.242 mesh location S4 = 265.996

+

270.000

₊

198.823 105.235

+

280.000

₊

197.649 115.165

+

290.000

₊

195.934 125.016

+

300.000

₊

193.768 134.777

+

310.000

₊

191.238 144.452

+

320.000

₊

188.431 154.049

+

330.000

₊

185.433 163.589

+

332.995

₊

184.510 166.439 mesh location Ss = 332.995

+

340.000

₊

182.322 173.093

+

350.000

₊

179.141 182.574

+

360.000

₊

175.908 192.036

+

370.000

₊

172.640 201.488

+

380.000

₊

169.357 210.933

+

390.000

₊

166.077 220.380

+

399.995

₊

162.819 229.829 mesh location S6 = 399.995 -1

Table 7. Results of the computation for 12 meshes List of curve points S X Y

+

0.000 0.000

₊

0.000 mesh location So 0

+

10.000

₊

9.997

₊

0.211

+

20.000

₊

19.981

₊

0.776

+

30.000

₊

29.947

₊

1.594

+

33.000

₊

32.934

₊

1.874 mesh location SI 33.000

+

40.000

₊

39.900

₊

2.564

+

50.000

₊

49.847

₊

3.586

+

60.000

₊

59.800

₊

4.559

+

66.000

₊

65.778

₊

5.078 mesh location S2 66.000

+

70.000

₊

69.766

₊

5.383

+

80.000

₊

79.749

₊

5.957

+

90.000

₊

89.746

₊

6.180

+

100.000

₊

99.742

₊

5.949

+

100.000

₊

99.742

₊

5.949 mesh location Sa = 100.000

+

110.000

₊

109.710

₊

5.166

+

120.000

₊

119.605

₊

3.737

+

130.000

₊

,

129.365

₊

1.572

+

133.000

₊

132.255

₊

0.767 mesh location S4 = 133.000

+

140.000

₊

138.907 1.409

+

150.000

₊

148.127 5.270

+

160.000

₊

156.901 10.058

+

165.999

₊

161.891 13.386 mesh location S5 = 165.999

+

1_I/V.VVV""7nfl"l\

₊

_{165.087 15.792}

+

180.000

₊

172.544 22.445

+

190.000

₊

179.171 29.927

+

199.998

₊

184.875 38.131 mesh location S6 =199.998

+

200.000

₊

184.886 38.125

+

210.000

₊

189.602 46.937

+

220.000

₊

193.327 56.213

+

230.000

₊

196.093 65.819

+

232.997

₊

196.741 68.744 mesh location S7 =232.997

+

240.000

₊

197.948 75.642

+

250.000

₊

198.951 85.588

+

260.000

₊

199.170 95.584

+

265.996

₊

198.954 101.575 mesh location Ss =265.996

+

270.000

₊

198.676 105.569

+

280,000

₊

197.548 115.504

+

290.000

₊

195.869 125.361

+

299.996

₊

193.726 135.123 mesh location S9 =299.996

+

300.000

₊

193.727 135.128

+

310.000

₊

191.201 144.803

+

320.000

₊

188.383 154.398

+

330.000

₊

185.363 163.931

+

332.996

₊

184.432 166.778 mesh location SlO=332.996

+

340.000

₊

182.226 173.426

+

350.000

₊

179.024 182.899

+

360.000

₊

175.790 192.362

+

365.996

₊

173.848 198.034 mesh location S11=365.996

+

370.000

₊

172.553 201.824

+

380.000

₊

169.333 211.291

+

390.000

₊

166.119 220.761

+

399.996

₊

162.907 230.226 mesh location Sl2 =399.996 -1

111. Approximation of measured data by piecewise bicubic polynomial functions

This Paper documents the use of piecewise bicubic polynomial functions for

smooth surface approximation of measured data. Piecewise polynomials are selected

for approximation because they combine the advantages of simple computation and

high flexibility. The programme system which is documented in this Paper, was so

far applied mainly to problems in geology and geodesy.

Statement of theproblem.

Measured values Z, are given in a set of discrete points (Xi'Yi;i = 1... n) in a

region G(a :::;;X:::;;b, c:::;;y:::;;d). The measured values Z, are erroneous. It is the

task to compute a continuous and smooth approximating surface u(x, y),so that

the differences ei =u(xi, yJ - Z, are in the order of the

pregiven measuring accuracy of the values Zi. (1)

Solution of the problem

As approximating function a piecewise bicubic polynomial function is chosen.

For this purpose the region G is subdivided into r-s rectangular regions

Hap (I;a :::;;X:::;;I;a+l' 11p:::;;y:::;;11p+1;1;0= a, I;r = b, 110= c,115= d);*

a = O ... r, ~ = O ••. s.

A piecewise bicubic polynomial function is then that function u(x, y) which has

continuous firstderivatives in G and which equals a polynomial of the form

~ ahk(x- I;a)h (y- 11p

t

(h,k, = 0, 1,2,3) within each ofthe intervalsHap.

The coefficients ahk may be expressed as function of the magnitudes u (1;,11),

u, (I;,11),uy(I;, 11),uxy(I;, 11)at the four mesh points of HaP which are chosen as

para-meters ofthe function u.Itholds

[ahk]= [dhrn(Ax)]. [brnn]·[dnk(Ay))T;h,k = O ... 3; m,n = 1, ... 4 (2)

where [ . ]stands for the symbol of a matrix,

* In practical computation the magnitudes rand sare, at this stage chosen empirically, by inspeetion

[bmn]

=

[U"_I'P_l , lly,"-I,P-l , U"-I,p , lly,"-I,p ] llx,,,-1,P-l' llxy,,,-1,p-1' Ux,,,- I.P -1> llxy."-I.P U",.P-l , lly'''.P-l , U"'P , Uy.".P llx.<x.P-l , llxy.".P-l , UX''''P , Uxy.".P [dhm(ilx)]

= [

~

- 3/ilx2 us,»

o

1 - 2/ilx l/ilx2

o

o

3/ilx2 - 2/ilx3

Itcan be readily proved, that with this definition of ahkthe function u is continuous

and has continuous first derivatives in G.

The magnitudes bmn are still unknown. Let us denote the total set of unknowns

by b. Those unknowns bare now computed from the condition (1) and from the

additional condition, that u is a smooth function in G. Consequently the following extremum problem may be formulated

I:

=

p.

7

(z,

-

u(xi,

v.

:

b»2

+

IIf(U, u., uy, u.., Uyy,UXy)dxdy _. min (3)

G

The first term in this expression I is a measure for the goodness of fit at the sample points and the second term is a measure for the smoothness ofthe function u. Various

possible choices exist for f; one can try to model with f the physical situation under

consideration or alternatively f can be chosen just to ensure a general smoothness ofthe surface, with e.g.

(4)

The weight factor Pis chosen empirically. The task is to select a value for Psuch

that the residuals ei correspond in magnitude to the measuring accuracy of the values

Z. If consequently after one run of the computation the values eiare too large, the

weight factor P is increased and the computation is repeated with th is new value of P.

A more detailed theory on the selection of P was given in Kubik (1969).

It is possible that the residualscannot be made sufficiently small merely by varying

P. If so, G must be subdivided into a larger number of subregions H"o, after which

the computation is repeated.

From the extremum problem (3) there results an equation system for the unknown parameters b, ofthe form:

81

-=0

8b (5)

It is recommended to choose f as aquadratic function in b, so that the corresponding equation system is already linear. Otherwise one of the known methods of numerical mathematics must be applied to solve(5).

In our applications the system(5) was sofar always linear and it is solved for the

unknowns b by the Gauss' elimination technique. The generation of the equation matrices is outlined in the following:

The functions u(S,TJ), u.(S,TJ), uy(S,TJ)are linear in b, and consequently for the smoothness condition of type(4),I (b) is aquadratic function in b, wich may be written to

I(b) = !bTMb

+

bTW-4 min

The matrices Mand w may be evaluated by substituting ahk,(2), into (3). The extremum problem is solved by equating the first derivatives of Iwith respect to b to zero, viz.

81

-- = Mb +w = 0 8bT

In the computer programme the arrays M,w are generated by accumulating the

contributions of the smoothness condition and of the individual data points. The contributions to be added to M due to the smoothness condition (4)for I meshHaP are

column of bm,n

3 ~

row Ofb;,j 1: dhm(~x)dnk(~Y)(Ahkh'k'

+

Ah'k'hk)dh,;(~x)dk'j(~Y) h,k,h',k'=O

bm•n., , the unknowns at the four mesh points of meshHaP

{

hk

Ahkh'k'= (h

+

k- 1)bh'

+

k '

+

I) if denominator >0 if denominator = 0

There is no contribution to w.

Contribution to M

column of bm,n

Contribution to w

In this way Mand ware generated.

The dimension of this equation system is four times the total number of mesh points involved,

The problem was programmed in Simple code on an Electrologica Xl computer configuration (10K care memory). The programme system has been used in production since June 1969.

Examples of application

The measured values Z, of lead concentration at the points (Xi,yi) are shown in Figure 1.Itis required to approximate these values by a trend surface.

In a first series of computation the function f was chosen as follows

The region G was subdivided into 5 rectangular regions H of equal size. Computa

-tions were carried out with the weight factors P = 100, P =1 and P = 0.01 respectively. The computed regression surfaces are shown in Figures 2, 4 and 6. The residuals ei were mapped for every regression surface, see Figures 3, 5 and 7. From those results it can be seen, that with smaller values of P smoother regression surfaces are obtained. The residuals ei thereby increase in magnitude, though in this particular example for an amount which is hardly significant.

The other parameters of the computation were identical to the case above. The computed regression surfaces and the map of residual errors ei are shown in the Figures 8 to 13. The same conclusions hold as above: smaller P values result in smoother approximation surfaces.

A comparison of the Figures 2 and 4 of the fi.rstseries of computation with Figures 8 and 10 from the second series demonstrates that there is no significant difference between those corresponding figures.

It follows that in this particular case with a proper value of P the problem is in-dependent of which of the two functions is chosen for f. Apossible explanation of this effect may be the dense arrangement ofthe sample points in the area.

Only for small values ofPthe solution depends on the chosen function f (cf. Figures 6 and 12).

References

K.Kubik (1969).

The estimation of the weights of measured quantities within the method of least squares.

Figure 1. The measured valuesof lead concentration at the data points 0" 036 0 3J 086 052 084068Ö0'04 0288 0124 0'29

ê

Z0156 1460 040 ~ 0125 0 0'95 0'50 4tl ·25 '00 0'05 0195 0'50 o °S5 0155 J<O o050~~: 100 0120 0200 1'50 0 01400 125

14

0

%

°0185 0200 '20 0'35 0225 0225 t23010 0250 0175 0220 080 0200075 0'00 0140 11000180 0170 0'00 24000210 0140 0115 0135 1250 185 o 40 00 200 ~ ttz ree o 0080 o 60 180 070 14500210 0'50 00°175 126 160 0" 0160 0012 " 060 035035 0'28 oao 40 o 0" 3480 0" 040 0" 080 0'80

8

"

" 0" 0" 0" 058 360 1140 116 0'48 05~ o 0120056 20400352 0" 4000128 0" 036 0" 0140 1720 090 0180 054 040 49

8

0,,4 040 42 064 O,?O'O 0" 052 020 O~240 0148 07942 13'_o 0.0280 °2032030116 030 136030 0'44 1040 058 40000

°

90022 048 80_1j44 048 022 05!! 060 Oso 044 067 960 35ct>192 0172 0" 050 048 011802050560 0" 2400.40 0" 012 028 0" 040 040 360 0" 040 022 ₀₄₀ 040 0" 0'6042 0" 48 48 036 000 028 036 56 048 40 034044°036 020 @zo." 020 0" 0811

8

220 " 0'16

Figure 2. Regression surface computed witht the weight factor P =100(contour intervalequals 25 units) -8000X -7000 X -6000 X >-o o o

'"

(") + 5555555 >-o o o

..,

(") + :H .7 :l.d.~ 1J 33.33.3.3 5555 77.7.7.7 97 " .d~.1.1 33.3 5555.5555555555555 .777.7.7.7 ~:l·!977~.d.d/l::.1~.3 55~5.5555 5.55.5.5

n

~·n}i

:

h

~4-97: :l/.1.133.3 ~55 77~~I~7.77}m:~~7.7.7 .7.n·u 97 :i.4-.1 3 55 77.7.7 .7777.7.7 7.7.7 ;r":i.~.1 .355 7.7.7 777.7 7.7 9:;1 fo:l.~ .1 .33 55 7.7 99.9 .777.7.7 555.5 .7

~

77

:

:l_/

1

:

:

3 ~ ;I

;

9~~~·~.~~9

.

~;-j:

;

5~~55~ :;

~7

ij ~

:

!:

:

3

-11

7;:;

·

gm

~

.

~

+

:

~

;;}~ ~~~ + ij:l H .1 3.3 55 77.7 .9999 77.7.7 5555555555 é :i.~ .1.3.3 55 77.7 7.7.7 55.5.555.5555 1B:i .4-.1.1.3.3 55 77.7.7 777.7 5555555555 p~:i.~ .1.1.33 55 .777.7.7 777.7.7 55555.55555 p:i~.!.4- .1::·1 53 5~55 .7H;:H7.7.77.7.n:~:n·7 ~~~~~ ~~ 1.1 333 555 77.7.7.777.7.77.7 555555555555 4- 1.1.1 33.3 5555 .777.7.7.7 55555 555555 1.11.1 33.3.3 55.5.55 55555 55555 111 11.1 33.3.33 555555 555555 55555 1111.1 333.3.3 5555555.55 55555555 55555 33333.33 555555555.5.55555555 55555 P333 ~~5m~g3.3 555.5f'~~.55 ~ P333333.3333.333333333.33.333 555.5.5 33333.33.3.3333.3333.33.3.33.33.3.3 55555 333333333.33333.3.33.333.33333.33 555.5.55 .5 333.333.3333.3.33.3333.3333333.33333.3 55555 555 3333333.333.333.33333333333333333 .55555 5555555 3333.333.33333333.33333333.3.3.3.3 55555 5555555.5.5 3.33.3333.33333.33.3.33.33333333 5555 5555555 3333.333.3.3 33.333333 5.555 555555 3333.3.33 33.333.3 555 77777.77.7.7.7 55555 3333.3.3 333.3.3 55 .777777.7.7.77.7.77.7 5555 3333.3 33333.3 5 777.7 77.7.77.7 .5555 3333.3 3333.3 ;;:; -f.m:~ ~~~5 ~53.3 1 355353 .77.7 .777.7 555 333 3.333.3 .77.7 77.7.7 555 33.3 333.333 77.7 777.7.7 555 .333.3 333333 .7777.7 .777.7.7 555 33.3.3 3333.3 ~ 7777.77777.7777 5555 33.3.3 33.333 ~ 7777777777 55.55 33333 .333.3.3.3 "'5 7777 555.5 3333 33.333.3 "'55 55555.5 :1·m33 3~3533 "'~~~~555555~~~~~555 33.3.3.3.3 33.33.3 333.3.33.3.3 33.333.3 33.333.3333.3 .3333.3 3333;uJ3.33.333.3 + 3333.3+ 33.333.33333.3 333.3 3333.3.3.3 333.3.3 p~53.3 111.1.::::::::: :.::: :.:1.1.1 ~m g3.3 .111.1.1.1 .111.1.11.1.11.1.1 33.33 p5 1.:: 11 1: : : : : ::1 .1 553 .3 111 ·H.H.H.H .1111.1.11 111.1.1 333 .3 1.1 -I-H-I-.H.H.H..I· 11.1.1.11.1 .1111.1.1.1 33.3 .3 1.1 -I-4-.H.H.H.H.-I-·H .111.1.11.1 .1.1111.11.1.1 33

33 1.1 ~-I-.-I--I-.H.H-I-.-I--I-.4--I-.-I-11.1.11.111.1.11.1.1.111.1.11.1.133.3

33 1.1 -I-H-I-.H.H-I-.H.-I--I-.4- 11.1.11.111.1.11.1.1.111.11.1.11.13.3

~ 3 1.1 ~-I-.H.-I--I-.H...~-I--I-~.-I-.4- 1.11.11.1.11.1.11.1.1 11.11 3 ~ 33 1.1 H"~.H.H.H.-I-.-I- 11.~ 1.11.1.11.1.1 .11.4 3 ~5 3.3.11.1 ~HHHH 11.1.1'1.1.11.1 11.1

5 .3.3 1.1.1 .11.1 ~~~.-I-.-I- .11.1

5 33 11.1.1 H..H-I-.-I-·H 11.1

7 55 .33 .111.1.1.1.1 HH~.~-I-.-I--I-.~-I-.~11

P77 77 5\5g33 3333.111.1::::: 11 H!!!H ~!-!.~ 99 77 555 333 11111 HH4- :l:ld.dd.d -I-U 1 999 777 55 333 111 HH :i:i:i:i:l.:l:ld.d:ld.:l .~ f,,,~~~1 99 n~7 o~5 ~ 3~" 11 H-I- ~~~666 .... ;id.d~'l'l >-o o o (") (") + >-o o o N (") + >-o o o (") + X0008- X OOOL- X 0009 -+

,

W o o o + W

'"

o o o -< + W

..

o o o , W W o o o + W N o o o

Figure 3. The residuals at the data points for the regression surface computed witht the weight FactorP = 100(lmm equals 19 units)

..

1iI. ...

..

? ... \I" \I

•

.,.

•.. ".,\1 Jo

.,

Figure 4. Regression surface computed with the weight factor P units)

1(contour interval equals 25

-8000 X -7000 X -6000 X 333.3 555555.5 5555555555.5555 .333.3 5555555 5555.555.5555.55 1 33.3.3 55555.5 777.7.7 5555555555.55 1 333.3 555555 7.77.7.7.777.7.7.7 555555555555 1 33.3.3 5555 77.7.77.7.77.7.77.7.77.7 55555555555 11 33.3.3 5555 7.77.7.7.777.7.777.7.7.7 55555555555 11 .33.3 555.55 .77 7.7.77.7.77.7.77.7.77 7.7 55555555555 1.1 33.3.3 555.55 .777.7.77.7.77.7.77.7.77.7.7 5555.5555555 11.1 .333.3 555.5 77.7.77.7.777.77.777.77.7 555555.5555.5 o- 111 .33.3.3 5555 77.7.77.7.77.7.777.7.77.7 5555555555.5 + 0 11.1 .333.3 55555 77.7.7.777.77.7.777.7.7.7 55555.55555 w 0 .1.1 .333.3 +5555 77.7.77·7.7~t777.7.7.7 5555555~ U' 0 ₀

'"

11.1 .333.3 55555 7 7.7.77.7.77..77.7.77.7 55.55555 0 (') 11.1 33.3.3 5555 77.7.77.7.777.7.777.7 55555.555.5 0 + 11.1 3.33.3 555.55 77.7.77.7.77.7.77.7 55555.5555 11.1 33.33.3 55555 77.7.77.7.77.7.7.7 55555555 11.1 .333.3.3 55555 77.7.77.7.7.7 55555555 1.1 33.3.3.3 555555 7.7 555555555 1 .3333.3 5.55555 555.55555 .3333.3 555555 555555555 3.3333.3 555555 5555555.5 33333.3 55555555 55555555 333333.3.3 5555555555 5555555.55 >- 333333.3.3 5555.555555.55555555555.555.5 ₊ 0 3.33333.3.333 55555.55555555.55555555 w 0 33333333.333~3.33 555555~5555.555.5555 ₊ C>. 0 ₀ -er 33333.333.333.3 33 55555 555555555 0 (') 33.3.33.3.33.3.33.3333 55555555 0 + 3.3.333.3.33.3.33.3.3 .33333.3 >-o o o (') (') + 5555555555555.5 333.3..33.3.33..33.33 5555555.55555555555 333.33..3.33.3.33.3.3.3.3 555 5.55555.555.5 33.3.33..3.33.3.333.3333.3 5555555 33.3.33.33.3.33.3.33.333.33.3 55555.5 33.3.333..3.33.3.33.3.333.333.3 5555.5 33.3.333..3.3 .3.33.33333

~~~

~3~~

.3~~~~~~

+

5555 3.33.3.3 33.3.333.3.3 5555 3.33.3.3 33.333..33.3 5555 33.3.3 .333.3.33.3.3 5555 33.3.3.3 33333.333.3 5555 33.33 33.3.33.333.33 555.5 333.33 33..333.3.33.3.3 5555 333.3 3333.3333.3.3.3 555.55 3333.3 333.33.333.333 55555 3333.3 333.33.3.33.3.3 555555.5 33.3.3.3 33.3.333.33.3.3 5555555555555.5 33.33.3.3 33..3.3.333.3.3 555555555555 .333.33.3 .1 333.333.3.3.3 5555.55555555+ 33.3.33.3 + 1.1.1 .333.3.339.3.3 3.333.3.3.3 .1111 33.3.33.33.3 .33333.3.3.3 .1111.1 33.3.33.3.3.3 .333333.33.3.3 1111.11.1 333.33.33 33333.3.33.3.3.3 111.1.11.1 33.3.33.3.3 3.33333.3.3.333.3 .1111.1.11.1.1.1 33.333.3 .33333.3.3.3 111.1.11.1.11.1.11.1.1 .333.3.3 3333.3 111.1.111.1.111.1.11.1.11.1 .333.3 3.3 111.1.11.1.111.1.11.1.11.1.111.1.1.1 .333.3 11.1.11.1.11.1.11.1.111.1.1111..1.11.11.1.1 3.3.3 11111.11.1.11.1.111.1.11.1.111.1.111.1.11.1.1 3.3 .11111.1.111.1.11.1.111.1.11.1.111.1.111.1.11.1.1 .3 11111.1.111.11.1.11.1.111.1.111.1.11.1.111.1.11.1 .3 1 1111.1.111.1.11.1.111.1.11.1111.11.1.111.1.11.1.1 .11 1111.1.11.1.\.1..1.1.11.1.11111.1.1.11.1.11.1.\.1..1.1.1 1.11 11.1.11.111.111.1.11.1.11.1.111.1.11.1.111.1.11.1.1 lil 11.11.1.111.1.11.1.11.1.111.1.11.1.11.1.111.1.11.1 1.11 1.1.11.1.11.1.11.111.1.111.1.11.1.11.1.111.1.11.1.1 1.11 11.11.1.11.1.111.1.111.1.11.1.111.1.111.11.1.1.11.1 111 11111111111111111111.11.1.111.1.11.1.111.1 lil 11111111111111111111111111111111111 111 1 1111111111111111 111 1111111111111 11 1111111" 11.1 .77777.7

7;

;m;r~;-1j

77777.77.7.77.7 77777.777.7.77.7 7.7777.7.77.7.77.7 77777.7.777.7 7777.7.777.7 7777.7.7 c-o o o N (') + >-o o o (') + X0008- X OOOL- X0009 -+ + w o o o -< + w w o o o + W N ei o o -<

Figure 5. The residuals at the data points for the regression surface cornputed with the weight

factor P = 1(Imrn equals 19 units)

..

-,

"

\

lil lil <,

"

"

"

..

~ "J,.

•

~lil

,

<l

"

,,

~

....

....

"

.

.."

..

•

Figure 6. Regression surface computed with the weight factor P 25 units)

0.01 (contour interval equals

-8000 X -7000 X -6000 X o o o '" M + + 5555555555555555555555555555 55555555555555555555555555~5 555555555555555555555555555 55555555555555555555555555 5555555555555555555555555 555555555555555555555555 55555555555555555555555 555555555555555555555 5555555555555555555 55555555555555555 55555555555555 55555.5fl5555 555555 + + w o o o -< + w (J' o o o -< >- ₊ 0 w 0 + "'" 0 + + 0 " 0 M 0 + _-< 3 3333.333.333 333333333.33333 33333333333.3333333 333333333333333333333 ~33333333333333333333333 333333333333333.333333333333 33333333333333333333333333333 33333333333333.3.3333333333333333 3333333333333333333.33333333333333 + 333333333333a333333333333333333353.3 3333333333333333333.33.333333333.333333 3333333333333333333333.333333333333333 33333.3333333333333333333333333333333333 333333.333333333333333.33.33333333333333333 333333333.3333.33333.333.3333.33333333333333333 3333333.33333.3333333.3.333.333.333.333333333333333 3

3333§~§§§§~3§§~

§

§~~~§~~§~~~~~~mm§

g~~g~~5gg~~g~

~

3333333.33.3333.333.333.333.33.3.3333.3.33.3.333.3.3.33.3.33333333333333 3 3333333333333333.33333333.33.3.3333333.3.3333.3.33.33333333333333 3333333.33.333333.333333333.33.3.33333333.33 333.3333333333333333 333333.33.33333333.33333333.33333.3.333.3333333333333333333333 3333333.33.333.3333333.3333.333333.3.3.3333.33.3.333333333333333333 ~3333333.333.33S3.3333.3.333.333.33333.3353333.3.33.333333333333a33 333333.333.33.33333.3.33.3.3.3333.3333.33.3.33.33333.3.33.33333333333333 33333333.33.33333.33.3.33.333.33.333.3.33.333.33.3.33333 33 333333333 33 ~~333333.3.333.33.3.3.33.3.33.3.333 3.333.3333 3.3333.3.33.33 333 33333.3 33 33 ~333333.33.3.33.3.3.333.3.333.3333.3.33.3.3.333333.33333.333333333333333

p~555§5~5§§~55~55m5555555555~5355n55555

~~~g3~~3TI~g~3

~3333333333333333333333333333333333333333333333333333333 ~3333333333333333333333333333333333333333333333333333333 ~333333~~~~~~~~~~~~~~~~~3333333333333333333333333~~~~~~ >-o o o M M + >-o o o N M + >-o o o M + + + X 0008- X OOOL- X0 009-+ + w w o o o + W N o o o -<

Figure 7. The residuals at the data points for the regression surface computed with the weight factor P =0.01 (1 mm equals 19 units)