UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Interactive spline approximation Merchant, Marian 1974

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1974_A6_7 M47.pdf [ 5.39MB ]
Metadata
JSON: 831-1.0051769.json
JSON-LD: 831-1.0051769-ld.json
RDF/XML (Pretty): 831-1.0051769-rdf.xml
RDF/JSON: 831-1.0051769-rdf.json
Turtle: 831-1.0051769-turtle.txt
N-Triples: 831-1.0051769-rdf-ntriples.txt
Original Record: 831-1.0051769-source.json
Full Text
831-1.0051769-fulltext.txt
Citation
831-1.0051769.ris

Full Text

INTERACTIVE SPLINE APPROXIMATION by MARIAN MERCHANT B.Sc, Simon Fraser U n i v e r s i t y , 1971 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science i n the Department of Computer Science We accept t h i s thesis as conforming tp the required standard THE UNIVERSITY OF BRITISH COLUMBIA January, 1974 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f the r e q u i r e m e n t s f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h Co lumb ia , I a g ree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s tudy . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s . It i s u n d e r s t o o d tha t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i thou t my w r i t t e n p e r m i s s i o n . The U n i v e r s i t y o f Br Vancouver 8, Canada Depa rtment - i i -ABSTRACT The use of s p l i n e basis functions i n s o l v i n g l e a s t squares approximation problems i s investigated. The question as to which are appropriate basis functions to use i s discussed along with the reasons why the f i n a l choice was made. The Householder transformation method for solving the f i x e d knot s p l i n e approximation problem i s examined. Descriptions of both an automatic procedure using function minimization and an i n t e r a c t i v e procedure using a graphics terminal for solving the v a r i a b l e knot s p l i n e approximation problem are given. In conclusion, numerical r e s u l t s using the i n t e r a c t i v e system are presented and analyzed. - i i i -TABLE OF CONTENTS Section 1 Splines i n I n t e r a c t i v e Approximation 1 Section 2 The Spline Representation Problem 3 2.1 Introduction 3 2.2 D e f i n i t i o n of a Spline Function 4 2.3 Piecewise Continuous Polynomial Spline 5 Representation 2.4 Mathematical Spline Representation 6 2.5 B-Spline Representation 9 2.6 Example of B-Spline Representation 12 2.7 Equivalence of Spline Representations 14 Section 3 The Linear Approximation Problem 18 3.1 Introduction 18 3.2 D e f i n i t i o n of the Linear Problem 19 3.3 Solution of the Linear Problem 20 3.4 Solution of the Fixed Knot Problem 2 7 Section 4 The Variable Knot Approximation Problem 30 4.1 Introduction 30 4.2 D e f i n i t i o n of the Variable Knot Problem 31 4.3 Solution of the Variable Knot Problem 32 4.4 Knot Optimization 34 Section 5 Numerical Results 37 5.1 Introduction 37 5.2 Use of the System 39 - i v -5.3 Titanium Heat Data 52 5.4 The Bug 62 Section 6 Summary and Future P o s s i b i l i t i e s 669, Bibliography 72 Appendix A Program L i s t i n g s 74 B Users' Guide 102 - v -TABLES I Test Data - 2 Non-uniform Knots 44 II Test Data -,2 Knots Optimized 45 I I I Test Data - 3 Uniform Knots 47 IV Test Data - 4 Non-uniform Knots 50 V Titanium Heat Data - 5 Uniform Knots 53 VI Titanium Heat Data - 5 Non-uniform Knots 56 VII Titanium Heat Data - 5 Knots Optimized 59 VIII The Bug 64 - v i -FIGURES 1 Elementary Cubic Spline Basis Functions 8 2. Cubic B-spline Basis Functions 15 3. Test Data - 2 Uniform Knots 43 4. Test Data - 2 Knots Optimized 46 5. Test D a t a - 3 Uniform Knots 48 6. Test Data - 4 Uniform Knots 49 7. Test .Data - 4 Non-uniform Knots 51 8. Titanium Heat Data - 5 Uniform Knots 55 9. Titanium Heat Data - 5 Non-uniform Knots 58 10. Titanium Heat Data - 5 Knots Optimized 61 11. The Dented Bug 6 3 12. The Bug 68 - v i i -ACKNOWLEDGEMENT There are many people who aided i n the development of t h i s t h e s i s , Since I cannot hope to include a l l i n d i v i d u a l s who provided help; I s h a l l mention only; those who provided major assistance and apolo-gize to anyone who may f e e l s l i g h t e d by being omitted. I would l i k e to thank Dr. J . M. Varah who provided academic assistance through our. discussions of the t o p i c . Also I am g r a t e f u l f o r the f i n a n c i a l aid which he provided i n the form of a reasearch a s s i s t a n t s h i p . I would l i k e to thank the Department of Computer Science and the Computing Centre, for the assistantships they provided. In p a r t i c u l a r I appreciated the IBM fellowship which I received for research i n Inte r a c t i v e Numerical Analysis. This provided the i n i t i a l ideas for t h i s t h e s i s . F i n a l l y , I would l i k e to thank my husband Peter for a l l the dishes he did during the w r i t i n g of t h i s t h e s i s . - v i i i -NOTATION The following common notation i s used throughout the thesis: [n], [n,pp. ] x € [a,b] x C(a,b) s CC m[a,b] s(x) dx J a) r(x) refe r s to reference n i n the bibliography; a <_ x <^  b; a < x < ,b; s has m continuous derivatives on [a,b]; i s the j - t h d e r i v a t i v e of s with respect to x evaluated at 6; i s the f i r s t d e r i v a t i v e of the function to with respect to x; A = {6.:i=l,2,...,k} i s the set A with elements 6,6,...,6 ; {(?£,y£):£=1,2,...,n} i s the set of orderd. pairs (x^,y^); £i i s the vector a; that i s , a = min f (a^ ) a. max {c,d} m i s the le a s t squares norm of aj i s the minimum of the function f with respect to a; i s the maximum of c and d; i s the transpose of the matrix Q . - 1 -Section 1 SPLINES IN INTERACTIVE APPROXIMATION Polynomials are frequently desired as a set of basis functions for approximation. Problems i n obtaining accurate r e s u l t s with,the 2 IB. standard set of basis functions {1 } lead to the development of orthogonal polynomials as the set of basis functions. The use of orthogonal polynomials i s a well-established technique i n approximation. Consequently, they have been studied i n depth and are not discussed further here. Polynomials do have one disadvantage i n approximation. That i s , t h e i r nature over the e n t i r e region of approximation i s determined by t h e i r behavior i n only a small area of t h i s region. Higher-order polynomials do not a l l e v i a t e t h i s problem but merely impose more o s c i l -l a t o r y behavior on the approximation. However, polynomial splines can counteract t h i s r e s t r i c t i v e nature of polynomials. Polynomial splines, consist of piecewise polynomials connected at points known as knots over the region of approximation. Each piecewise polynomial determines the shape of the approximation i n a small area r e l a t i v e l y independently of surrounding areas. The amount of dependence i s determined by imposing co n t i n u i t y requirements at the knots. It i s p r e c i s e l y the involvement of these knots that makes splines i d e a l l y suited f or i n t e r a c t i v e approximation. Previously, automatic pro-cedures in v o l v i n g the minimization of the l e a s t squares error were used - 2 -to determine the best locati o n s for the knots. As the knots occur non-l i n e a r l y , appropriate i n i t i a l guesses had to be made and there ex i s t s the p o s s i b i l i t y that the procedure would converge to a l o c a l minimum rather than the global minimum. Such a l o c a l minimum might not have been a suit a b l e approximation. Graphical i n t e r a c t i o n allows immediate contact with the s o l u t i o n procedure. The locations of the knots can be chosen v i s u a l l y and the approximation attempted with t h i s set can be observed. The i n t e r a c t i v e procedure enables poor i n i t i a l estimates for the knot locations to be eliminated. Also i t allows for the manipulation of the knots u n t i l a s a t i s f a c t o r y approximation i s obtained. Using the knots obtained from the i n t e r a c t i v e procedure as i n i t i a l guesses, an automatic procedure should converge r a p i d l y to a sui t a b l e minimum. - 3 -THE SPLINE REPRESENTATION PROBLEM 2.1 Introduction There i s no u n i v e r s a l l y accepted representation f o r s p l i n e basis functions. There are several known^representations each with i t s own p a r t i c u l a r advantages and disadvantages. Carasso and Laurent 14] discuss three methods of numerical con-s t r u c t i o n of splines - a proj e c t i o n method, a method of d i r e c t r e s o l u t i o n and a method using a bas i s . Of these, they recommend the use of a method inv o l v i n g a bas i s . With the choice of a reasonable b a s i s , Carasso and Laurent conclude that t h i s method provides more accurate.results than the proje c t i o n method and three times l e s s computation than the method of d i r e c t r e s o l u t i o n . G r e v i l l e [8] provides a comprehensive overview of basis functions for s p l i n e s . From the d e f i n i t i o n of a s p l i n e function, he develops a representation-using truncated power functions (discussed i n Section 2.4) and one using B^splines (discussed i n Section 2.5). de Boor and Rice [5] summarize these representations and also include a representation i n v o l v i n g piecewise continuous polynomials (discussed i n Section 2.3). Schultz [17] gives a general basis f o r B^-splines. In [18], Schultz describes the representation f o r cubic B-splines i n more d e t a i l . The bas i s function r e s u l t i n g from applying the set of cubic B-rsplines to the s p e c i a l case of uniformly spaced knots i s stated. The d e r i v a t i o n of t h i s r e s u l t i s given i n d e t a i l i n Section 2.6. - 4 -2,2 Definition, of a Spline Function Although splines e x i s t i n engineering and d r a f t i n g as a device for curve smoothing, the'.basic mathematical f ormulationcof a s p l i n e function comes from piecewise continuous polynomials. The mathematical d e f i n i t i o n formalizes the engineering concept. D e f i n i t i o n 1;' A (polynomial) s p l i n e function of degree m on [a,b] m-1 i s a polynomial of degree, m which Is i n C [a,b] . Although t h i s d e f i n i t i o n incorporates the basic notion of a s p l i n e function, i t does not provide the e s s e n t i a l components needed for the use of splines i n numerical problem-solving. For t h i s purpose, the following, more constructive d e f i n i t i o n of splines i s better. D e f i n i t i o n ,2;' Given a p a r t i t i o n a = 8n < 6_ < ....'< 6.. < 6\ ,. = b —•• — — 0 1 k k+1 then a (polynomial) s p l i n e function of degree m with k i n t e r n a l knots S-^ j ' '"' ^k ° n fa»^-' ^ s a function S(x) with the following properties; 1. S(x) i s a polynomial of degree m or le s s i n 2, S(x) and i t s derivatives of orders 1, 2, m^l are continuous everywhere. Let A = {S^ ; i = 0, 1,. k+1} be the set of knots and S(x) = (s^(x) ; i = 0, 1, ..., k} be the set of polynomials such that s.(x) i s i n [ 6 . , 6 . . - ] for i = 0,1, k . The set S(x) must x i ' x+1 • ' s a t i s f y the following continuity conditions at the knots: 7j Sl-l(x) - 5 -= — . s . (x) 6 . dx3 1 ' 6 \ l l for i = 1, 2, . .."k; j = 0, 1, ..., m - 1 . 2.3 Piecewise Continuous Polynomial Spline Representation The piecewise continuous polynomial d e f i n i t i o n (2) can be formulated into an approximation problem. Suppose that on each i n t e r v a l [ 6 \ , 6 ^ + ^ ] ' - for i = 0, 1, k; the data i s approximated by a poly-nomial of degree m or l e s s . Given the set of c o e f f i c i e n t s {c_^ _. : i = 0, 1, k; j =1, 1, .. . , m} the problem i s to f i n d the c. . ' s where m S(x) = s.(x) = J c.Cx-S.)3 for 6 . < x < <$.,.. where i = 0, 1, .... k . .1 — — l + l This system r e s u l t s i n (m + 1)(k + 1) = mk + m + k + 1 unknowns c.. which i s mk more than are required for a non-redundant s p l i n e i j representation. Therefore i t i s only necessary to compute the set { c _ : i = 0 and j = 0, 1, ..., m and i = 1, 2, ..., k and j = m} . The remaining c o e f f i c i e n t s can be computed from the constraints derived from the continuity conditions; that i s 1 d J , . c . . = —, — . s . . (x) 1 J j ! dx 3 for i = 1, 2, k; j = :0, 1, m - 1 6. l - 6 -This system i s useful for approximation as i t gives an e x p l i c i t representation for each piecewise polynomial between each knot p a i r . But despite the f a c t that the basis functions are defined over a p a r t i c u l a r knot i n t e r v a l , they must be computed over the e n t i r e i n t e r v a l . Also, the system of equations formed tend to be i l l - c o n d i t i o n e d ; that i s , the r e s u l t i n g s o l u t i o n i s not accurate. Since t h i s representation i s i d e n t i c a l to the mathematical representation between each knot p a i r , i l l - c o n d i t i o n i n g occurs for the same reason (as described i n Section 2.4). 2.4 Mathematical Spline Representation The standard representation of splines i s that of elementary splines ( a l i a s truncated power functions). This representation i s used mainly i n mathematical analysis. Most theorems i n v o l v i n g s p l i n e functions are derived and proved using elementary splines as they are easy to manipulate a n a l y t i c a l l y . D e f i n i t i o n 3: An elementary s p l i n e function of degree m, y™ , i s defined by y™ for y > 0 m y - 7 J y+ = \ 0 f o r y < 0 Elementary splines give r i s e to a set of basis functions for s p l i n e s . In p a r t i c u l a r the set {1, x, x 2, x m, (x-6^)™ , (x-6 k>™ } forms a set of basis functions for a s p l i n e of degree m . An example of these basis functions f o r a cubic s p l i n e with four uniformly spaced i n t e r n a l knots i s given i n Figure 1 . - 7 -These basis functions can also be formulated into an approximation problem. Given the set of c o e f f i c i e n t s {a^ : i = 1, 2, . .., m + k + 1} the problem i s the determination of the a j _ ' s where m+1 ._. k S(x) = T a.x 1 L + I • a . j . - C x - f i ' . ) " . , L n l . L n x+m+1 l + 1=1 1=1 This system r e s u l t s i n m + k + 1 unknowns - exactly the number needed for a unique representation. Therefore a l l the c o e f f i c i e n t s must be computed. Despite i t s s i m p l i c i t y the mathematical representation of splines should never be used for computational purposes. Splines computed by t h i s method w i l l produce i l l - c o n d i t i o n e d systems of equations as m + k + 1 increases. I n t u i t i v e l y , a reason for t h i s can be seen from the example plotted i n Figure 1. Notice that the l a s t basis function 3 ( x - 6 ^ ) + i s zero nearly everywhere as w e l l as being extremely small r e l a t i v e to the other basis functions. Consequently i t i s possible to produce a l i n e a r combination of these basis functions which i s almost zero; that i s , the set of basis functions i s almost l i n e a r l y dependent. Using t h i s set w i l l produce a system of l i n e a r equations for the approx-imation problem whose corresponding matrix i s nearly singular. This matrix w i l l be i l l - c o n d i t i o n e d i n most cases. Hence i t i s better to choose basis functions which are more d i f f i c u l t to conceive a n a l y t i c a l l y but are more stable computationally. ELEMENTARY CUBIC SPLINE BASIS FUNCTIONS FIGURE 1 - 9 2.5 B-spline Representation In deriving a set of basis functions which produce w e l l ^ • conditioned systems; one would l i k e to s a t i s f y these two c r i t e r i a ; 1. That the support (region i n which the fu n c t i o n a l value i s nonszero) of the splines i s f i n i t e ; 2. That the number of knot i n t e r v a l s involved i n the support i s minimal (as small as p o s s i b l e ) . To t h i s end a set of basis functions i s derived on the concept of divided differences. D efinition^; The fu n c t i o n a l D defined by Df = f ( 6 ± F « 1 + 1 > 6 i + m + 1 ) • l + l i+m+1 l l+m i+m+1 i i s the divided d i f f e r e n c e of f_ of order m + 1 . The divided d i f f e r e n c e depends l i n e a r l y on f(x) . Also, and more important, the divided d i f f e r e n c e of order m + 1 i s zero for any polynomial of degree m . Now, for any function g, by the Lagrange i n t e r p o l a t i o n formula, - 10 -w h e r e co(x) = v^-^) * ( x - S i + ^ ) • • • ( x - 6 ^ ^ ^ ) . I n o r d e r t o d e f i n e B - s p l i n e s l e t t h e f u n c t i o n g(o"., ) b e t h e (m+ 1 )st d i v i d e d d i f f e r e n c e i + n . o f C ^ ^ ^ ^ x ) ™ . S u b s t i t u t i n g f o r g , t h e L a g r a n g e i n t e r p o l a t i o n f o r m u l a g i v e s m+1 C 6 . + n - x ) m g C 6 i ?6 i + 1,...,6 i + m + 1) = I T~ n=0 i + n w h e r e w ( x ) i s a s p r e v i o u s l y d e f i n e d . L e t t i n g i r a n g e o v e r -m t o k g i v e s e x a c t l y t h e n u m b e r o f B - s p l i n e s r e q u i r e d t o f o r m a s e t o f b a s i s f u n c t i o n s ; t h a t i s m + k + 1 f u n c t i o n s . T h u s , t h e s e t ( s ^ ( x ) : i = -m,-m + 1 , k } w h e r e m+1 (6 -x)? r \ V i + n + S . (x) = > , F A t — n=0 i + n w i t h w ( x ) = (x-6 .) • (x-6 . .. ) • • • (x-6 ) i s p r e c i s e l y t h e s e t o f l i + I XTTirrl B - s p l i n e b a s i s f u n c t i o n s . To c o m p l e t e t h e d e f i n i t i o n , t h i s s e t r e q u i r e s t h e a d d i t i o n o f 2m s u p p l e m e n t a r y k n o t s t o t h e o r i g i n a l k n o t s e t A = '^ 1 '"'"'^k+ 1 ^ " T h e s e k n o t s m u s t b e e x t e r n a l t o t h e o r i g i n a l k n o t s e t w i t h m k n o t s l e s s t h a n 6^ a n d m k n o t s g r e a t e r t h a n $^+1 * ^ n e P o s s l D l e m e t h o d o f c h o o s i n g t h e s e e x t r a k n o t s i s t h e f o l l o w i n g : Vol: 1 * - i = 6o kTT~ 6 = 6 + % + I ^ D I * 1 1 + k + 1 k + 1 k + 1 f o r i = 1 , 2 , ..., m - 11 -For B-splines i t can-Be shown that: 1. s . Cx) i s s t r i c t l y p o s i t i v e i n [6\. ,8.., ...]:• i J I x+m+1 ' 2. The support of s^(- x) i s f i n i t e and r e s t r i c t e d to the i n t e r v a l 1 ^ , 6 . ^ ] ; 3> Any s p l i n e S(x) can be uniquely represented as a l i n e a r combination of (s^(x) : 1 = *-m, -m + 1, ..., k} . I t i s simple to formulate an approximation problem from B-spline functions. Given the set of c o e f f i c i e n t s {a^ : i = -m k} the problem i s to f i n d the a ^ ' s i n k S(x) = T a.s.(x) x=-m where {s_^(x) : i = -m, -m + 1, k} i s the set of basis functions for B-splines. This system r e s u l t s i n m + k + 1 unknowns as i n the mathematical system. There are no redundant parameters and hence a l l the c o e f f i c i e n t s must be computed. The system of equations derived using B-splines remains wel'l-conditioned as m + k + 1 increases. In f a c t , one can produce numerical upper bounds on the condition number of the matrix of normal equations for a uniformly spaced knot set following the method described i n Schultz [18, pp. 70-72] . Also, because the basis functions give minimal support the systems produced are banded with the band-width dependent on the degree of the s p l i n e . - 12 -2.6 Example of B-spllne Representation To give a concrete example of what B-spline basis functions look l i k e , consider the p a r t i c u l a r case of cubic splines on a uniform p a r t i t i o n . An e x p l i c i t representation f o r the basis functions (s^(x) : i = -3, -2, k} can be developed i n the following manner: Given a uniform p a r t i t i o n , the mesh length i s h = T~r and therefore the i - t h knot i s 6 . = ^ k+1 i k+1 for i = -3, -2, k + 3, k + 4 . Thus s ± ( x ) n=0 »'<W for i = -3, -2,...., k with co'(x) = 4 n (x j=o Substituting f o r the knots gives i+n :fk+l - x) 3~ + n=0 C J ' i+n ^k+lj j (i+n - (k+l)x)^ n=0 (k+1)' i+n k+1 - 13 -Now 4 n j=0 (i+n-i-j) k + 1 4 n j=o k+1 Substituting back into s^Cx) 1 • 4 O + i j * ^ 0 C n _ j ) Jr^V s.(x) = 4 I n=0 (i+ri-(k+l)x)' (k + l ) 3 •(t'+l) 4 n (n-j) 3=0 = (k+1) 4 I n=0 (i+n-(k+1)x)' 4 n j=o (n-j) Expanding for n and j : s . ( x ) . ( k + 1 ) { C i - ( k + l ) x ) 3 _ ( i + l - ( k + l ) x ) 3 , Cl+2-(k+l)x) 3 _ (i+3-(k+l)x) 3 + i r — + — — - + , (i+4-(k+l)x) 3 , . + ^ - + i for i = -3, -2, ..., k . - 14 -Le t t i n g x' = (k+l)x - i -• 2 gives and defining S(x*) = s ^ k + D x - i - 2) r o S(x') 5 (k+l)< q + x ' ) 3 ( x * ) 3 q - x ' ) 3 (2-x') 3 6 4 ~ 6 24 ( x ' ) 3 d - x ' ) 3 , (2-x') 3 4 ~ 6 24 (1-x*) 3 , (2-x') 3 6 24 (2-x') 3 . 24 0 ' x' <_ -2 -2 <_ x' £ -1 -1 £ x' £ 0 0 <_ x* £ 1 1 £ x' £ 2 2 < x' . When the e x p l i c i t representation i s used on a. uniform p a r t i t i o n , of [0,1] with four i n t e r n a l knots; the set {s^(x) : i = -3, -2, ..., 4} i s as shown i n Figure 2. Equivalence of Spline Representations Although B-spliries provide a s u i t a b l e method f or solv i n g s p l i n e approximation problems the coef f i c i e n t s . o b t a i n e d are not extremely u s e f u l . In p a r t i c u l a r , . i t i s preferable to know the c o e f f i c i e n t s of the piece-wise continuous polynomials between adjacent knots than to know the CUBIC B-SPLINE BASIS FUNCTIONS FIGURE 2 - 16 -c o e f f i c i e n t s of B-spline basis functions. To t h i s end, i t i s appropriate to derive the unknown set of the piecewise continuous polynomial c o e f f i c i e n t s {c.. : i - 0,1, ..., k, j = 0, 1, ..., m}. from the known values of the B-spline c o e f f i c i e n t s {a^ • 1 " 1> 2, m + k +1} . Consider the functional value of the s p l i n e for each knot; that i s , at x = 6. . ' x Case I': For piecewise continuous polynomials c. . 1 dJ " dx 1 X J < dx J for i = 1, 2,..., k + 1 ; j =0, 1, m Case II: For B-splines -4 S(x) dx J d j v 6\ dx J £=-m * ^ x k d j I a„ — s„(x) £=-m dx 3 ^ 6. x for i - l , 2, Now —, s„(x) dx J k + 1; j = 0, 1, ..., m m 6. x o£ *f %+n - X ) + dxj riio U ' C W 6". X m+1. i ,J _ t 1 •- d"1 .m io-^W. d^  CV'X)+ Differentiating with respect to x ~ i C6-e+n - x )+ dx X to) ! C ^ - < - J 0 By equating the terms of the derivatives: i * £ for i = 1, 2, k + 1; j = 0/1, m - 18 -Section 3 THE LINEAR APPROXIMATION PROBLEM 3.1 Introduction Several methods are known for solv i n g l i n e a r approximation problems. These methods can be applied to problems involving general sets of basis functions. Spline approximation with a fi x e d knot set i s a p a r t i c u l a r a p p l i c a t i o n of the general problem/ de Boor and Rice [5] describe an approximation method invo l v i n g orthogonal projection. ' The basic idea i s to minimize the error ||y-u|| of approximating y by u by the orthogonal projection Py of y . Py i s best calculated using an orthonormal bas i s . Therefore, given a general set of basis functions for the approximation, an orthonormal set of basis functions must be derived, de Boor and Rice use a modified Gram-Schmidt process to generate such a set of basis functions. The most common technique used f o r . s o l v i n g l i n e a r approximation problems i s the method of normal equations (described i n Section 3.3). Patent [13] discusses the general l i n e a r l e a s t squares problem and l i n e a r l e a s t squares problem using splines i n d e t a i l giving r e s u l t s concerning the uniqueness of the s o l u t i o n and the symmetric and p o s i t i v e d e f i n i t e prop-e r t i e s of the associated l e a s t squares matrix. Patent also includes a program solving the s p l i n e approximtion problem with fixed knots. The basis functions used i n generating the system,of normal equations were B-splines. - 19 -Smith [19, pp. 110-119] develops the method of noraal equations for the pa r t i c u l a r case of spline approximation. The pa r t i c u l a r set of basis functions used are the mathematical representation. The.least squares matrix derived using these basis functions i s given i n f u l l . Golub [7] develops a method using Householder transformations Cdescr ibed i n Section 3.3). An Algol program based on this procedure for a general set of basis functions i s given i n Businger and Golub [3]. Although a suitable set of basis functions i s available which prevents i l l - c o n d i t i o n i n g as tbe number of knots increases, there i s s t i l l the problem of preventing i l l - c o n d i t i o n i n g as the knots become non-uniformly spaced. The orthogonal projection method counteracts this problem because the basis functions are orthonormalized before solution. The method of normal equations frequently produces i l l - c o n d i t i o n e d systems as shown i n the example cited i n Golub [7]. Solving the spline approx-imation problem using normal equations on a non-uniformly knot set i s a prime example of this i l l - c o n d i t i o n i n g . However, the method of Householder transformations counteracts t h i s problem because of the orthogonality of the transformations. Consequently, this method i s necessary for a stable solution to the spline approximation problem. 3.2 D e f i n i t i o n of the Linear Problem A l l aspects of the li n e a r problem are incorporated i n the following d e f i n i t i o n : - 20 -D e f i n i t i o n 5: Given a di s c r e t e set of data { (x£,y_£) : £ ~ 1, 2, . .. , n} and a function m S(x) = T a.s.(x) i = l where {s^(x) : i = 1, 2, ..., m} • i s a set of basis functions and {a^ : i = 1, 2, ..., m} i s a set of unknown c o e f f i c i e n t s occurring l i n e a r l y ; the l i n e a r  approximation problem i s to determine values for the a^'s to produce the "best f i t " of S(x) to the set of data. 3.3 Solution of the Linear Problem Most approximation problems consider the minimization of the lea s t squares error as s a t i s f y i n g the "best f i t " c r i t e r i a . The part-i c u l a r form of the le a s t squares error used i n t h i s case i s the square of the le a s t squares norm where: D e f i n i t i o n 6: Given a vector v = > the l e a s t squares norm of v, | v , i s Thus, given the function • S(x) and the data set {Cx^y^) .: £.= 1, 2, .. the s o l u t i o n of the l i n e a r problem becomes the. minimization of the l e a s t squares error by the appropriate choice of the unknowns ^a^ : i = 1, 2, That i s , by.finding n m 2 min I [yp - £ a s (x„)] : . {a ±} l«X- *- i - l 1 1 The usual method of s o l u t i o n , that of normal equations, i s developed i n the following way:. In order to f i n d the minimum n m 2 {a±} £=1 i = l 1 1 *-d i f f e r e n t i a t e the summation with respect to each of the parameters {a. : i = 1, 2, m} and set to zero. Thus l „ n m 9 VM: [y£- X a i s i ^ > ] } 9 a j i=i *- i = i n m = -2 I, [y» - I a s (x„).] • s (x») 1=1 L . i - l 1 1 * 3 *-= 0 for j = 1, 2, m . Rearranging the terms gives L=l ^ J l ' SJ ( X£ ) ] = J x ^ S J ( X £ } - 22 -Using the following changes: 1. Denote by. S the m x m matrix n S = S.. := for 1 = 1 , 2, m and j = 1, 2,..., m ; 2. Denote by y_ the m-dimension vector n f o r j = l , 2, ..., m ; and 3. Denote by &_ the m-dimension vector of unknown coef-i f i c i e n t s a_ = {a_^  : i = 1, 2, ..., m} ; the problem becomes one of f i n d i n g the s o l u t i o n to the system of equations Sa = y_ . Another method of s o l v i n g n m 2 min I [yp ~ I a . s . ( x o ) ] {a.} 1=1 *- i = l 1 1 X' l i s by using orthogonal transformations. This requires the following changes: 1. Denote by S the n x m matrix of function values; that i s , S = sl± = s i ( x £ } for £ = 1, 2, n and i = l , 2 m; 2. Denote by y_ the n-dimension vector of ordinates y_ = {y£ • & = 1? n) ; and 3. Denote by a. the m-dimension vector of unknown coef-ficients a, = (a^ : i - = 1, 2, m} ; then the problem becomes to find 2 min | |_y_ - SaJ | a. Consider multiplying the previous equation by an orthogonal T matrixs Q . Because multiplying by an orthogonal matrix does not change the norm; the linear least squares problem remains the same. Thus the problem becomes to find min | QTSa| | 2 . a_ T Now consider Q to be a series of orthogonal transformations T which transforms Q S into an upper triangular matrix R . If such a series can be found then the linear least squares problem reduces to finding min | | Q,Ty_ - Ra | | 2 . a. Since the zero part of R is independent of a,; i t is only necessary • T to. solve the system Ra = b_ where b. = (Q j ) f o r i = 1, 2, m . - 24 -The remainder of Q y_ contains the le a s t squares error; that i s , n IZ " S a l I (Q Ty)' i=m+l There exists a serie s of orthogonal transformations Q which w i l l reduce S to an upper t r i a n g u l a r matrix known as Householder trans-formations. They can be constructed as follows: Given a vector v construct a symmetric orthogonal matrix P such that Pv = w where w i s a unit vector.whose f i r s t element i s ±||v|| and whose remaining elements are zero. T T Householder showed.that for any two vectors v, w with v v = w w T there e x i s t s a symmetric orthogonal matrix P = I - 2uu . such that w = Pv . The symmetry and orthogonality of P i s proven i n Acton [1, P. 327] . The problem now i s to determine the required vector _u i n P . The method i s described i n Acton [1, pp. 324-329] with s l i g h t modifications and the derived u i s K v n "2 2 where K = 2 . v | ± 2v^| v| Two computational considerations come into e f f e c t when using Householder transformations. The f i r s t i s which sign to choose i n - 25 -computing _u . The choice i s to pick K such that K 2 = max{2| |v| | 2 + ^ 1 |v| | , 2 | | v | | 2 - 2v± | | v| | } v . / s • = max! K, K„ i n order to avoid c a n c e l l a t i o n . Thus i f v^ >^  0 choose K^; i f v^ < 0 choose IC, . The second consideration i s the computation of Pv . Rewriting Pv. = (I - 2uu T)v as T P v = Iv - 2 mi v T = v - _u(2_u v) T the s c a l a r 2u_ v i s computed f i r s t followed by the vector subtraction. Hence the matrix P need never be formed e x p l i c i t l y . This method i s far more e f f i c i e n t than forming P and performing a matrix m u l t i p l i c a t i o n . To manipulate Householder transformations to form the upper t r i a n g u l a r matrix consider applying P to the matrix S . This i s equivalent to applying. P to each column i n S . That i s , there i s a P^ such that P^S reduces column 1 of S, w^ , to w„ - 26 -where w| i s the transformed column 1 of S . Note that, a l l the other columns of S are alt e r e d by t h i s transformation. S i m i l a r l y , there i s a . such that. P^S reduces column 2 of S, w_2, to •±11*2 I where w^  i s , t h e transformed column 2 of S . However a l l other columns of S are alt e r e d also. In p a r t i c u l a r , column 1 reverts back to non-zero status which i s not desirable. Therefore, i n order to preserve the zeroes i n column 1 l e t P_ be the Householder transformation which reduces column 2 of S, w , ±11*2 I where w^  i s the transformed column 2 of S and w^  i s the f i r s t element of w' , This.transformation w i l l leave column 1 unchanged but "2 w i l l a l t e r a l l the remaining columns of S . Continuing the process m times;. S can be reduced to an upper t r i a n g u l a r matrix of the form - 27 -P x : % P.P-S = m 2 1 + n II ! ±11^11 • w2 0 0 w w„ w m m-1 w m ±||w | Thus i n the l i n e a r l e a s t squares problem using Q = p m **' P 2 P 1 and T T applying Q to S, Q reduces S to the upper t r i a n g u l a r matrix R . 3.4 Solution of the Fixed Knot Problem The f i x e d knot problem for splines i s a p a r t i c u l a r a p p l i c a t i o n of the general l i n e a r l e a s t squares problem. In the d e f i n i t i o n 5 for the l i n e a r problem allow the functions s^(x) to be the B-spline basis functions with a fixed knot set. This creates the f i x e d knot l e a s t squares approximation problem f o r s p l i n e s . This problem can be solved exactly l i k e the general problem using Householder transformations. However because a p a r t i c u l a r basis set i s being used a few computational considerations come into e f f e c t . The f i r s t of these i s that of basis function evaluation. Because the basis functions are evaluated often i t i s important that an e f f i c i e n t method of computation be a v a i l a b l e . In the f i x e d knot case the values of the basis function denominator vw1 (x) remain"'.constant throughout a l l computations s Therefore i t i s advantageous to c a l c u l a t e oo'(x) i n i t i a l l y and r e t a i n the values for future use. Also, because the basis functions have minimal support i t i s only necessary to compute a basis functional value i f i t i s within the range of support. Otherwise the functional value i s 0 . It i s possible to c a l c u l a t e the function at a l l points without t e s t i n g for the region of support. This has two disadvantages. F i r s t i t i s more time-consuming to c a l c u l a t e a value than to test f or the range. Second because of round-off error, the summation w i l l be the order of machine accuracy rather than zero. The second computational consideration i s the reduction of the r e s u l t i n g l e a s t squares matrix S . Because the basis functions have l i m i t e d support the matrix S has a banded structure. This banded structure i s not the familar handedness usually associated with matrices but rather a s p e c i a l structure dependent on the locationvof the abscissa point x^ . If the data point x^ i s outside the support region, that i s , x» £ [6., 6.. ... ] then the s». element of S w i l l be 0' . Hence the matrix w i l l have blocks of non-zero elements along the diagonal where each non^zero block of the matrix i s formed from the' values of X n where 6-. < x„ < 6., ,.. . This s p e c i a l banded structure can be Z i — Z — x+m+1 r - 29 -used i n the fixed knot approximation. Since the zeroes are unchanged by manipulation by Householder transformations i t i s only necessary to reduce the banded part of S . This causes a considerable reduction of computational e f f o r t . The t h i r d computational consideration i s the condition number of the matrix. Since S i s i n i t i a l l y a non-square matrix i t i s d i f -f i c u l t to state anything about i t s condition number. However the ortho-gonal transformations do not a f f e c t the norm (and hence the condition number) of S i n any way. Thus, by determining the condition number of the square upper t r i a n g u l a r part of the matrix R the condition of S has e f f e c t i v e l y been determined. The s t i p u l a t i o n for a well-conditioned matrix i s that the condition number be les s than in + k + 1 and that the condition number for a given problem should remain independent of the l o c a t i o n of the knots. - 30 -Section 4 THE VARIABLE KNOT APPROXIMATION PROBLEM 4.1 Introduction The f i x e d knot problem has been investigated thoroughly and adequate methods e x i s t for i t s s o l u t i o n . The v a r i a b l e knot problem involves choosing locations f o r the knots so as to provide the "best approximation" possible. Because the knots occur non-1inearly,. t h i s problem i s more complex. One a l t e r n a t i v e to sol v i n g the .variable knot problem i s to minimize the l e a s t squares error with respect to the knots. Another alternative, i s to p o s i t i o n the knots v i s u a l l y to a r r i v e at an approximation which i s reasonable to the eye although not ne c e s s a r i l y the "best" i n the le a s t squares sense. This requires graphical i n t e r a c t i o n with the sp l i n e approximation problem. By i n t e r a c t i o n a reasonable approximation can be derived manually, and good i n i t i a l guesses for an automatic minimization technique can be obtained. de Boor and Rice [6] solve the v a r i a b l e knot problem by an auto-matic technique. They minimize the l e a s t squares error i n i n t e g r a l form fb . , 2 \ . { [y - S ( x , A ) ] V ^ a over a l l splines of degree , m with k knots. The trapezoidal r u l e i s used to obtain an approximation to this, i n t e g r a l . A d i s c r e t e Newton's method i s applied to minimize each knot i n d i v i d u a l l y while the re s t of - 31 -the knots remain stationary. The knots are optimized by sweeping through the knot set from r i g h t to l e f t and re-evaluating the fi x e d knot problem each time. Smith [19, p. 110^119] presents the use of splines i n i n t e r -a c t i v e data f i t t i n g . Although he does not allow f or the p o s s i b i l i t y of res p e c i f y i n g c e r t a i n knot locations; he does allow the p o s s i b i l i t y of resp e c i f y i n g the e n t i r e knot set and attempting the fi x e d knot problem again. 4.2 D e f i n i t i o n of the Variable Knot Problem Whereas the l i n e a r l e a s t squares problem can be generalized to any set of basis functions; the d e f i n i t i o n of the non-linear approximation problem i s r e s t r i c t e d to s p l i n e functions and i s ref e r r e d to as the v a r i -able knot problem. D e f i n i t i o n 7: Given a d i s c r e t e set of data {(x^,y^) : I •= 1, 2 n} ; a set of knots i n s t r i c t l y increasing order A = : 1 = 1 , 2, ..., k> ; a set of s p l i n e basts functions: {s.(x) : i =» -m, -m .+ 1, k} > - 32 -and a set of l i n e a r c o e f f i c i e n t s {a.^ : i =" -m, l , ••... j k} ' k with S(x,A) = y a.s.Cx.A.) •1' « i x x - - m the v a r i a b l e knot problem i s to determine values f o r the set A to produce the "best f i t " of S to the set of data. Note that i n t h i s case the unknown parameters are the knots - not the l i n e a r c o e f f i c i e n t s {a^ : i = -m, -m + 1 , ..., k} . To be completely correct both the c o e f f i c i e n t s and the knots should be evaluated simultan-eously to produce the "best f i t " . This problem i s much more d i f f i c u l t . However, a f e a s i b l e a l t e r n a t i v e i s f o r each knot set A, to pick the c o e f f i c i e n t s {a^ : i = -m,_-m + 1, k} by sol v i n g the fi x e d knot problem. 4.3 Solution of the Variable Knot Problem In the introduction to t h i s section two a l t e r n a t i v e s were proposed f o r so l v i n g the v a r i a b l e knot problem: an automatic . procedure and an i n t e r a c t i v e procedure. I t i s best to discuss these i n the reverse order since the f i r s t follows inherently from the second. The use of i n t e r a c t i o n f o r sol v i n g s p l i n e approximation problems can be summarized i n the following steps: - 33 -1. Read i n the data set and graph i t on a graphics terminal. 2. Specify an i n i t i a l set of knots and overlay t h e i r l o c a t i o n on the x-axis. 3. Solve the fixed knot problem using t h i s knot set and over-lay the r e s u l t i n g s p l i n e curve. 4. Allow r e s p e c i f i c a t i o n of the l o c a t i o n of any of the knots. 5. Recalculate the fixed knot problem using the new knot set and graph the r e s u l t i n g spline approximation curve i f i t i s wanted. 6. Allow the options of r e s p e c i f y i n g knots, changing the number of knots, or optimizing the knot set. This i n t e r a c t i v e procedure produces a reasonable f i t much fa s t e r than an automatic procedure. For example, i f the o r i g i n a l knot set gives a poor approximation, the s i t u a t i o n can immediately be remedied by manipulating the knot.set d r a s t i c a l l y as opposed to the more cautious procedures of automatic techniques. This approach allows an extremely fast i n i t i a l approach to a good f i t . Then automatic refinement could use the r e s u l t i n g knot set to reach an optimal knot set quickly. There are several possible approaches for p o s i t i o n i n g the knots, de Boor and Rice [6, p. 12-18] present some of these for an automatic procedure but v a r i a t i o n s seem sui t a b l e for i n t e r a c t i v e placement. One possible approach suggested i s that a d d i t i o n a l knots be placed near the l o c a t i o n of the maximum error. This seems reasonable i n i t i a l l y as i t i s usually desired to produce a better f i t i n that area. Eventually, ihowever, the data i n that region w i l l become i n t e r -- 34 -polated which i s not the desired phenomenon. Also t h i s procedure does not compensate for the appropriate placement of fewer knots over the ent i r e range rather than the concentration of knots i n one place. Another p o s s i b i l i t y i s to place knots at positions of rapid change i n the data. This allows the polynomial to determine i t s own shape i n the i n t e r v a l nearly independently of the surrounding i n t e r v a l . This i s because at positions of rapid change, the highest order term of the piecewise polynomial dominates. I t i s p r e c i s e l y t h i s term that i s not included i n the. continuity constraints. This helps overcome one of the basic problems with .polynomials - that t h e i r o s c i l l a t o r y nature makes i t d i f f i c u l t for .them to adequately approximate data. With a graph of the data within reach i t i s possible to make reasonable predictions about.knot placement. This i s the greatest value of graphical i n t e r a c t i o n - the data and the a b i l i t y to manipulate the approximation are d i r e c t l y at hand. 4.4 Knot Optimization Despite the fact that reasonable approximations can be made to a set of data i n t e r a c t i v e l y , often a more formal f i t c r i t e r i a i s desired. This can be achieved by automatically r e f i n i n g the e x i s t i n g knot set l o c a l l y by minimizing the l e a s t squares error. Given the function S(x,A) and the data set {( x£> v£) : Z - 1,- 2, .. ., n} the s o l u t i o n of the va r i a b l e knot problem becomes the minimization of the le a s t squares error over the knot set A . That i s , f i n d n k min min £ ' [y„ - £ a s (x„,A)] A {a f} Z=l i=-m 1 1 L - 35 -I t i s not too c r u c i a l to obtain the global minimum i n the i n t e r a c t i v e case as a reasonable estimate of the knot set already e x i s t s . The purpose of optimization i s merely to r e f i n e t h i s estimate. The minimization method chosen i s known as COMPLEX which essen-t i a l l y involves! r e f l e c t i n g the function around i t s centroid. The d e t a i l s are not discussed here but are adequately described i n Box [2], Reasons for choosing COMPLEX involve the fact that i t does not require derivatives and that i t w i l l converge f a i r l y r a p i d l y towards the minimum. COMPLEX contains one a d d i t i o n a l feature. This i s that con-s t r a i n t s can be imposed on the function. These constraints can be e i t h e r e x p l i c i t - meaning that the independent v a r i a b l e can be bounded by some function or constant; or i m p l i c i t - meaning that the f u n c t i o n a l value can be bounded by some function or constant. In approximation using s p l i n e functions i t i s only necessary to have e x p l i c i t constraints to prevent the knots from coalescing. These constraints involve keeping the knots separated by a c e r t a i n distance. How t h i s distance i s determined i s p a r t i a l l y dependent on the machine pr e c i s i o n and hence the matrix used for s o l v i n g the fixed knot problem. Because of machine p r e c i s i o n the knots must be separated by a distance as l e a s t as great as the machine accuracy. Otherwise the l e a s t squares matrix w i l l be singular. More important, however, i s the condition number of the matrix used i n the f i x e d knot problem. As two knots converge towards each other, the two corresponding rows of the least.squares matrix become more l i n e a r l y dependent causing the condition number to r i s e . Therefore adequate con-- 36 -s t r a i n t s must be put on the knots to prevent them from causing numerical i n s t a b i l i t i e s . For these reasons, the following constraint was placed on each knot: set h = 6, ,- - 6_ and constrain each knot S. by k+1 0 l 6. . + .0001-h < 6. < - .0001-h . i - l — l — l + l for i = 1, 2, k . - 37 -NUMERICAL RESULTS 5.1 Introduction There are three main areas where l e a s t squares approximation can be used. These are: the approximation of mathematical functions; the approximation of experimental data; and computer-aided, design. In the f i r s t of these splines do not provide s u f f i c i e n t accuracy to warrant t h e i r use as fu n c t i o n a l approximations. In the. second arear.splines give good r e s u l t s . In the t h i r d area splines are able to f i t the contours of a design extremely well because of t h e i r piecewise nature. In order to demonstrate the p o s s i b i l i t i e s of i n t e r a c t i v e s p l i n e approximation three examples are given. The f i r s t example.fits cubic splines to an i n t e r e s t i n g set of data mainly to demonstrate the p o s s i b i l i t i e s of the method. The second example f i t s a cubic s p l i n e curve to a data set given i n de Boor and Rice [5], [6] involving data from a Titanium heat experiment. The f i n a l example i s the approximation of the o u t l i n e of a Volkswagen. This r e s u l t i s merely intended to demonstrate t h e . p o s s i b i l i t i e s of splines i n computer-aided design rather than having any p r a c t i c a l importance. A l l examples were run on an Adage Graphics Terminal connected to an IBM 360/67 duplex operating under MTS (Michigan Terminal System) located at the Un i v e r s i t y of B r i t i s h Columbia. Hardcopy plo t s were, obtained on a Calcomp p l o t t e r . Program l i s t i n g s and a user's guide are presented i n Appendices A and B r e s p e c t i v e l y . - 38 -The i n t e r a c t i v e system produces the following output: 1. Questions regarding what option the user would l i k e are printed on a conversational terminal. An example of one such session with the i n t e r a c t i v e system i s given i n Appendix B. 2. A graph of the data points, knots and f i t t e d curve are produced, on the graphics terminal. This i s i d e n t i c a l to the hardcopy that can be produced from i t as, for example, that of Figure 3. 3. A hardcopy pl o t of the graph can be produced upon.request. This p l o t i s i d e n t i f i e d by a t i t l e s p e c i f i e d as input.and a run number 'n' which indicates that i t i s the n-th hardcopy of the current terminal session. 4. A hardcopy printout as given i n Table I which corresponds to the p l o t . It can be matched to the plot by the t i t l e and run number. The hardcopy printout contains the following information: the abscissae and ordinates of the data points (the o r i g i n a l input to the system); the f i t t e d ordinates (the approximation to the ordinates by the system); the residuals (the diffe r e n c e between the ordinates and.the f i t t e d ordinates); the l e a s t squares error (the square root of the sum of the squares of the r e s i d u a l s ) ; the l o c a t i o n of the knots 6 , 6 , .. ., <S and the c o e f f i c i e n t s of the piecewise continuous polynomial ( s ^ x ) : i = 0, 1, ..., k} between each knot. That i s , between knot pair 0 ^ , 6 ^ ] the polynomial c o e f f i c i e n t i s the value f o r c.. i n m .(x) = I c (x - 6 ) j . A — A J -L j=0 5.2 Use of the System To demonstrate the use of the i n t e r a c t i v e system a simple set of 11 data points was chosen. As can be seen from the p l o t i n Figure 3 (ignoring the f i t t e d curve for the moment) the eye tends to approximate the data with the following curve: .9 .7 .6 .5 .4 .3 .2 .1 0 .1 .2 .3 .4 ,9 1, One would l i k e to manipulate splines so that they also approximate t h i s data with the above curve. - 40 -The f i r s t attempt was made to approximate the data with two uniformly spaced i n t e r n a l knots over '[0,1] . As can be seen from the p l o t i n Figure 3 and. l e a s t squares error of .1104, the r e s u l t wasn.not s u i t a b l e . A second attempt was made by moving the two i n t e r n a l knots further apart to .25 and .75 res p e c t i v e l y . This resulted i n the approximation given i n Figure 5, and Table I. The r e s u l t s are somewhat worse than the previous approximation (the l e a s t squares error was -.1574 as compared to .1104). Consequently t h i s approximation was eliminated and the previous uniformly spaced knot set retained. Further attempts were made to produce a better approximation by moving the two knots closer together. When i t became apparent that the l e a s t squares error had been reduced adequately with the knots located at 14 and .6 res p e c t i v e l y (plot and r e s u l t s not given), these values were given to the minimization procedure to f i n d the optimal knot loca t i o n s . The r e s u l t s are given i n Figure 4 and Table I I . The optimi-zation was terminated by the constraints on the knots. However, the lea s t squares,error was.reduced s i g n i f i c a n t l y (the f i n a l error was .0544). Presumably - t h i s i s the best approximation possible with two knots. The next step i n the procedure would then be to increase the number of knots to three. This was done and the i n i t i a l r e s u l t s of three uniformly spaced i n t e r n a l knots are shown in.Figure .5. and Table I I . The i n t e r e s t i n g thing to note i s that the le a s t squares error and approx-imation are i d e n t i c a l to that i n Table I. This i s because the t h i r d - 41 -knot located at the point of symmetry Is i n e r t and does not contribute to the approximation at a l l . Consequently the piecewise polynomial between .5 and .75 i s merely the polynomial between , .25 and .75 s h i f t e d . Since i t was useless to continue with three knots the number was increased to four. The i n i t i a l approximation with four uniform i n t e r n a l knots gave the r e s u l t i n Figure 6. This r e s u l t i s better than the optimized r e s u l t with 2 knots (least squares error of .0247 as compared .with .0544). But the curve i n the end regions, although smaller, contains more o s c i l l a t i o n s . The next step i s to vary the knots i n t e r a c t i v e l y . F i r s t , adjustment of the f i r s t , and l a s t knots towards the boundaries produced s i g n i f i c a n t l y better results, which terminated around .25 and .75 . Next, adjustment of the second and t h i r d knots produced be t t e r r e s u l t s continuously. The movement of the two knots was f i n a l l y terminated at .49999 and .50001 as i t was f e l t that they were coalescing too much (although the condition number of the l e a s t squares matrix was s t i l l only 5.96 and remaining f a i r l y constant). The r e s u l t s of t h i s approximation are spectacular. As can be seen from Table VI and Figure 7 the l e a s t squares error was .0000043 and the pl o t resembles the one expected. ( I t i s i n t e r e s t i n g to note the s i m i l a r i t y of the knot locations i n Figures 5 and 7. Although Figure 7 represents two nearly equal knots at the center the difference i n the:approximations obtained indicates that i t i s not n e c e s s a r i l y a good strategy to replace two coalescing knots by one knot and reducing the order of the system by one. 0 - 42 -This example demonstrates the power of the i n t e r a c t i v e system to approximate data. In p a r t i c u l a r , the r e s u l t s of the i n t e r a c t i v e procedure on the four knot approximation produced such impressive r e s u l t s that automatic refinement was unnecessary. FIGURE 3 - 44 -TEST DATA - 2 NON-ABSCISSAE 0.0 1.000000E-01 2.000000E-01 3.000000E-01 4.000000E-01 5.000000E-01 6.000000E-01 7.000000E-01 8.000000E-01 9.000000E-01 LOOOOOOE 00 THE LEAST SQUARES KNOT LOCATION 0 0.0 UNIFORM KNOTS ORDINATES FITTED OROINATES RESIDUALS 0.0 -7.871840E-0.0 3. 829566E-0.0 -4.988915E-0.0 -2. 5950 19E-L000000E-01 1. 877502E-5.000000E-01 5.000004E-9.000000E-01 8.122501E-1.000000E 00 L025949E LOOOOOOE 00 L049888E LOOOOOOE 00 9. 617022E-LOOOOOOE 00 L007872E ERROR IS L574225E-01 POLYNOMIAL POWER 03 7.8718U0E-03 02 -3.829566E-02 02 4.988915.E-02 02 2.595019E-02 01 -8.775020E-02 01 -4.619360E-07 01 8.774978E-02 00 -2. 594873E-02 00 -4.988807E-02 01 3. 829774E-02 00 -7.872522E-03 POLYNOMIAL COEFFICIENT 2,500000E-01 7.500000E-01 0 -7.871840E-03 1 L979496E 00 2 -1.940852E 01 3 4.230307E 01 0 -6.504732E-02 1 2.070669E-01 2 L231876E 01 3 -1.642503E 01 0 L065044E 00 1 2.070355E-01 2 -1.231879E 01 3 4.230357E 01 LOOOOOOE 00 TEST DATA - 2 NON-UNIFORM KNOTS Table I - -45 -TEST DATA - 2 KNOTS OPTIMIZED Oo 1. 2. 3. 4. 5. 6. 7. 8. 9. 1. THE L KNOT 0 ABSCISSAE 0 ORDINATES 0.0 0.0 0.0 0.0 1.000000E-01 5.000000E-01 9.000000E-01 000000E 00 •000O00E 00 000000E 00 .00O000E 00 FITTED ORDINATES RESIDUALS -7.379014E-03 7.379014E-03 000000E-01 000000E-01 OCOOOOE-01 000000E-01 000000E-01 O00O0OE-O1 000000E-01 000000E-01 00000OE-O1 000000E 00 EAST SQUARES ERROR IS 5. i,43563E-02 LOCATION POLYNOMIAL POWER 0.0 2.207708E-02 1.452189E-02 1. 511128E-02 1.223733E-01 4.999992E-01 8.776275E-01 015110E 00 014520E 00 779216E-01 007378E 00 •2. 207708E-02 1.452189E-02 1. 51 1 128E-02 •2.237328E-02 8.030709E-07 2.237248E-02 •1. 51 1062E-02 •1 .452056E-02 2.207839E-02 •7.378042E-03 POLYNOMIAL COEFFICIENT 4.999000E-01 5.001000E-01 0 -7.379014E-03 1 9.650481E-01 2 -8.405960E 00 3 1.701073E 01 0 4.994677E-01 1 5.313705E 00 2 1.710503E 01 3 -5.700724E 04 0 5.005385E-01 1 r- ~\ * '\ s r- "A .1 r\ ry 3.3 IOOCOI1 UU 2 -1.710493E 01 3 1.701070E 01 1.000000E 00 TEST DATA - 2 KNOTS OPTIMIZED Table II FIGURE 4 - 47 -TEST DATA - 3 UNIFORM KNOTS AESCISSAE ORDINATES FITTED ORDINATES RESIDUALS 0. 0 0 .0 -7. 872656E- 03 7. 872656E- 03 1. 000000E- 01 0 .0 3. 829633E- 02 -3. 829633E- 02 2. OOOOOOE-01 0 .0 -4. 988962E- 02 4. 988962E- 02 3. OOOOOOE-01 0 .0 -2. 595067E- 02 2. 595067E- 02 4. OOOOOOE-01 1 .OOOOOOE-01 1. 877483E- 01 -8. 774823E- 02 5. OOOOOOE-01 5 .OOOOOOE-01 4. 999984E- 01 1. 564622E- 06 6. OOOOOOE-01 9 oOOOOOOE-01 8. 122493E- 01 8. 775061E- 02 7. OOOOOOE-01 1 •OOOOOOE 00 1. 025949E 00 -2. 594928E- 02 8. OOOOOOE-01 1 .OOOOOOE 00 1. 049888E 00 -4. 988800E- 02 9. OOOOOOE-01 1 .OOOOOOE 00 9. 617013E- 0 1 3„ 829866E- 02 1- OOOOOOE 00 1 .OOOOOOE 00 1 . 007872E 00 -7. 871866E- 03 THE LEAST SQUARES ERROR IS 1.574225E-01 KNOT 0 LOCATION 0.0 2.500000E-01 5.000000E-01 7.500000E-01 P O L Y N O M I A L P O W E R 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 1.OOOOOOE 00 POLYNOMIAL COEFFICIENT -7.872656E-03 1.979532E 00 -1.940877E 01 4.230345E 01 -6.504667E-02 2.070408E-01 1.231880E 01 -1.642502E 01 4.999984E-01 3.286754E 00 5.340576E-05 -1.642529E 01 1.065044E 00 2.070377E-01 -1.231889E 01 4.230394E 01 TEST DATA - 3 UNIFORM KNOTS Table I I I FIGURE 5 FIGURE 6 - 50 -TEST DATA - 4 NON-UNIFORM KNOTS ABSCISSAE 0. 0 1. 000000E- 01 2. 000000E-01 3. OOOOOOE-01 4. OOOOOOE-01 5. O000OOE-01 6. OOOOOOE-01 7. 000000E-01 8. 000000E-01 9. 0 00000E-01 1. 0O0OOOE 00 THE LEAST SQUARES ORDINATES 0.0 0.0 0.0 0.0 1.000000E-01 5.000000E-01 9.000000E-01 1.000000E 00 1.000000E 00 1.OOOOOOE 00 1.000000E 00 ERROR IS 4 FITTED ORDINATES 3.96103 1E--1.974404E-2. 697052E-2.083834E-9.999895E-5.000019E-8. 999965E-1.OOOOOOE 9.99999 1E-9.999999E-9.999999E-266889E-06 KNOT 0 LOCATION 0.0 POLYNOMIAL POWER 2.500000E-01 4.999900E-01 5.000100E-01 7.500000E-01 RESIDUALS 08 -3.961031E-08 07 1.974404E-07 07 -2. 697052E-07 07 -2.083834E-07 02 1.032065E-06 01 -1.941880E-06 01 3.502471E-06 00 -4.593458E-07' 01 8.458737E-07 01 7.450581E-08 01 1.192093E-07 POLYNOMIAL COEFFICIENT 0 3.961031E-08 1 -1.450448E-02 2 2.175139E-01 3 -7.249289E-01 0 -1.357395E-03 1 -4.168792E-02 2 -3.261279E-01 3 3.405853E 01 0 4.999400E-01 1 6. 180718E 00 2 2.521675E 01 3 -8.394170E 05 0 5.000489E-01 1 6.180779E 00 2 -2.521672E 01 3 3.405827E 01 0 1.001356E 00 1 -4.164546E-02 2 3.259857E-01 3 -7.245331E-01 1.OOOOOOE 00 TEST DATA - 4 NON-UNIF0RM KNOTS' Table IV FIGURE 7 - 52 -5.3 Titanium Heat. Data de Boor and Rice [5], [6] presented a set of data obtained from a heat experiment i n v o l v i n g Titanium. This data was used to demonstrate the f i x e d knot program using a knot set with f i v e equally spaced knots [5, p. 18] and the v a r i a b l e knot program which involved the optimization of t h i s f i x e d knot set. Further studies involved the optimization of a non-uniform knot set containing f i v e knots. One of the best ways to test a system i s to t r y i t on previously done r e s u l t s and check the answers. I n i t i a l l y the i n t e r a c t i v e system was done with the uniform knot set s p e c i f i e d i n de Boor and Rice [5, p. 18]. The r e s u l t s and pl o t are given i n Table V and Figure 8. These r e s u l t s compare favorable with those i n de Boor and Rice (the r e s i d u a l was 1.16 as compared to t h e i r 1.24). However, as can be seen from the p l o t , the r e s u l t i n g approximation i s none too s a t i s f a c t o r y . It i s here that the power of the i n t e r a c t i v e system comes into e f f e c t . A f t e r moving the knots so that the pl o t produced i s more accurate, better r e s u l t s are obtained. The r e s u l t s given i n Table VI and Figure 9 indic a t e s u b s t a n t i a l improvement. The r e s u l t i n g knot set was used i n an i n i t i a l guesssfor knot optimization giving f i n a l r e s u l t shown i n Table VII and Figure 10. The reduction i n the l e a s t squares error was s u b s t a n t i a l (.092 as compared to 1.16) l a r g e l y because the i n t e r a c t i v e knot place-ment allowed for the deriving of accurate s t a r t i n g values. -.53 " TITANIUM HEAT DATA - 5 UNIFORM KNOTS 1 ABSCISSAE ORDINATES FITTED ORDINATES RESIDUALS 5. 950000E 02 6. 440000E- 01 6. 235025E- 01 2. 049747E-02 6. 050000E 02 6. 220000E- 01 6. 425120E-01 -2. 051 199E-02 6. 150000E 02 6. 3 800 00E-01 6. 516960E-01 -1. 369604E-02 6. 250000E 02 6. 490000E- 01 6. 534591E-0 1 -4. 459172E- 03 6. 350000E 02 6. 520000E- 01 6. 502057E- 01 1. 794243E-03 6. 450000E 02 6. 390000E-01 6„ 443403E-01 -5. 340301E- 03 6. 550000E 02 6. 460000E- 01 6. 382672E- 01 7. 732764E-03 6. 650000E 02 6. 570000E-01 6. 343911E- 0 1 2. 260887E- 02 6. 750000E 02 6. 520000E- 01 6. 351163E- 01 1. 688367E-02 6. 850000E 02 6. 550000E-01 6. 420718E-01 1. 292809E- 02 6. 950000E 02 6. 640000E- 01 6. 537848E- 01 1. 021518E- 02 7. 050000E 02 6. 630000E-01 6. 680065E-0 1 -5. 006507E- 03 7. 150000E 02 6. 630000E- 01 6. 824884E- 01 -1. 948840E-02 7. 250000E 02 6. 680000.E-01 6. 949821E- 01 -2. 698214E- 02 7. 350000E 02 6. 760000E- 01 7. 032389E-01 -2. 723893E-02 7. 450000E 02 6. 760000E-01 7. 050105E- 01 -2. 901048E- 02 7. 550000E 02 6. 860000E- 01 6. 980482E-01 -1. 204824E-02 7. 650000E 02 6. 790000E-01 6. 815894E- 01 -2. 589412E- 03 7. 750000E 02 6. 780000E- 01 6. 608155E- 01 1. 718443E-02 7. 850000E 02 6. 830000E-01 6. 4 2 39 3 3E-01 4. 060671E- 02 7, 950000E 02 6. 940000E- 01 6. 329899E-01 6. 101005E-02 8. 050000E 02 6. 990000E-01 6. 392726E-0 1 5. 972740E- 02 8. 150000E 02 7. 100000E- 01 6. 679080E-01 4. 209192E-02 8. 250000E 02 7. 300000E-01 7. 255636E- 01 40 436404E- 03 8o 350000E 02 7. 630000E- 01 8. 189063E-01 -5o 590630E-02 8. 450000E 02 8. 120000E-01 9. 505556E-01 -1. 385556E- 01 8. 550000E 02 9. 070000E- 01 1. 106939E 00 -1. 999393E-01 8.650000E 02 1. 0440O0E 00 1. 270436E 00 - 2 , 264375E- 01 8. 750000E 02 1. 336000E 00 1. 423430E 00 -8. 743048E-02 8. 850000E 02 1. 881000E 00 1. 548299E 00 3. 326998E- 01 8. 950000E 02 2. 169000E 00 1. 627424E 00 5. 415753E-01 9. 050000E 02 2. 075000E 00 1. 643184E 00 4. 318160E- 01 9. 150000E 02 1. 598000E 00 1. 583930E 00 1. 406908E-02 9. 250000E 02 1. 21 1000E 00 1. 461898E 00 -2. 508975E- 01 9. 350000E 02 9. 160000E- 01 1. 295290E 00 -3. 792901E- 01 9. 450000E 02 7. 460000E-01 1. 10231 1E 00 -3. 563114E- 01 9. 550000E 02 6. 720000E- 01 9. 011649E- 01 -2. 291650E- 01 9, 650000E 02 6. 270000E-01 7. 100576E-01 -8. 305758E- 02 9. 750000E 02 6„ 150000E- 01 5. 471910E-01 6. 780899E-02 9. 850000E 02 6. 070000E-01 4. 307705E-01 1. 762294E- 01 9, 950000E 02 6. 060000E- 01 3. 790006E-01 2, 2&9993E-01 1. 005000E 03 6. 090000E-01 4. 022449E-01 2. 067550E- 01 1. 015000E 03 6. 030000E- 01 4. 795058E- 01 1. 234941E-01 1. 025000E 03 6. 010000E-01 5. 819451E-0 1 1. 905493E- 02 1. 035000E 03 6. 030000E- 01 6. 807239E-01 -7. 772392E-02 1. 045000E 03 6. 010000E-01 7. 470052E-01 -1. 460052E- 01 1. 055000E 03 6. 1 10000E-01 7. 519506E-01 -1. 409506E-01 1. 065000E 03 6. 010000E-01 6. 667216E-01 -6o 572163E- 02 1. 075000E 03 6. 080000E- 01 4. 624799E- 01 . 1. 455200E- 01 THE LEAST SQUARES ERROR IS 1.157334E 00 KNOT LOCATION POLYNOMIAL POWER POLYNOMIAL COEFFICIENT 0 5.950000E 02 - 54 -6.750000E 02 7.550000E 02 8.350000E 02 9.050000E 02 9.950000E 02 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 6.235024E-01 2.472371E-03 -6.1 14945E-05 4.007443E-07 6.351163E-01 3.827438E-04 3.502920E-05 -3.747538E-07 6. 980482E-01 -1.207871E-03 -5.49116;E-05 1.1 11178E-06 8.189063E-01 1. 134087E-02 2.117710E-04 -2.936617E-06 1.643183E 00 -2.179519E-03 -4.049190E-04 3.034045E-06 3.790006E-01 -1.337662E-03 4.142732E-04 -4.806359E-06 1.075000E 03 TITANIUM HEAT DATA - 5 UNIFORM KNOTS Table V - 56 -5 NON-UNIFORM KNOTS TITANIUM HEAT DATA ABSCISSAE 5.950000E 02 6.050000E 02 6.150000E 02 6.250000E 02 6.350000E 02 6.450000E 02 6.550000E 02 6.650000E 02 6.750000E 02 6.850000E 02 6.950000E 02 7.050000E 02 7.150000E 02 7.250000E 02 7.350000E 02 7.450000E 02 7.550000E 02 7.650000E 02 7.750000E 02 7.850000E 02 7.950000E 02 8.050000E 02 8.150000E 02 8.250000E 02 8.350000E 02 8.450000E 02 8.550000E 02 8.650000E 02 8.750000E 02 8.850000E 02 8.950000E 02 9.050000E 02 9.150000E 02 9.250000E 02 9.350000E 02 9.450000E 02 9.550000E 02 9.650000E 02 9.750000E 02 9.850000E 02 9.950000E 02 1.005000E 03 1.015000E 03 1.025000E 03 1.035000E 03 1.045000E 03 1.055000E 03 1.065000E 03 1.075000E 03 THE LEAST SQUARES ERROR IS 1.142650E-01 KNOT LOCATION POLYNOMIAL POWER 0 5.950000E 02 ORDINATES FITTED ORDINATES RESIDUALS 6 .440000E -01 6 .252043E -01 1.879563E -02 6 .220000E -01 6 . 338547E -01 -1. 185482E -02 6 .380000E -01 6 .407921E -01 -2. 792177E--03 6 .490000E -01 6 .462172E -01 2.782739E -03 6 .520000E -01 6 .503309E -0 1 1. 669 100E-03 6 . 390000E -01 6 .533335E -01 -1.433358E -02 6 .460000E -01 6 .554263E -0 1 -9. 426299E--03 6 .570000E -01 6 . 568095E -01 1.904927E -04 6 . 520000E' -01 6 -576843E -01 -5. 684260E--03 6 .550000E--01 6 .582511E -01 -3.251 180E--03 6 .640000E--01 6 .587108E -01 5.2891 15E--03 6 . 630000E--01 6, . 592643E--01 3. 735680E-03 6 -630000E--01 6 .601120E--01 2. 888002E--03 6 . 680000E--01 6, . 614549E--01 6. 545052E--03 6 .760000E--01 6 .634936E--01 1.250640E--02 6, .760000E--01 6, » 664289E-01 9. 571020E--03 6 . 860000E--01 6, . 704616E--01 1.553836E--02 6. . 790000E--01 6, . 757923E-01 3.207672E--03 6. .780000E--01 6. .826219E--01 -4.621979E--03 6, • 830000E-01 6, . 9 115 11E--01 -8. 151092E--03 6, .940000E--01 7, .015807E--01 -7.580712E--03 6, .990000E--01 7. , 14 1112E--01 - 1. 511123E-02 7, , 100000E--01 7. .289435E--01 -1.894355E--02 7, . 300000E-01 7, • 462784E-01 -1. 627839E-•02 7, > 630000E-01 7» .6 6 3165E-- A 1 v/ 1 • ~J I u J *f I U . A -1 " vj J 8. . 120000E-•01 7o 908105E--01 2.118939E-•02 9. .070000E-•01 8, •572057E--01 4.979425E-•02 1. .044000E 00 1. 038639E 00 5.359702E-•03 1.336000E 00 1. • U02930E 00 -6.693059E-•02 1. 881000E 00 1. 859838E 00 2.116160E- 02 2. . 169000E 00 2. 161066E 00 7.933423E-•03 2. 075000E 00 2. 064366E 00 1.063279E- 02 1. 598000E 00 r. 624727E 00 -2.672801E-•02 1. 211000E 00 1. 185819E 00 2.518102E- 02 9. 160000E- 01 9. 078093E-•01 . 8. 190691E-•03 7. 460000E- 01 7. 544307E-•01 -8, 430697E- 03 6. 720000E- 01 6. 808758E-•01 -8.875843E- 03 6. 270000E-01 6. 432318E- 01 - 1. 623188E-02 6 o 150000E-01 6. 181573E- 01 -3.157331E- 03 6. 070000E-01 6. 028832E-•Q 1 4.116789E- 03 6- 060000E- 01 5. 955343E- 01 1.046569E- 02 6. 090000E-01 5. 942356E-0 1 1. 476442E-02 6. 030000E- 01 5. 971119E-01 5.888034E- 03 6, 0 10000E-01 6. 022885E-01 -1. 288563E- 03 6. 030000E- 01 6. 078902E-01 -4.890207E-03 6. 010000E-01 6. 120420E-01 - 1. 104198E-02 6. 110000E- 01 6. 128688E-01 -1.868859E-03 6. 010000E-01 6. 084959E-01 -7.495925E- 03 6. 080000E- 01 5. 970479E-01 1.095206E- 02 POLYNOMIAL COEFFICIENT I { - 57 -8.400000E 02 8.700000E 02 9.OOOOOOE 02 9.200000E 02 9.600000E 02 0 6.252043E-01 1 9.573922E-04 2 -9 .568968E-06 3 3.345709E-08 0 7.774121E-01 1 2.293389E-03 2 1.502197E-05 3 1.244829E-05 0 1.195837E 00 1 3.680510E-02 2 1o135369E-03 3 -4 .252815E-05 0 2.173558E 00 1 -9 .898979E-03 2 -2 .692169E-03 3 6.085630E-05 0 1.385562E 00 1 -4 .455817E-02 2 9.592120E-04 3 -7 .467897E-06 0 6.600300E-01 1 -3 .667146E-03 2 6.306361E-05 3 -3 .125027E-07 1.075000E 03 TITANIUM HEAT DATA - 5 NON-UNIFORM KNOTS Table VI - 58 -O - 59 -TITANIUM HEAT DATA - 5 KNOTS OPTIMIZED ABSCISSAE ORDINATES FITTED ORDINATES RESIDUALS 5. .950000E 02 6, .440000E--01 6. 225190E--0 1 2. 148103E--02 6. , 050000E 02 6, .220000E--01 6. .327376E--01 -1. 073763E-•02 6. , 150000E 02 6. , 380000E-01 6. 408629E--01 -2. 862922E-•03 6. , 250000E 02 6. .490000E--01 6. .471329E--01 1. 867094E-•03 6. .350000E 02 6. • 520000E-01 6. 517856E--01 2. 144140E-•04 6. .450000E 02 6. .390000E--01 6, .550590E--01 -1. 605903E-•02 6. 550000E 02 6. , 460000E-01 6. 571912E--01 -1. 1 191 16E-•02 6. . 650000E 02 6. , 570000E-01 6. 584202E--01 -1. 420241E-•03 6. 750000E 02 6. 520000E--01 6. 58984 1E--01 -6. 984055E- 03 6. 850000E 02 6. , 550000E-01 6. 591210E--01 -4. 121054E-•03 6. 950000E 02 6. 640000E-•01 6. 590686E-•01 4. 931286E- 03 7. 050000E 02 6. .630000E--01 6. 590654E-•01 3. 934622E- 03 7. 150000E 02 6. 630000E--01 6. 593490E-•0 1 3. 65094 1E-03 7. 250000E 02 6. 680000E-•01 6. 601578E-•01 7. 842168E-•03 7. 350000E 02 6. 760000E-•01 6. 617296E-•01 1. 427035E- 02 7. 450000E 02 6. 760000E-•01 6. 643026E-•01 1. 169730E-•02 7. 550000E 02 6. 860000E-•01 6. 681147E-•01 1. 788527E- 02 7o 650000E 02 6. 790000E-•01 6. 734042E-•01 5. 595822E-•03 7. 750000E 02 6. 780000E-•01 6. 804088E-•01 -2. 408907E- 03 7o 850000E 02 6. 830000E-•01 6. 893667E-•01 -6 = 366715E-•03 7. 950000E 02 6. 940000E-•01 7. 005159E-•01 -6. 515928E- 03 8. 050000E 02 6. 990000E-•01 7. 1409U6E-•01 -1. 509461E- 02 8. 150000E 02 7. 100000E-•01 7. 303407E-•01 -2. 034071E- 02 8. 250000E 02 7. 300000E-•01 7. 494920E-•01 -1. 949206E-02 8, 350000E 02 7 S 630000E-•0 1 7^  717870E-•0 1 -8. 786 9 95E-03 8. 450000E 02 8. 120000E-•01 7. 992001E-•01 1. 279991E- 02 8. 550000E 02 9. 070000E- 01 8. 665996E-•01 4. 040042E- 02 8. 650000E 02 1. 044000E 00 1. 038002E 00 5. 996864E-03 8. 750000E 02 1. 336000E 00 1. 378159E 00 -4. 215827E- 02 8. 850000E 02 1. 881000E 00 1 . 851600E 00 2. 939950E-02 8. 950000E 02 2. 169000E 00 2. 178109E 00 -9. 109914E- 03 9, 050000E 02 2. 075000E 00 2. 067245E 00 7. 754855E-03 9. 150000E 02 1. 598000E 00 1. 616635E 00 -1. 863577E- 02 9. 250000E 02 1. 21 1000E 00 1. 193341E 00 1. 765861E-02 9. 350000E 02 9. 160000E-01 9. 150418E-01 9. 581815E- 04 9. 450000E 02 7. 460000E- 01 7. 509134E- 01 -4. 913419E-03 9. 550000E 02 6. 720000E-01 6. 683230E-01 3. 676936E- 03 9. 650000E 02 6. 270000E- 01 6. 346373E-01 -7. 637367E-03 9. 750000E 02 6. 150000E-01 6. 188952E-0 1 -3. 895170E- 03 9. 850000E 02 6. 070000E- 01 6. 086131E- 01 -1. 613097E-03 9. 950000E 02 6. 060000E-01 6. 027005E-01 3. 299463E-03 1. 005000E 03 6. 090000E- 01 6. 002449E-01 8. 755121E-03 1. 015000E 03 6. 030000E-01 6. 00 3 335E-01 2. 666444E- 03 1. 025000E 03 6. 010000E- 01 6, 020537E- 01 -1. 053706E-03 1. 035000E 03 6. 030000E-01 6. 044927E-01 -1. 492694E-03 1. 045000E 03 6. 010000E- 01 6. 067377E-01 -5. 737711E-03 1. 055000E 03 6. 1 10000E-01 6. 078761E-01 3. 123831E- 03 1. 065000E 03 6. 010000E- 01 6. 069952E-01 -5. 995244E-03 1. 075000E 03 6. 0800O0E-01 6. 031824E-01 4. 817545E-03 THE LEAST SQUARES ERROR IS 9.286332E-02 KNOT LOCATION POLYNOMIAL POWER POLYNOMIAL COEFFICIENT 0 5.950000E 02 - 60 -8.395486E 02 8.733201E 02 8.989514E 02 9.179270E 02 9.681765E 02 0 6.225189E-01 1 1.134473E-03 2 -1.165708E-05 3 3.967317E-08 0 7.830166E-01 1 2.551379E-03 2 1 .744409E-05 3 1.084204E-05 0 1.306674E 00 1 4.082611E-02 2 1.115898E-03 3 -5.281300E-05 0 2.196890E 00 1 -6.059237E-03 2 -2.945114E-03 3 6.665658E-05 0 1.476898E 00 1 -4.5826 13E-02 2 8.494323E-04 3 -5.438899E-06 0 6.288878E-01 1 -1.658810E-03 2 2.952565E-05 3 -1. 521178E-07 1.075000E 03 TITANIUM HEAT DATA - 5 KNOTS OPTIMIZED Table VIII FIGURE 10 - 62 -5.4 The Bug The use of s p l i n e approximation i n computer-aided design i s a r e l a t i v e l y unexplored area. The Volkswagen data presented i n t h i s section demonstrates possible applications of splines i n the area. Also the approximation demonstrates the use of a large knot set. The. i n i t i a l approximation of 18 equally spaced knots produced a Volkswagen with a dented hood as can be seen i n Figure 11. Inter-active manipulation of the knots produced the more reasonable resemblance to a Volkswagen shown i n Figure 12 and Table VIII.. Since the desired r e s u l t s were based on representing the data accurately by sight rather than minimizing the le a s t squares error the knots were not optimized. - 64 -THE BUG ABSCISSAE ORDINATES 3. 500000E-•01 3 .OOOOOOE-•01 3. 525000E- 01 4 .OOOOOOE-•01 3. 550000E-•01 5.500000E-•01 4. 300000E-•01 6 . 500000E-•01 4. 700000E-•01 8 .OOOOOOE--01 5. 000000E-•01 9 .500000E-•01 5. 500000E-•01 1 .060000E 00 6. OOOOOOE-•01 1 .200000E 00 7. OOOOOOE-•01 1 .250000E 00 8. OOOOOOE-•01 1 •290000E 00 9. OOOOOOE-•01 1 .3.00000E 00 1. OOOOOOE 00 1 „430000E 00 1. 100000E 00 1 .500000E 00 1. 200000E 00 1 .560000E 00 1, 300000E 00 1 .620000E 00 1. 400000E 00 1 .675000E 00 1. 500000E 00 1 •720000E 00 1. 600000E 00 1 .750000E 00 1 . 700000E 00 1 .780000E 00 1. 800000E 00 1 .810000E 00 1. 900000E 00 1 .840000E 00 2. OOOOOOE 00 1 .860000E 00 2. 100000E 00 1 „890000E 00 2. 200000E 00 1 .915000E 00 2. 300000E 00 1 .930000E 00 2, 4U0O0OE 00 1 .950000E 00 2. 500000E 00 1 .960000E 00 2. 600000E 00 1 .970000E 00 2. 700000E 00 1 .980000E 00 2. 800000E 00 1 .990000E 00 2. 9C0000E 00 2 .OOOOOOE 00 3. OOOOOOE 00 2 .005000E 00 3. 100000E 00 2 .005000E 00 3. 200000E 00 2 .005000E 00 3. 350000E 00 2 .005000E 00 3. 400000E 00 2 . 100000E 00 3. 500000E 00 2.250000E -00 3. 600000E 00 2 .430000E 00 3. 700000E 00 2 .600000E 00 3. 800000E 00 2 .750000E 00 3. 900000E 00 2 .890000E 00 4. OOOOOOE 00 2 .920000E 00 4. 100000E 00 2 .950000E 00 4. 2C0000E 00 2 .970000E 00 4. 300000E 00 3 .OOOOOOE 00 4. 400000E 00 3 .020000E 00 4. 500000E 00 3 .040000E 00 4. 6O0000E 00 3 .050000E 00 4. 700000E 00 3 .055000E 00 4. 800000E 00 3 .065000E 00 4. 900000E 00 3 .075000E 00 5. OOOOOOE 00 ' 3 .080000E 00 ORDINATES RESIDUALS 4 .035096E-•01 -1. 035095E-01 4.113406E-•0 1 -1. 134065E-02 4 .192352E-•01 1. 307648E-•01 6 .758505E-•01 -2. 585047E-02 8 .169750E-•01 -1. 697499E-02 9 .183491E-•01 3. 165080E-02 1 .068680E 00 -8. 680243E-03 1 . 18203 1E 00 1. 796837E-02 1 .275474E 00 -2. 547405E-02 1 .276578E 00 1. 342188E-02 1 .304118E- 00 -4. 118513E-03 1 .421022E 00 8. 977752E-03 1 .509212E 00 -9. 212483E-03 1 .564281E 00 -4. 281554E-03 1 .616894E 00 3. 105875E-03 1 .668297E 00 6. 702900E-03 1 .713498E 00 6. 501570E-03 1 .752367E 00 -2. 367944E-03 1 .785924E 00 -5. 924519E-03 1 •815188E 00 -5. 188584E-03 1 .841183E 00 -1, 182909E-03 1 .864781E 00 -4. 782557E-03 1 .887129E 00 2. 871351E-03 1 .908055E 00 6. 944049E-03 1 .927378E 00 2. 621699E-03 1 .944908E 00 5. 091667E-•03 1 .960458E 00 -4. 580617E-04 1 .973843E 00 -3. 842760E-03 1 .984876E 00 -4. 87637 1E-03 1 .993369E 00 -3. 369596E-03 1 .999138E 00 8. 617891E-04 2 .002012E 00 2. 987819E-03 2 .002145E 00 2. 855293E-03 2 .001301E 00 3. 698572E-03 2 .03421 1E 00 -2. 921 108E-02 2 .079496E 00 2. 050392E-02 2 .238708E 00 1. 129068E-02 2 .428134E 00 1. 865786E-03 2 .612069E 00 -1. 206875E-02 2 .767555E 00 -1. 755604E-02 2 .871636E 00 1. 836338E-02 2 .916960E 00 3. 039550E-03 2 .946946E 00. 3. 053293E-03 2 .973553E 00 -3. 552493E-03 2 .996966E 00 3. 032722E-03 3 .017381E 00 2. 619695E-03 3 .034980E 00 5. 019724E-03 3 .049955E 00 4. 446507E-05 3 .062496E 00 -7. 496823E-03 3 .072790E 00 -7. 791314E-03 3.081028E 00 -6. 028723E-03 3 .087397E 00 -7. 396754E-03 - 65 -5. 100000E 00 3. 090000E 00 3. 092069E 00 -2. 069191E- 03 5. 200000E 00 3. 100000E 00 3. 095146E 00 4. 853774E-03 5. 300000E 00 3. 10000OE 00 3. 096716E 00 3. 283732E- 03 5. 400000E 00 3. 100000E 00 3. 096868E 00 3. 132608E-03 5. 500000E 00 3. 100000E 00 3. 095680E 00 4. 319567E- 03 5. 600000E 00 3o 100000E 00 3. 093245E 00 6. 755099E-03 5. 700000E 00 3. 100000E 00 3. 089649E 00 1. 035092E- 02 5. 800000E 00 3. 090000E 00 3. 084974E 00 5. 025774E-03 5. 900000E 00 3. 080000E 00 3. 079309E 00 6. 905058E- 04 6. OOOOOOE 00 3. 070000E 00 3. 072740E 00 -2. 739966E-03 6. 100000E 00 3. 060000E 00 3. 065260E 00 -5. 259581E- 03 6. 200000E 00 3. 050000E 00 3. 056500E 00 -6. 500497E-03 6. 300000E 00 3. 040000E 00 3. 046000E 00 -5. 999859E- 03 6. 400000E 00 3. 030000E 00 3. 033295E 00 -3. 295366E- 03 6. 500000E 00 3. 010000E 00 3. 017928E 00 -7. 928025E- 03 6, 600000E 00 3. OOOOOOE 00 2. 999435E 00 ' 5. 641840E- 04 6. 700000E 00 2. 980000E 00 2. 977354E 00 2. 645336E- 03 6. 800000E 00 2. 950000E 00 2. 951224E 00 -1. 224987E- 03 6. 900000E 00 2. 920000E 00 2. 920587E 00 -5. 867698E- 04 7. OOOOOOE 00 2. 890000E 00 2. 884996E 00 5. 003192E-•03 7. 100000E 00 2. 850000E 00 2. 844051E 00 5. 948193E- 03 7. 200000E 00 2. 800000E 00 2. 797825E 00 2. 174579E-•03 7, 300000E 00 2. 750000E 00 2. 746466E 00 3. 533684E- 03 7. 400000E 00 2. 690000E 00 2. 690125E 00 -1. 264103E-•04 7, 500000E 00 2. 630000E 00 2. 628954E 00 1. 045644E- 03 7. 600000E 00 2. 560000E 00 2. 563101E 00 -3. 100801E-•03 7. 700000E 00 2. 490000E 00 2. 492720E 00 -2„ 719902E- 03 "7 / « Q A A A A A T? VJ O W V W SJ .LJ 00 O L « fi T n n n n w ~T t- W \J \J \J U A A £. o / i n n c c u nr\ * r\ it i n n c o_ \J T - T 4. \J _ / i_ -A O vy ~J 7. 900000E 00 2. 340000E 00 2. 338963E 00 1. 036749E- 03 8, O OOOOE 00 2. 250000E 00 2. 255831E 00 -5. 881384E-•03 8. 100000E 00 2. 160000E 00 2. 169028E 00 -9. 028815E-•03 8.200000E 00 2 © 080000E 00 2. 079279E 00 7. 202886E-•04 8. 300000E 00 2. OOOOOOE 00 1, 987667E 00 1. 233252E-•02 8. 400000E 00 1. 900000E 00 1. 895224E 00 4. 775889E-•03 So 500000E 00 1. 805000E 00 1. 802980E 00 2. 018948E-•03 8.600000E 00 1. 7050C0E 00 1. 710693E 00 -5. 694207E-•03 8. 700000E 00 1. 605000E 00 1. 613008E 00 -8. 008420E-•03 8. 800000E 00 1. 500000E 00 1. 503284E 00 -3. 283981E--03 8.900000E 00 1. 375000E 00 1. 374889E 00 1. 102020E-•04 9. OOOOOOE 00 1. 240000E 00 1. 221188E 00 1. 881186E--02 9. 100000E 00 1. 030000E 00 1. 038939E 00 -8. 940164E-•03 9. 200000E 00 8. 300000E--01 8. 384878E--01 -8. 487869E--03 9. 300000E 00 6. 400000E--01 6. 335645E--01 6. 435510E-•03 9, 400000E 00 4. 450000E--01 4. 379095E--01 7, 090468E--03 9. 500000E 00 2. 500000E--01 2. 652552E--01 -1. 525517E--02 9. 550000E 00 2. 000000E--01 1. .918465E--01 8. , 153468E-03 THE LEAST SQUARES ERROR IS 1.897547E-01 KNOT LOCATION POLYNOMIAL POWER POLYNOMIAL COEFFICIENT 0 3.500000E-01 0 4.035096E-01 1 3.119434E 00 2 5.245779E 00 3 -2.106841E 01 1 6.000000E-01 - 66 -2 9.000000E-01 3 1.OOOOOOE OO 4 1.2000OOE OO 5 1.400000E OO 6 2.OOOOOOE OO 7 3.OOOOOOE OO 8 3.200000E OO 9 3.400000E OO 10 3.500000E 00 11 3.900000E 00 0 1.182031E 00 1 1.792031E 00 2 -1.055562E 01 3 1.979568E 01 0 1.304118E 00 1 8.034888E-01 2 7.260439E 00 3 -3.605061E 01 0 1.421021E 00 1 1.174047E 00 2 -3.554706E 00 3 6.331032E 00 0 1.564281E 00 1 5.119148E-01 2 2.439298E-01 3 -1.019521E 00 0 1.668296E 00 1 4.870993E-01 2 -3.677540E-01 3 1.702099E-01 0 1.864781E 00 1 2.299367E-01 -6. 161 657E-02 3 -3.110615E-02 0 2.002012E 00 1 1.324511E-02 2 -1.546146E-01 3 3.537191E-01 0 2.001300E 00 1 -6.130829E-03 2 5.760336E-02 3 9.639282E 00 0 2.079494E 00 1 1.173615E 00 2 5.841224E 00 3 -1.655806E 01 0 2.238708E 00 1 1.845140E 00 2 8.737149E-01 3 -3.827115E 00 0 2.871635E 00 1 7.071088E-01 2 -3.718792E 00 3 1.180118E 01 - 67 -12 13 14 15 16 17 18 19 4.OOOOOOE OO 5.OOOOOOE OO 6.OOOOOOE OO 7.OOOOOOE OO 8.OOOOOOE OO 8.500000E OO 9.OOOOOOE OO 9.550000E OO 0 2.916960E OO 1 3. 173848E-01 2 -1.784239E-01 3 3.148174E-02 0 3.087396E OO 1 5.497074E-02 2 -8.398247E-02 3 1.435216E-02 0 3.072739E OO 1 -6.993294E-02 2 -4.092169E-02 3 -7.690954E-02 0 2.884996E OO 1 -3.825288E-01 2 -2.716408E-01 3 2.506399E-02 0 2.255880E OO 1 -8.505980E-01 2 -1.964577E-01 3 1.721067E-01 0 1.802979E 00 1 -9. 179758E-0 I 2 6. 170338E-02 3 -1.105849E 00 0 1.221187E 00 1 -1.685660E 00 2 -1.597072E 00 3 2.289328E 00 THE BUG Table VIII - 69 -Section 6 SUMMARY AND FUTURE POSSIBILITIES The r e s u l t s shown i n the previous section provide a b r i e f out-look into the p o s s i b i l i t i e s of -interactive l e a s t squares approximation using s p l i n e functions. There remain, however, many problems and pos-s i b i l i t i e s for further research i n t h i s area. With regards to the basis functions; a p o s s i b i l i t y i s to f i n d more appropriate methods of choosing the external knots. Schultz [18, p. 73] suggests choosing the external knots so that 6.., =.6. l + l l - 6. l + l l For the case of a uniform knot set t h i s i s equivalent, to the method described i n Section 2. However, for the non-uniform knot set the knot spacing i s d i f f e r e n t . Schultz's choice i s poor when the boundary knot and f i r s t i n t e r i o r knot are coalescing as i t causes high peaks i n the basis functions i n v o l v i n g the external and boundary knots. This causes a sharp increase i n the.condition number of the l e a s t squares matrix - a s i t u a t i o n which should be avoided. The second p o s s i b i l i t y (interconnected with the f i r s t ) i s to f i n d optimal weights.for the b a s i s functions i n v o l v i n g the external knots so as to reduce the condition number. At present the f i r s t and l a s t basis function are a r b i t r a r i l y m u l t i p l i e d by (m + 1) so as to obtain an adequate portion of these basis functions when d e r i v i n g the l e a s t squares matrix. An i n t e r e s t i n g sutdy could be made of the 6 1 - <S0 -m < i < -1 ; 6, ,, - 6, k + l < i < k + m k+1 k — — - 70c-t h e o r e t i c a l and p r a c t i c a l importance of such weights for producing good bounds for the norm of the l e a s t squares matrix. A t h i r d p o s s i b i l i t y i s to produce a more e f f i c i e n t method of computing the B-spline basis functions using recursive formulae. These methods are described i n Lyche and Schumaker [TO]. They provide sub-s t a n t i a l reductions i n the time involved i n evaluating the B-splines as i d e n t i c a l terms are not recomputed. Since the B-spline basis function evaluation uses the most s i g n i f i c a n t amount of time i n the system; de r i v a t i o n of an improved, more e f f i c i e n t algorithm f or spl i n e basis function evaluation would be highly j u s t i f i e d . With regards to the f i x e d knot l e a s t squares problem; a better estimate of the l e a s t squares norm might aid i n the approximation process. Both de Booraand Rice and Patent use an i n t e g r a l norm i n t h e i r approximation as opposed to the d i s c r e t i z e d norm (sum of the squares of the residuals) used i n t h i s system, de Boor and Rice [5], [6] estimate the norm by the trapezoidal r u l e ; hence, i n e f f e c t , obtaining a l i n e a r i n t e r p o l a t i o n between f i t t e d data points. Patent [13] used F i l o n quadrature (an in t e r p o l a r y quadrature), the theory of which i s discussed i n his the s i s . Since the main i n t e r e s t of t h i s system was i n t e r a c t i v e approximation rather than l e a s t squares approximation, the d i s c r e t i z e d norm was adequate. With regards to the v a r i a b l e knot problemj better algorithms are needed for the non-linear l e a s t squares problem. In p a r t i c u l a r , an algorithm for the s p e c i f i c case of knot optimization would be u s e f u l . - 71 -Also the p o s s i b i l i t y of keeping knots stationary while optimizing other knots could be u s e f u l i n the i n t e r a c t i v e case. Interactiveetechniques are a v i t a l asset to numerical compu-ta t i o n . For too long, numerical analysts have been developing 'black box* algorithms i n which numbers are read i n and (hopefully) correct answers are printed out. By i n t e r a c t i n g with the procedure, what i s happening to the s o l u t i o n can be observed immediately. Interaction gives the power to interrupt the s o l u t i o n process and change the s t a r t i n g conditions of the problem. Part of the future of numerical computation depends on the development of i n t e r a c t i v e methods. Hopefully, the use of an i n t e r -a c t i v e approach rather than a 'black box' approach w i l l be valuable i n solving many numerical problems. BIBLIOGRAPHY Acton, F. S., Numerical Methods that (Usually) Work, New York: Harper and Row, 1970. Box, M. J . , "A New Method of Constrained Optimization and a Comparison with Other Methods". The Computer Journal 8. 1965, pp. 42-52. Businger, P. and Golub, G., "Linear Least Squares Solutions by Householder Transformations". Numerische Mathematik 7, 1965, pp. 269-276. Carasso, C. and Laurent, P. J . , "On the Numerical Construction and the P r a c t i c a l Use of Interpolating Spline Functions". Proceedings of the IFIP Congress, 1968, pp. A6-A9. de Boor C. and Rice, J. R., "Least Squares Cubic Spline Approx-imation I - Fixed Knots". Purdue Un i v e r s i t y : Technical Report CSD TR 20, 1968. de Boor, C. and Rice, J. R., "Least Squares Cubic Spline Approx-imation II - Variable Knots". Purdue U n i v e r s i t y : Technical Report CSD TR 21, 1968.. Golub, G., "Numerical Methods for Solving Linear Least Squares Problems". Numerische Mathematik .7,, 1965, pp. 206-216. G r e v i l l e , T. N. E., "Introduction to Spline Functions". Theory . and A p p l i c a t i o n of Spline Functions, 1969, pp. 1-35. G r e v i l l e , T. N. E. (ed.), Theory and A p p l i c a t i o n of Spline Functions, New York: Academic Press, 1969. Lyche, T. and Schumaker, L., "Computation of Smoothing and Inter-polating Natural Splines v i a Local Bases". U n i v e r s i t y of Texas Center for Numerical Analysis, 1971. - 73 ~ [11] O'Reilly, D., ' "Introduction to Non-linear Function Optimization". U n i v e r s i t y of B r i t i s h Columbia Computing Centre, 1971. [12] O^Reilly, D., 1 "Non-linear Function Optimization - The Complex Method of M. J . Box". U n i v e r s i t y of B r i t i s h Columbia Computing Centre, 1971. [13] Patent, P. D., "Least Square Polynomial Spline Approximation". C a l i f o r n i a I n s t i t u t e of Technology: t Ph.D. D i s s e r t i o n , 1972. [14] Rice, J...R. , The Approximation of Functions - Vol. I I , Reading, Massachusetts; Addison-Wesley, 1969. [15] Richardson, J . A. and Kuester, J. L., "The Complex Method for Constrained Optimization". Comm. A.CM. 8, 1973, pp. 487-489. [16] Schoenberg, I. J . ( e d l ) , Approximations with Special Emphasis  on Spline Functions, New York: Academic Press, 1969. [17] Schultz, M., " M u l t i v a r i a t e Spline Functions and E l l i p t i c Problems". Approximations with Special Emphasis on Spline  Functions, 1969, pp. 279-347. [18] Schultz, M., Spline Analysis, New York: P r e n t i c e - H a l l , 1973. [19] Smith, L., The Use of Man-Machine Interaction i n D a t a - F i t t i n g Problems". Stanford U n i v e r s i t y : Technical Report CS 131, 1969. - 74 -Appendix A . PROGRAM LISTINGS The following pages contain l i s t i n g s of the majority of the program used i n the i n t e r a c t i v e s p l i n e approximation. The subroutines consist of a main co n t r o l routine ..which c a l l s the various subroutines i n the proper order. As i s the. problem with most programs written f o r s p e c i f i c hardware configurations most of the input/output subroutines are machine dependent. This program i s no exception. However the subroutines should not be d i f f i c u l t to change to a system c o n s i s t i n g of correspond-ing hardware; that i s , a conversational terminal, a graphics terminal and a p l o t t e r . The routines i n question are: 1. INFREE - which i s a free format input routine. A l l c a l l s to INFREE can be replaced by FORTRAN READ and FORMAT statements. 2. The plo t subroutines PLOT, SCALE, AXIS, SYMBOL and PLOTND need only minor modifications from one plot system to another. Note, however, that PSEND i s a s p e c i f i c subroutine f o r the graphics terminal which simply p l o t s a l l currently a v a i l a b l e p l o t information. 3. AGPLOT creates a hardcopy of a plo t currently being d i s -played on a graphics terminal. I f such a routine i s unavailable t h i s routine can be omitted. However t h i s causes the loss of the a b i l i t y to keep the plo t information permanently. - 75 " There i s also one matrix manipulation subroutine which, i s not given i n the program. The subroutine INVERT i s used to obtain the condition number of the matrix. However a s i m i l a r subroutine i s a v a i l -able almost everywhere. The code for the function minimization i s also not given as i t was not written by the author. However a s i m i l a r subroutine i s av a i l a b l e i n Richardson and Kuester [15]. c* C* C O N T R O L R O U T I N E C * C* DIMENSION X(200), Y(200), YF (200) , RES (200) DIMENSION A (200), DELTA (50) INTEGER YES /'YES'/, NO /'NO'/, ANSWER /» '/ C C *** ENTER AND PLOT THE DATA POINTS 207 CALL DATAIN (X, Y, N) 205 CALL DATAPL (X, Y, XMIN, XMN, XMAX, EX, YMIN, EY, N) C C *** ENTER AND PLOT THE ORIGINAL KNCT SET 206 ERRORP = 1000. CALL KNOTIN (X, DELTA, M, N, MK, S206) CALL KNOTPL (DELTA, XMIN, DX, YMIN, DY, M, HK) C C *** COMPUTE THE INITIAL SPLINE APPROXIMATION CALL FIXED (X, Y, YF, RES, A, DELTA, ERROR, ERRCBF, M, N, MK) C c *** PLOT THE SPLINE APPROXIMATION CALL SPLNPL (A, DELTA, XMIN, XMN, XMAX, DX, YMIN, DY, M, MK) C C *** ASK IF HARDCOPY OUTPUT IS WANTED CALL SPLOUT (X, Y, YF, RES, A, DELTA, ERROR, N, M, MK) C *** ASK WHETHER OR NOT TO CONTINUE ITERATING 202 WRITE (6,600) 600 FORMAT ('SDO YOU WISH TO VARY THE KNOTS') CALL INFREE (2058, ANSWER, 4) IF (ANSWER.EQ.YES) GO TO 200 IF (ANSWER.EQ.NO) GO TO 201 WRITE (6,601) 601 FORMAT (* SPLEASE ANSWER YES OR NO') GO TO 202 C C *** VARY THE KNOTS 200 CALL CONTRL (X, Y, YF, RES, A, DELTA, ERROR, N, M, MK) GO TO 202 C C *** ASK WHETHER OR NOT TO OPTIMIZE THE KNOT SET 201 WRITE (6,602) 602 FORMAT ('SDO YOU WISH TO OPTIMIZE THE KNOT SET') CALL INFREE (2058, ANSWER, 4) IF (ANSWER.EQ.YES) GO TO 203 IF (ANSWER.EQ.NO) GO TO 204 WRITE (6,601) GO TO 201 C C *** OPTIMIZE THE KNOT SET 203 CALL OPT (X, Y, YF, RES, A, DELTA, ERROR, N, H , H K ) C - 77 -C *** ASK WHETHER OR NOT TO RESTART WITH A NEW KNOT SET 204 WRITE (6,603) 603 FORMAT ('&DO YOU WISH TO RESTART WITH A NEW KNOT SET') CALL INFREE (2058, ANSWER, 4) IF (ANSWER.EQ.YES) GO TO 205 IF (ANSWER.EQ.NO) GO TO 209 WRITE (6,601) GO TO 204 C C *** ASK WHETHER OF NOT TO RESTART 209 WRITE (6,604) 604 FORMAT (»6D0 YOU WISH TO RESTART WITH A NEW DATA SET') CALL INFREE (2058, ANSWER, 4) IF (ANSWER.EQ.YES) GO TO 207 IF (ANSWER.EQ.NO) GO TO 208 WRITE (6,601) GO TO 209 208 CALL PLOTND STOP END SUBROUTINE DATAIN (X, I, H ) C £***«**************************** ************************************** C* C* I N P U T O F D A T A P O I N T S C* C* C********************************************************************** DIMENSION X (1) , Y (1) C C *** READ IN THE NUMBER OF DATA POINTS CALL INFREE (4, N) C C *** READ IN THE ABSCISSA AND ORDINATE VALUES DO 100 L=1,N CALL INFREE (20, X (L) , Y (L) ) 100 CONTINUE C C /// DEBUGGING INFORMATION WRITE (3,301) N 301 FORMAT {•OTHE ', 13, * DATA POINTS ARE: 1) DO 400 L=1,N WRITE (3,300) X ( L ) , Y (L) 400 CONTINUE 300 FORMAT (' 1P2E15.6) RETURN END - 79 -SUBROUTINE KNOTIN (X, DELTA, M, N, HK, *) C c * C* I N P U T O F K N O T S E T C* DIMENSION X(1), DELTA (1) C C *** REQUEST THE DEGREE OF THE SPLINE WRITE (6,600) 600 FORMAT ('8DEGREE OF THE SPLINE») CALL INFREE (10, M) C C *** READ IN THE NUMBER OF KNOTS WRITE (6,601) 601 FORMAT ('SNUMBER OF KNOTS') CALL INFREE (10, K) C C *** CALCULATE THE TOTAL NUMBER OF KNOTS MK = M + K + 1 C C *** CHECK FOR AN UNDERDETERMINED SYSTEM IF (MK.LE.N) GO TO 200 WRITE (6,603) 603 FORMAT ('OSYSTEM IS UNDERDETERMINED:'/ * ' REDUCE DEGREE OF SPLINE OR NUMBER OF KNOTS') -RETURN 1 C C *** READ IN THE POSITION CF THE KNOTS 200 DO 100 1=1,K IM = I • H + 1 WRITE (6, 602) I CALL INFREE (26, DELTA (IM)) 100 CONTINUE 602 FORMAT ('8P0SITI0N OF KNOT (', 1 2 , » ) • ) C C *** CALCULATE THE POSITION OF THE SUPPLEMENTARY KNOTS H = (X(N) - X(1)) / (K + 1) M1 = H + 1 DO 101 1=1,M1 DELTA (I) = X(1) + (I - M1) * H DELTA (MK + I) = X (N) + ( 1 - 1 ) * H 101 CONTINUE RETURN END - 8 0 -SUBROUTINE DATAPL (X, Y, XMIN, XMN , XMAX, EX, YMIN, DY, N) C Q************************************************************ ********** c* C* P L O T O F D A T A P O I N T S C* '• C* Q**************************** ****************************************** DIMENSION X (1) , Y (1) DIMENSION XS (200) , YS (200) C C *** RESET THE PLOTTER CALL PLOT (15., 0., -3) C C *** PROTECT THE ORIGINAL DATA POINTS XMN = X (1) XMAX = X(N) DO 100 L=1,N XS(L) = X(L) YS(L) = Y (L) 100 CONTINUE C C *.** SCALE THE ABSCISSA AND ORDINATE VALUES CALL SCALE (XS, N, 10., XMIN, DX, 1) CALL SCALE (YS, N, 10., YMIN, EY, 1) C c *** PLOT THE DATA POINTS DO 101 L=1,N CALL SYMBOL (XS(L)-.035, YS(I)-.07, .14, • + •, 0., 1) 101 CONTINUE C *** PLOT THE AXES ORY = -YMIN / DY IF (ORY.LT.O.) ORY = 0. CALL AXIS (0., ORY, • -1, 10., 0., XHIN, DX) ORX = -XMIN / DX IF (ORX.LT.O.) ORX = 0. CALL AXIS (ORX, 0., » 1, 10., 90. f YMIN, DY) • C C *** DISPLAY THE PLOT CALL PSEND RETURN END - 81 -SUBROUTINE SPLNPL (ft, DELTA, XMIN, XMN, XMAX, DX, YMIN, # DY, M, MK) C c * C* P L O T O F S P L I N E C U R V E C* C* DIMENSION A(1), DELTA (1) COMMON /DW/ DOMEGA(50,5) C C *** CALCULATE THE SPLINE FUNCTION DENOMINATOR CALL OMEGA (DELTA, M, MK) C C *** RAISE THE PEN TO REPOSITION IT IPEN = 3 X = XMN XSTEP = .02 * DX NP = 50 * INT ( (XMAX - XMN) / DX + .5) DO 100 J=1,NP C C *** CALCULATE THE SPLINE FUNCTION VALUE AT THE CURRENT X VALUE Y = 0. DO 101 1=1,MK Y = Y + A (I) * SPLINE (DELTA, X, M, MK, I, 0) 101 CONTINUE C C *** SCALE THE X AND Y VALUES XP = (X - XMIN) / DX YP = (Y - YMIN) / DY C c *** IF THE POINT IS WITHIN THE RANGE IF ( (YP.GT. 10.5) .OR. (YP.LT.-.5) ) GO TO 200 C C *** THEN PLOT IT CALL PLOT (XP, YP, IPEN) IPEN = 2 GO TO 201 C C *** OTHERWISE RAISE THE PEN 200 IPEN = 3 201 X = XMN + J * XSTEP 100 CONTINUE C C *** DISPLAY THE SPLINE CURVE CALL PSEND RETURN END - 82-SUBROUTINE KNOTPL (DELTA, XMIN, DX, YMIN, DY , M, MK) C C* C* P L O T O F K N O T S E T C* C* Q * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DELTA(1) C C *** SCALE THE KNOT ORDINATE ONTO THE X-AXIS YP = -YMIN / DY IF (YP.LT.0.0) YP = 0.0 M2 = M + 2 DO 100 I=M2,HK C C *** SCALE THE KNOT XP = (DELTA (I) - XMIN) / DX C C *** PLOT THE KNOT CALL SYMBOL (XP-.0525, YP-.115, .21, «X», 0., 1) 100 CONTINUE C C *** DISPLAY THE KNOTS CALL PSEND RETURN END - 83 -SUBROUTINE FIXED (X, Y, YF, RES, A, DELTA, ERROR, ERRORP, M, N, MK C C********************************************************************** c * C* F I X E D K N O T L I N E A R C* L E A S T S Q U A R E S A P P R O X I M A T I O N C* c * Q********************************************************************** DIMENSION X(1), Y(1), YF(1), RES{1), A(1), E ELT A (1) DIMENSION S (200,50) C C *** CALCULATE THE SPLINE FUNCTION DENOMINATOR CALL OMEGA (DELTA, M, MK) C C *** FORM THE MATRIX OF SPLINE BASIS FUNCTION VALUES DO 100 1=1,MK DO 101 L=1,N S(L,I) = SPLINE (DELTA, X (L) , M, MK, I, 0) 101 CONTINUE 100 CONTINUE C C /// DEBUGGING INFORMATION WRITE (3,300) 300 FORMAT (»0AT », 13X, "THE HOUSEHOLDER MATRIX IS:*) DO 400 L=1,N WRITE (3,301) X (L) WRITE (3,302) ( S ( L , I ) , 1=1,MK) 400 CONTINUE 301 FORMAT (« », 1P8E15.6) 302 FORMAT (» «, 15X, 1P7E15.6) C C *** TRANSFER Y INTO A DO 102 L=1,N A (L) = Y(L) 102 CONTINUE C C /// DEBUGGING INFORMATION WRITE (3,303) 303 FORMAT (•OTHE HOUSEHOLDER VECTOR IS:») WRITE (3,301) (A (L) , L=1,N) C C *** PERFORM THE HOUSEHOLDER TRANSFORMATIONS CALL TRANS (S, A, N, MK) C C *** CALCULATE THE RESIDUALS AND LEAST SQUARES ERROR ERROR = RESERR (X, Y, YF, RES, A, DELTA, N, M, MK) CALL ITER (ERRORP, ERROR) RETURN END - 84 -FUNCTION RESERR (X, Y, YF, RES, A, DELTA, N, M, HK) C C******************************* **************** * ^ « ^ « 4 0 H 4 4 * * 4 < > » * * 4 * * * C* C* R E S I D U A L S A N E C* L E A S T S Q U A R E S E R R O R C* C* c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION X(1), Y(1), YF(1), RES(1), A(1), DELTA (1) DOUBLE PRECISION DY, DRES, DERR COMMON /DW/ DOMEGA(50,5) C C *** CALCULATE THE SPLINE FUNCTION DENOMINATOR CALL OMEGA (DELTA, M, MK) DERR = 0.D0 DO 100 L=1,N C C *** CALCULATE THE SPLINE FUNCTION VALUES AT TEE DATA POINTS DY = 0.D0 DO 101 1=1,MK DY = DY + A (I) * SPLINE (DELTA, X (L) , M, MK, I, 0) 101 CONTINUE YF (L) = DY C C *** CALCULATE THE RESIDUALS DRES = Y(L) - DY RES (L) = DRES C C *** CALCULATE THE LEAST SQUARES ERROR DERR = DERR + DRES * DRES 100 CONTINUE RESERR = DSQRT (DERR) C C /// DEBUGGING INFORMATION WRITE (3,300) RESEBR 300 FORMAT ('OTHE LEAST SQUARES ERROR IS •, 1PE15 .6 ) RETURN END - 85 -SUBROUTINE CONTRL (X, Y, YF, RES, A, DELTA, ERRORP, N, M, MK) C c * c * C* V A R I A B L E K N O T C O N T R O L C* C* DIMENSION X(1), Y(1), YF(1), RES(1), A(1), D ELT A (1) DIMENSION AP{200), DELTAP (50) , RESP(200), YFP(200) INTEGER YES /'YES'/* NC /'NO'// ANSWER /• •/ MMK = MK • M + 1 C C *** TRANSFER THE PREVIOUS KNOT SET DO 100 1=1,MMK DELTAP (I) = DELTA (I) 100 CONTINUE 601 FORMAT ('SPLEASE ANSWER YES OR NO') C C *** OBTAIN WHICH KNOT TO VARY 200 CALL VARYIN (I) C C *** REPOSITION THAT KNOT IM = I + M + 1 CALL VARYKT (DELTAP (IM), I) C C *** ASK WHETHER OR NOT TO REPOSITION OTHER KNOTS 204 WRITE (6,602) 602 FORMAT ('SDO YOU WISH TO VARY ANOTHER KNOT') CALL INFREE (2058, ANSWER, 4) IF (ANSWER.EQ.YES) GO TO 200 IF (ANSWER.EQ.NO) GO TO 203 WRITE (6,601) GO TO 204 203 CONTINUE C C /// DEBUGGING INFORMATION WRITE (3,300) 300 FORMAT ('OTHE NEW KNOT SET IS:') WRITE (3,301) (DELTAP ( I ) , 1=1,MMK) 301 FORMAT (• «, 1P8E15.6) C C *** RECALCULATE THE SPLINE APPROXIMATION CALL FIXED (X, Y, YFP, RESP, AP, DELTAP, ERROR, ERRORP, M, N,*MK) C C *** ASK WHETHER OR NOT TO PLOT THE NEW APPROXIMATION 207 WRITE (6,603) 603 FORMAT ('SDO YOU WISH TO SEE THE NEW APPROXIMATION') CALL INFREE (2058, -ANSWER, 4) IF (ANSWER. EQ. YES) .GO TO 205 IF (ANSWER.EQ.NO) GO TO 201 WRITE (6,601) GO TO 207 - 86 -2 0 5 C A L L D A T A P L ( X , Y , X H I N , X H N , X M A X , E X , Y M I N , E Y , N ) C A L L K N O T P L ( D E L T A P , X M I N , D X , Y M I N , D Y , M , M K ) C A L L S P L N P L ( A P , D E L T A P , X M I N , X M N , X M A X , D X , Y M I N , E Y , M , M K ) C C * * * D E C I D E W H E T H E R O R N O T T O K E E P N E W A P P R O X I M A T I O N 2 0 6 W R I T E ( 6 , 6 0 4 ) 6 0 4 F O R M A T ( « S D O Y O U W I S H T O C O N T I N U E W I T H T H E NEW A P P R O X I M A T I O N • ) C A L L I N F R E E ( 2 0 5 8 , A N S W E R , 4 ) I F ( A N S W E R . E Q . Y E S ) G O T O 2 0 8 I F ( A N S W E R . E Q . N O ) G O T O 2 0 9 W R I T E ( 6 , 6 0 1 ) G O T O 2 0 6 C C * * * I F S O , R E T A I N T H E N E W V A L U E S 2 0 8 D O 1 0 1 1 = 1 , M K A ( I ) = A P ( I ) 1 0 1 C O N T I N U E D O 1 0 2 1 = 1 , M M K D E L T A ( I ) = D E L T A P ( I ) 1 0 2 C O N T I N U E D O 1 0 3 L = 1 , N Y F ( L ) = Y F P ( L ) R E S ( L ) = R E S P ( L ) 1 0 3 C O N T I N U E C C * * * A S K I F A H A R D C O P Y I S W A N T E D C A L L S P L O U T ( X , Y , Y F , R E S , A , D E L T A , E R R O R , N , M , M K ) E R R O R P = E R R O R GO T O 2 0 1 , • C C * * * R E P L O T T H E O L D S P L I N E 2 0 9 C A L L D A T A P L ( X , Y , X M I N , X M N , X M A X , D X , Y M I N , D Y , N ) C A L L K N O T P L ( D E L T A , X M I N , D X , Y M I N , E Y , M , M K ) C A L L S P L N P L ( A , D E L T A , X M I N , X M N , X M A X , D X , Y M I N , D Y , M, M K ) 2 0 1 R E T U R N E N D - 87 -S U B R O U T I N E V A R Y I N ( I ) C* C* C* C* 6 0 0 I N P U T O F V A R I A E I E K N O T W R I T E ( 6 , 6 0 0 ) F O R M A T ( ' S K N O T T O B E V A R I E D 1 ) C A L L I N F R E E ( 1 0 , I ) R E T U R N END S U B R O U T I N E V A R Y K T ( D E L T A I , I ) C* C* C* C* V A R I A B L E K N O T V A L U E W R I T E ( 6 , 6 0 0 ) I 6 0 0 F O R M A T ( ' S N E W P O S I T I O N F O R KNOT ( ' , 1 2 , F ) ' ) C A L L IN F R E E ( 2 6 , D E L T A I ) R E T U R N E N D - 88 -SUBROUTINE ITER (ERRORP, ERROR) C C * * * * * * * * * * * * * * * * * * * « * * * * * * * 4 $ * * * * ^ c * C* I T E R A T I O N O U T P U T C* C* WRITE (6,600) ERRORP 600 FORMAT (•OTHE LEAST SQUARES ERROR OF THE PREVIOUS «, * 'APPROXIMATION WAS 1PE15.6) WRITE (6,601) ERROR 601 FORMAT (• THE LEAST SQUARES ERROR OF THE NEW «, * 'APPROXIMATION IS •, 1PE15.6) RETURN END - 89 -SUBROUTINE SPLOUT (X, Y, YF, RES, A, DELTA, ERROR, N, M, HK) C c * C* O U T P U T C * C* DIMENSION X{1), Y(1), Y F ( 1 ) f RES(1), A(1), DELTA (1) INTEGER YES /'YES'/* NO /'NO'/, ANSWER /• •/ INTEGER TITLE (10), BLANK /' •/, WUM /1/ C C *** ASK WHETHER OR NOT TO RETAIN RESULTS 200 WRITE (6,600) 600 FORMAT ('SDO YOU WISH TO RETAIN THE APPROXIMATION RESULTS*) CALL INFREE (2058, ANSWER, 4) IF .(ANSWER. EQ. YES) GO TO 201 IF (ANSWER.EQ.NO) GO TO 202 WRITE (6,601) 601 FORMAT ('&PLEASE ANSWER YES OR NO') GO TO 200 C C *** BLANK OUT THE TITLE 201 DO 101 1=1,10 TITLE (I) = BLANK 101 CONTINUE C C *** ASK FOR THE TITLE OF THE PLOT WRITE (6,602) 602 FORMAT (•STITLE FOR PLOT') CALL INFREE (2058, TITLE, 40) C C *** HARDCOPY PRINTOUT WRITE (8, 603) TITLE, NUM 603 FORMAT ('1', 10A4, 20X, 15) WRITE (8, 604) 604 FORMAT ('0', 5X, 'ABSCISSAE', 5X, 'ORDINATES', 3X, • 'FITTED ORDINATES', 3X, 'RESIDUALS') C C *** PRINT OUT THE DATA, FITTED RESULTS AND RESIDUALS DO 100 L=1,N WRITE (8, 605) X (L) , Y (L) , YF (L) , RES (L) 100 CONTINUE 605 FORMAT (» ', 1P5E15.6) C C *** PRINT OUT THE LEAST SQUARES ERROR WRITE (8,606) ERROR 606 FORMAT ('OTHE LEAST SQUARES ERROR IS', 1PE15.6) C C *** CALCULATE THE POLYNOMIAL COEFFICIENTS CALL COEFF (X, Y, A-, DELTA, N, M, M K ) C C *** PLOT A HARDCOPY OF THE APPROXIMATION CALL SYMBOL (1., 10., .14, TITLE, 0., 40) - 90 -F N U M = N U M C A L L N U M B E R ( 9 . , 1 0 . , . 1 4 , F N U M , 0 . , - 1 ) C A L L P S E N D C C A L L A G P L O T ( 1 0 . 0 , 6 0 0 0 ) N U M = N U M + 1 2 0 2 R E T U R N END - 91 -SUBROUTINE NORM (S, MK) C C* C* M A T R I X C O N D I T I O N N U M B E R C* C* C«**************************«*^ DIMENSION S(200,50), ST(50,50) C C *** SAVE THE ORIGINAL MATRIX DO 100 1=1,MK DO 100 L=1,MK ST(L,I) = S(L,I) 100 CONTINUE C C *** INVERT THE MATRIX AND CALCULATE ITS CONDITION NUKEER CALL INVERT (ST, MK, 50, DET, COND) C C /// DEBUGGING INFORMATION WRITE (3,300) 300 FORMAT ('OTHE INVERTED TRANSFORMED HOUSEHOLDER MATRIX IS:') DO 400 1=1,MK WRITE (3,301) (ST(L, I ) , L=1,MK) WRITE (3,302) 400 CONTINUE 301 FORMAT (• ', 1P8E15.6) 302 FORMAT (' ») C C *** PRINT OUT THE CONDITION NUMBER WRITE (6,600) COND 600 FORMAT {•OTHE CONDITION NUMBER IS', 1PE15.6) RETURN END r- 92 -SUBROUTINE TRANS (S, B f N, MK) C c* C* H O U S E H O L D E R T R A N S F O R M A T I O N C* C O N T R O L R O U T I N E C* C* C****************************************** * * * * * * * * * * * * * * * * * * 4 * * * * * * * * * DIMENSION S (200,50), V (200) DIMENSION B(1) REAL KSQ C C *** AUGMENT THE MATRIX S BY B MK1 = MK + 1 DO 100 L=1,N S (L,MK 1) = B (L) 100 CONTINUE C C *** PERFORM THE HOUSEHOLDER TRANSFORMS DO 101 1=1,MK C C *** TRANSFER CURRENT COLUMN OF HOUSEHOLDER MATRIX DO 102 L=1,N V(L) = S(L,I) 102 CONTINUE C C *** DERIVE HOUSEHOLDER TRANSFORMATION VECTOR CALL UVEC (V, N, I , KSQ) C C *** PERFORM HOUSEHOLDER TRANSFORMATION ON THE VECTOR CALL HOUSE (S, V, N, MK1, I, KSQ) 101 CONTINUE C C /// DEBUGGING INFORMATION WRITE (3,300) 300 FORMAT ('OTHE TRANSFORMED HOUSEHOLDER MATRIX IS:') DO 400 L=1,N WRITE (3,301) ( S ( L , I ) , 1=1,MK) WRITE (3,302) 400 CONTINUE 301 FORMAT (« • , 1P8E15.6) 302 FORMAT (» ') C C *** CALCULATE THE CONDITION NUMBER OF THE MATRIX CALL NORM (S, MK) C C *** PERFORM THE BACKWARDS SUBSTITUTION CALL SOLVE (S, B, N, MK) RETURN END SUBROUTINE UVEC (V, N # I, KSQ) C Q****************************************************** **************** C* C* H O U S E H O L D E R V E C T O R C* C* Q************************ ********************************************** DIMENSION V(1) DOUBLE PRECISION DNORH REAL KSQ C C *** COMPUTE THE NORM OF V DNORM = 0.D0 DO 100 L=I,N DNORM = DNORM + V (L) * V (L) 100 CONTINUE DNORM = DSQRT (DNORM) C C *** CALCULATE THE COEFFICIENT NSGN = 1 IF (V(I).LT.O.) NSGN = -1 KSQ = 2 * DNORM * (DNORM • NSGN * V(I)) C C *** CREATE THE HOUSEHOLDER VECTOR V(I) = V(I) + NSGN * DNORM RETURN END - 94 -SUBROUTINE HOUSE (S, U, N, HK1, I, KSQ) C C* C* H O U S E H O L D E R T R A N S F C R M A T I C N C* c* Q*************************************** ****************** ************* DIMENSION S(200,S0), U(1) REAL KSQ C *** PERFORM HOUSEHOLDER TRANSFORMATION ON REQUIRED COLUMNS OF MATRIX DO 100 J=I fMK1 C C *** MULTIPLY THE TRANSPOSE OF THE VECTOR C C BY THE COLUMN OF THE MATRIX UT = 0. DO 101 L=I,N UT = UT + U (L) * S (L, J) 101 CONTINUE UT = 2. * UT / KSQ C C *** PERFORM THE HOUSEHOLDER TRANSFORMATION DO 102 1=1,H S (L, J) = S (L,J) - UT * U (L) 102 100 CONTINUE CONTINUE RETURN END i - 95 -SUBROUTINE SOLVE (S, B, N, HK) C C***************4*******4***4**4**** C* C* B A C K W A R D S S U B S T I T U T I O N C* C* DIMENSION S (200, 1) , B (1) , MK1 = MK + 1 C C *** COMPUTE THE FINAL ELEMENT B (MK) = S(MK,MK1) / S (MK , MK) IF (MK.EQ.1) GO TO 200 MKM1 = MK - 1 DO 100 I=1,MKM1 MKI = MK - I -B(MKI) = 0. MKI1 = MKI • 1 DO 101 J=MKI1,MK B (MKI) = B (MKI) + B (J) * S(MKI,J) 101 CONTINUE B (MKI) = (S(MKI,MK1) - B (MKI) ) / S(MKI,MKI) 100 CONTINUE C C *** TRANSFER THE REMAINING ELEMENTS INTO B 200 IF (N.EQ.MK) GO TO 201 DO 102 L= MK1 , N B(L) = S(L,MK1) 102 CONTINUE 201 RETURN END - 96 -SUBROUTINE OPT (X, Y, YF, RES, A, DELTA, EBRCRP, N, M, MK) C C********************************************************************** c* C* O P T I M I Z A T I O N O F K N O T S C* c* Q********************************************************************** DIMENSION X {1) , Y(1), YF(1), RES(1), A(1), BELT A (1) DIMENSION AP(200), DELTAP (50) , YFP(200), RESP (200) DIMENSION P (50,75) INTEGER YES /•YES»/» NO /'KOV, ANSWER /• »/ EXTERNAL RESERR, GDELTA, HDELTA C C *** TRANSFER KNOTS FOR OPTIMIZATION MMK = MK + M + 1 DO 300 1=1,MMK DELTAP (I) = DELTA (I) 300 CONTINUE C C *** TRANSFER PARAMETERS FOR OPTIMIZATION DO 301 1=1,MK AP(I) = A (I) 301 CONTINUE C C *** TRANSFER THE KNOTS TO THE SIMPLEX DO 102 J=1,MMK P ( J , 1) a DELTAP (J) 102 CONTINUE C C *** COMPUTE THE NUMBER OF POINTS IN THE SIMPLEX IF (MMK.GE.10) NPLUS = (3 * MMK + 1) / 2 IF (MMK.LT.10) NPLUS = 2 * MMK C C *** MINIMIZE THE LEAST SQUARES ERROR CALL COMPLX (X, Y, YFP, RESP, AP, DELTAP, EBRORP, # P, 50, MMK, NPLUS, MMK, M, N, MK, 1.3, -1, # 100, 100, 1, -.0001, RESERR, GDELTA, # HDELTA, G.DELTA, 5200, S200) C C *** RECALCULATE THE FIXED KNOT APPROXIMATION CALL FIXED (X, Y, YFP, RESP, AP, DELTAP, ERROR, EBRORP, M, N, MK) GO TO 207 C C *** IF KNOT OPTIMIZATION FAILED, PRINT AN ERROR MESSAGE 200 WRITE (6,801) 801 FORMAT (» 0OPTIMIZATION FAILED. *) C C *** ASK WHETHER OR NOT TO PLOT THE NEW APPROXIMATION 207 WRITE (6,603)' 603 FORMAT (•SDO YOU WISH TO SEE THE NEW APPROXIMATION') CALL INFREE (2058, ANSWER, U) IF (ANSWER.EQ.YES) GO TO 205 IF (ANSWER.EQ.NO) GO TO 201 - 97 -WRITE (6,601) 601 FORMAT (•SPLEASE ANSWER YES OR NO') GO TO 207 205 CALL DATAPL (X, Y, XMIN, XMN, XMAX, DX, YMIN, DY, N) CALL KNOTPL (DELTAP, XMIN, DX, YMIN, DY, M, MK) CALL SPLNPL (AP, DELTAP, XMIN, XMN, XMAX, DX, YMIN, DY, M, MK) C C *** DECIDE WHETHER OR NOT TO KEEP NEW APPROXIMATION 206 WRITE (6,604) 604 FORMAT ('SDO YOD WISH TO CONTINUE WITH THE NEW APPROXIMATION') CALL INFREE (2058, ANSWER, 4) IF (ANSWER.EQ.YES) GO TO 208 IF (ANSWER.EQ.NO) GO TO 209 WRITE (6,601) GO TO 206 C C *** IF SO, RETAIN THE NEW VALUES 208 DO 101 1=1,MK A (I) = AP(I) 101 CONTINUE DO 302 1=1,MMK DELTA (I) = DELTAP (I) 302 CONTINUE DO 103 L=1,N YF (L) = YFP(L) RES (L) = RESP (L) 103 CONTINUE C • . C *** ASK IF A HARDCOPY IS WANTED CALL SPLOUT (X, Y, YF, RES, A, DELTA, ERROR, N, M, MK) ERRORP = ERROR GO TO 201 C C *** REPLOT THE OLD SPLINE 209 CALL DATAPL (X, Y, XMIN, XMN, XMAX, EX, YMIN, DY, N) CALL KNOTPL (DELTA, XMIN, DX, YMIN, DY, M, MK) CALL SPLNPL (A, DELTA, XMIN, XMN, XMAX, EX, YMIN, DY, M, MK) 201 RETURN E N D ' - 98 -FUNCTION GDELTA (T, H, MMK, J) C C * * * * * * * * * * * * * * * * * * * * * * * 4 * 4 * 4 4 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C* C* L O W E R B O U N D C O N S T R A I N T C* c * £********************************************************************** DIMENSION T ( 1 ) C C *** SET THE LOWER BOUNDS ON THE KNOTS IF ( (J.LE. M+1) .OR. (J.GT. (MMK-(M + 1) )) ) GOTO 100 H = T (MMK- (M + 1) ) - T (M + 1) GDELTA = T(J-1) + .0001 * H GO TO 999 100 GDELTA = T (J) 999 RETURN END FUNCTION HDELTA (T, M, MMK, J) C C********************************************************************** C* C* U P P E R B O U N D C O N S T R A I N T C* C* c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION T (1) C c *** SET THE UPPER BOUNDS ON THE KNOTS IF ( (J.LE.M + 1) .OR. (J.GT. (MMK-(M + 1)))) GOTO 100 H = T (MMK- (M+1) ) - T (M + 1) HDELTA = T ( J + 1) - .0001 * H GO TO 999 100 HDELTA = T(J) 999 RETURN END - 99 -SUBROUTINE COEFF (X, Y, A, DELTA, N, M, MK) P O L Y N O M I A L C O E F F I C I E T S C* C* C* C* **** DIMENSION X(1), Y ( 1 ) r A(1), DELTA (1) DIMENSION C(50,4) COMMON /DW/ DOMEGA(50,5) C C *** CALCULATE THE SPLINE DENOMINATOR CALL OMEGA (DELTA, M, MK) C C *** DETERMINE THE NUMBER OF PIECEWISE POLYNOMIALS K = MK - M M1 = M + 1 C c #** DETERMINE THE COEFFICIENTS OF THE PIECEWISE POLYNOMIALS JFACT = 1 DO 100 J=1,M1 J1 = J - 1 DO 101 1=1,K IM = I + M C(I,J) = 0.0 DO 102 L=1,MK C(I,J) = C(I,J) + A (L) * * SPLINE (DELTA, DELTA (IM) , H, MK, L,\ J1) 102 CONTINUE C(I,J) = C(I,J) / JFACT 101 CONTINUE JFACT = JFACT * J 100 CONTINUE C C *** PRINT OUT THE COEFFICIENTS WRITE (8,800) FORMAT («0 KNOT', 5X, 'LOCATION', 5X, 'POLYNOMIAL POWER *, * 5X, 'POLYNOMIAL COEFFICIENT') DO 103 1=1,K 11 = I - 1 IM = I + M WRITE (8,801) 11, DELTA (IM) DO 104 J=1,M1 J1 = J - 1 WRITE (8,802) J1, C(I,J) CONTINUE CONTINUE MK1 = MK + 1 WRITE (8,801) K, DELTA (MK1) FORMAT (' *, 14, 2X, 1PE15.6) FORMAT (• ', 30X, 11, 15X, 1PE15.6) RETURN END 800 104 103 801 802 - 100.-F U N C T I O N S P L I N E ( D E L T A , X , H, M K , I , J ) C********************************: c * S P L I N E F U N C T I O N A N D D E R I V A T I V E S C * C * C * C * C c * * * C c * * * 1 0 0 C c * * * 1 0 1 2 0 0 D I M E N S I O N D E L T A ( 1 ) C O M M O N / D W / D O M E G A ( 5 0 , 5 ) D O U B L E P R E C I S I O N D S P L N , D X , D D E L T A , B Y M C O N = 1 M J = N - J D S P L N = 0 . D 0 C H E C K I F T H E P O I N T I S I N T H E R E G I O N O F S U P P O R T M 1 = M • 1 I F ( ( X . L T . D E L T A ( I ) ) . O R . ( X . G T . D E L T A ( I + M 1 ) ) ) G O T O 2 0 0 E V A L U A T E T H E S P L I N E A T T H E G I V E N P O I N T D X = X M2 = M + 2 " DO 1 0 0 L = 1 , M 2 L L = I + L - 1 D D E L T A = D E L T A ( L L ) D Y = D D E L T A - D X I F ( D Y . L E . O . D O ) G O T O 1 0 0 D S P L N = D S P L N + D Y * * M J / D O M E G A ( I , L ) C O N T I N U E I F ( ( I . E Q . 1 ) . O R . ( I . E Q . M K ) ) D S P L N = M 1 * D S P L N C O M P U T E T H E C O N S T A N T T E R M I F ( J . E Q . O ) G O T O 2 0 0 M J 1 = M J + 1 DO 1 0 1 L = M J 1 , M M C O N = - M C O N * L C O N T I N U E S P L I N E = M C O N * D S P L N R E T U R N E N D - 101 -SUBROUTINE OMEGA (DELTA, M, MK) C* C* S S P L I N E F U N C T I O N D E N O M I N A T O R C* C* DIMENSION DELTA (1) DOUBLE PRECISION DDELTA COMMON /DW/ DOMEGA(50,5) C C *** CALCULATE THE MATRIX OF VALUES FOR THE DENOMINATOR C *** OF THE SPLINE FUNCTION M2 = M + 2 DO 100 1=1,MK DO 100 J=1,M2 DOMEGA(I,J) = 1. J J = I + J - 1 DO 100 1=1,M2 IF (L.EQ.J) GO TO 100 LL = I + L - 1 DDELTA = DELTA (JJ) - DELTA (LL) DOMEGA(I,J) = DOMEGA(I,J) * DDELTA 100 CONTINUE C C /// DEBUGGING INFORMATION WRITE (3,300) 300 FORMAT (*OTHE BASIS FUNCTION DENOMINATOR VALUES ARE:») DO 400 1=1,MK WRITE (3,301) (DOMEGA ( I , J ) , J=1,M2) 400 CONTINUE 301 FORMAT (' •, 1P5E15.6) RETURN END - 102 -Appendix B INTERACTIVE SPLINE APPROXIMATION Purpose This system enables the use. to perform l e a s t squares approx-imations using s p l i n e functions. These approximations are performed i n t e r a c t i v e l y with the aid of a graphics terminal. The approximations are displayed immediately a f t e r computation and theuuser can respecify knot locations and r e c a l c u l a t e the f i t . Features are included to optimize the knot set and change the number of knots i n the knot set. Type of Routine This i s a self-contained program written i n FORTRAN IV. How to Use To run th i s program under MTS at the Adage Graphics Terminal; enter the command: $RUN IRAM:SPLINE+AGT:BASIC 3=debugfile 4=datafile 8 = p r i n t f i l e 9 - p l o t f i l e where debugfile contains the debugging information. Unless the system runs into problems with the approximation t h i s should be set to *DUMMY* . d a t a f i l e i s the f i l e containing the input data. The format i s described i n Section A: Data Input Format. - 103 -p r i n t f l l e i s the f i l e . c o n t a i n i n g the printed output from an approx-imation corresponding to a p l o t . p l o t f i l e i s the f i l e . c o n t a i n i n g the pl o t information to be retained for a hardcopy. Upon completion of the : program a hardcopy of the printouts and plots which were saved can be obtained by i s s u i n g the following commands: $C0PY p r i n t f i l e *PRINT* , $R PL0T:Q EAR=plotfile where ' p r i n t f i l e ' contains the p r i n t output described above and ' p l o t f i l e ' contains the p l o t information. Description A. Data Input Format The data f i l e has the following structure: The f i r s t l i n e states the number of points i n the f i r s t data set followed by data l i n e s with the sequence ABSCISSA ORDINATE. The remaining data sets follow with an i d e n t i c a l format. The data i t s e l f i s i n free-format providing that at l e a s t one blank delimits each entry. - 104 -B. Displays Because of the i n t e r a c t i v e c a p a b i l i t i e s of the system, displays play an i n t e g r a l part i n the structure. Inoorder to best describe the flow of control through the system a sample run follows with comments inserted to augment the display information. Example; Demonstration run $SIG IRAM 'DEMONSTRATION RUN1, PASSWORD Signon information $C0MMENT LOAD GRAPH IF NECESSARY $C0PY AGT:GRAPH > AGTI $COMMENT SPLINE APPROXIMATION PROGRAM $RUN IRAM:SPLINE+AGT;BASIC 3=*DUMMY* 4=DATA 8=-PRINT 9=-PL0T EXECUTION BEGINS (One data set from unit 4 i s immediately read i n . A pl o t of t h i s set appears on the graphics display.) DEGREE OF THE SPLINE? 3 (The degree of the piecewise polynomials wanted i s requested, most frequently used degree i s 3; that i s , a cubic s p l i n e approximation.) NUMBER OF KNOTS? 4 (The number of i n t e r n a l knots wanted i s requested. Boundary knots and supplementary knots are computed by the program.) - 105 -POSITION OF KNOT ( 1)? .2 POSITION OF KNOT ( 2)1 .4 POSITION OF KNOT ( 3)? .6 POSITION OF KNOT ( 4 ) ? .8 (The l o c a t i o n of the abscissa of each knot i s requested. The knots must be entered i n a s t r i c t l y increasing sequence. Immediately following the input of, a l l the knots the knot set i s plotted along the x-axis.) (A f i x e d knot l e a s t squares approximation i s now computed. The r e s u l t i n g s p l i n e function i s o v e r l a i d on the graphics d i s p l a y ) . THE LEAST SQUARES ERROR OF THE PREVIOUS APPROXIMATION WAS 1.OOOOOOE 03 (The i n i t i a l l e a s t squares error i s a r b i t r a r i l y set to 1000 since no previous approximation was done.) THE LEAST SQUARES ERROR OF THE NEW. APPROXIMATION IS 1.537876E-02 DO YOU WISH TO RETAIN THE APPROXIMATION RESULTS? NO (A message asking whether or not to keep a hardcopy of the present approximation set i s printed.) DO YOU WISH TO VARY THE KNOT SET? YES (A message asking whether or not to re p o s i t i o n any of the knots i s printed.) KNOT TO BE VARIED? 1 NEW POSITI'ONFFOR KNOT ( 1)? .3 DO YOU WISH TO VARY ANOTHER KNOT? YES KNOT TO BE VARIED? 4 NEW POSITION FOR KNOT ( 4)? .7 DO YOU WISH TO VARY ANOTHER KNOT? NO (A f i x e d knot l e a s t squares, approximation(is computed using the new knot - 106 -set. The r e s u l t i n g l e a s t squares error information i s printed.) THE LEAST SQUARES.ERROR OF THE PREVIOUS APPROXIMATION WAS 1.537876E-03 THE LEAST SQUARES ERROR OF THE NEW APPROXIMATION.IS 9;576172E-03 DO YOU WISH TO SEE THE NEW APPROXIMATION? YES (If the answer to t h i s question i s a f f i r m a t i v e the new knot set and spl i n e curve replaces the one currently being displayed on the graphics terminal.) DO YOU WISH TO CONTINUE WITH. THE NEW' APPROXIMATION? YES (The values of the new approximation are transferred to the proper arrays.) DO YOU WISH TO RETAIN THE APPROXIMATION RESULTS? YES TITLE FOR PLOT? TEST DATA'' (A t i t l e i s printed on each graph and printout. These are also l a b e l l e d with a number which i s the number of hardcopies obtained thus far during the run. This method enables corresponding plots and printouts from a p a r t i c u l a r run to be uniquely i d e n t i f i e d . ) DO YOU WISH TO VARY THE KNOT SET? NO DO YOU WISH TO OPTIMIZE THE KNOT SET? NO (A request i s printed as to whether oranot to optimize the current knot set. WARNING: optimization i s time-consuming and expensive. It is. best done non-interactively.) DO YOU WISH TO RESTART WITH A NEW KNOT SET? NO (A request, i s issued whether or not to r e s t a r t with an expanded or contracted knot set. If the answer i s a f f i r m a t i v e the data set i s replotted and the program r e s t a r t s from the beginning.) DO YOU WISH TO RESTART WITH A NEW DATA SET? NO (If the answer i s a f f i r m a t i v e , the program reads another data set from - 107 -unit 4 and returns to the beginning of the program). PLOTTING WILL TAKE APPROX. 2:MINS 32 SECONDS STOP 0 EXECUTION TERMINATED $C0PY -PRINT *PRINT-* $RUN PL0T:Q PAR=^PL0T $SIG C. Output Approximation information corresponding to a plo t i s put i n a f i l e . The format of each approximation printout i s as follows: TITLE NUMBER ABSCISSAE ORDINATES FITTED ORDINATES RESIDUALS for one data set THE LEAST SQUARES, ERROR IS KNOT LOCATION POLYNOMIAL POWER POLYNOMIAL COEFFICIENT This p r i n t s the c o e f f i c i e n t s of every power for each piece-wise polynomial between each knot p a i r . Plots obtained for the p l o t t e r correspond to this, output. They can be matched by the t i t l e and l a b e l number. Re s t r i c t i o n s The number of data points must be les s that 200 . The degree of the sp l i n e msut be les s than 3 . The number of knots must be l e s s that 42 . 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051769/manifest

Comment

Related Items