NASA Contractor Report 178001 

MINIMAL NORt1 CONSTRAINED INTERPOLATION 

Larry Dean Irvine 

OLD DOMINION UNIVERSITY 

Norfolk, Virginia 

Grant NAGl-363 

September 1985 

NI\SI\ 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665 

NASA-CR-178001 

19860004487 

'"if"" --I"( ,/" -I' 

,i _ ;_ '. --

NOV 2 6 1985 

Lh:~L.~_' I\_:"LM.I)"='/ ~L .. feR 

L'CP' ';(Y NAS,\ 

'-AMPTQN VI~GIN',\ 

111111111111111111111111111111111111111111111 

NFOI022 
Chapter 1. 

Chapter 2. 

TABLE OF CONTENTS 

THE NATURAL SPLINE INTERPOLANT . . . 

A MINIMAL NORM INTERPOLATION PROBLEM 

IN THE L [a,b] SPACES ..•• 

p 

Chapter 3. THE CONVEX SPLINE INTERPOLANT 

Chapter 4. THE SHAPE-PRESERVING SPLINE INTERPOLANT 

Chapter 5. CONSTRAINED MINIMIZATION IN A DUAL SPACE 

LIST OF REFERENCES • . . . • . . • . 

APPENDIX A A PROGRAM FOR CONSTRUCTING THE NATURAL CUBIC 

Page 

1 

16 

20 

38 

55 

76 

SPLINE INTERPOLANT TO GIVEN DATA . • . • • . •. 77 

APPENDIX B A PROGRAM FOR CONSTRUCTING THE SHAPE-PRESERVING 

CUBIC SPLINE INTERPOLANT TO GIVEN DATA • • • •• 82 

1. 
This Page Intentionally Left Blank 
Acknowledgements 

I wlsh to express my appreclatlon to Dr. Phl11P Smlth, Dr. Robert 

Smlth, and Dr. John Swetlts for thelr advlce, asslstance, and encouragement durlng my doctoral work and research. I wlsh to thank Dr. 

Surendra N. T1WarJ for the flnanclal support WhlCh I recelved from the 

Instltute for Computatlonal and Applled Mechanlcs. 

I also wlsh to thank Barbara Jeffrey for typlng the dlssertatl0n. 

Lll 
This Page Intentionally Left Blank 
ABSTRACT 

MINIMAL NORM CONSTRAINED INTERPOLATION 

Larry Dean Irvine 

Old Dominion Unlversity, 1985 

Director: Dr. PhlllP W. Smith 

In computatlonal fluid dynamlcs and ln CAD/CAM a physlcal boundary, 

usually known only discreetly (say, from a set of measurements), must 

often be approximated. An acceptable approximation must, of course. 

preserve the salient features of the data (convexity, concavlty. etc.) 

In thlS dlssertation we compute a smooth interpolant which is locally 

convex where the data are locally convex and is locally concave where 

the data are locally concave. 

Such an interpolant is found by posing and solving a minimlzation 

problem. The Solutlon is a piecewise CUblC polynomlal. We actually 

solve this problem indirectly by using the Peano kernel theorem to 

recast this problem into an equivalent minimization problem havlng the 

second derivatlve of the interpolant as the solution. 

We are then led to solve a nonllnear system of equations. We 

show that with Newton's method we have an exceptlonally attractlve and 

efflcient method for solvlng thlS nonlinear system of equations. 

We dlsplay examples of such lnterpolants as well as convergence 

results obtalned by uSlng Newton's method. We llSt a FORTRAN program 

to compute these shape-preserving interpol ants. 

Next we consider the problem of computing the interpolant of 

mlnlmal norm from a convex cone in a normed dual space. This lS an 

extension of de Boor's work on minimal norm unconstralned interpolation. 

v 
1. The Natural Spl1ne Interpolant 

We cons1der the problem of comput1ng an 1nterpolant to glven data. 

Throughout our d1SCUSS10n we shall denote the data 

1=1,2, ••• ,n 

where a = tl < t2 < ••• < tn = band 1n th1S chapter we place no 

restr1ctl0ns on the numbers y. There are, of course, many such 

1 

lnterpolants wh1ch we can form. For example, we can calculate the 

un1que polynomlal p of order n (degree n-1 or less) wh1ch 1nterpolates 

the data. However, as pOlnted out In [deB(l), chapter 2], for large 

n (and espec1ally for equally spaced p01nts t ) the po1ynom1al 1nter1 

polant lS notorl0US for large changes In 1tS flrst der1vatlve near the 

endpolnts. Flgure (1.1) dlsplays the polynomlal lnterpolant to the 

functlon 

f(t) 1 - sln(7 7T t) 

2 

at the p01nts tl = (1-1)/10 for 1 = 1,2 ..... 11. Slnce 0 ~ Yl ~ 1 for 

each 1, we expect a good lnterpolant to remaln reasonably close to 

these bounds. However, because of ltS behav10r near the endpolnts, 

the polynom1al lnterpo1ant falls to model the data well. ThlS behavlor 

lS tYPlcal of hlgh-order polynom1al 1nterpo1ants. 

In order to decrease the unnaturally large changes In the f1rst 

derlvatlve characterlstlc of the polynom1al 1nterpo1ant, we w1sh to 

calculate the lnterpolant whlch "bends" the least over all sUltable 

lnterpolants. The norm of the second derlvatlve of an lnterpolant 

wlll furnlsh a measure of the bend1ng of the lnterpolant so we pose a 

mlnlmlzatlon problem on L

2(2)[a,b], the Sobolev space of functlons w1th 
2 

~o 

Flgure (1.1): The Polynomlal Interpolant. 
'>ccond derlvat1ves 1n the normed llncar space L2la,b]. Let A delloLe 

the set of all lnterpolant& 1n the Subolev &pace. We conslder the 

m1n1m1zat1on problem .. 

F1nd f"E A such that Ilfi~(2) 112 ~ IIf(2) II 2 for all fE A. (1.1) 

We shall see that the solut10n to (1.1) lS p1ecew1se CUblC w1th two 

cont1nuous der1vat1ves; that lS 

p (t) 

1 

for 1 = 1,2, ... ,n-1 where p lS a CUb1C po1ynom1a1 and f* lS 1n 

1 

2 C [a,b]. We follow the pattern 1n [deB(l), chapter 5], tak1ng ad vantage of the fact that L2[a,b] 18 not only a normed llnear space, but 

also a Hllbert space wlth an lnner product deflned by 

(f,g) = Ibf(t)g(t)dt 

a 

for any two elements f and g In L2[a,b]. 

Assume f lS an element of A. (The set A lS nonempty Slnce lt 

contalns the polynomlal lnterpolant.) We shall use the Peano kernel 

theorem to obtaln a set of equat10ns for f(2). By the Fundamental 

Theorem of Calculus we have 

t 

fCt) = fCa) + 1 f(1)Cs)ds 

a 

(1. 2) 

We lntegrate J rt f (1) (s)ds by parts not1ng that f udv = uv - I vdu. 

a 

Let 

u(s) f(l)(s) and dyes) = ds 

so that 

dues) f(2)(s)ds and yes) = -(t-s) 
4 

where t lS a constant. Hence 

l~(l)(S)dS = (t-a)f(l)(a) + !.<t-S)f(2)(S)dS 

a a 

and so (1.2) becomes 

f(t) ( 1.3) 

where ql(t) = f(a) + f(l)(a)(t-a). (Th1S lS actually a Taylor's ser1es 

w1th 1ntegral rema1nder.) 

To acqU1re constant llm1t of 1ntegrat1on we can wr1te (1.3) as 

b 

f(t) = q (t) + l(t-s) f(2)(s)ds 1 + 

a 

where (h)+, the pos1t1ve part of the funct10n h, lS def1ned by 

h( t) ~ 0 

h(t) ;:;; O. 

(1.4) 

Now we cons1der the d1v1ded d1fference operator. G1ven a funct10n 

g and a set of p01nts {T,T 1, ... ,T }, the m-th d1v1ded d1fference 

1 1+ l+m 

of g - denoted by [T,T 1, ... ,T ]g(e) - lS the lead1ng coeff1c1ent 1 1+ l+m 

of the polynom1al of order m+l wh1ch 1nterpolates g at T,T 1, ..• T 1 1+ l+m 

(and hence lS a funct10n of T ,T +l, ... ,T +). The recurS1ve relat10ns 

1 11m 

[T ]g(e) = geT ) p p 

[T T ]g(e)-[T T ]g(e) l+l, ... ,l+m 1, ... , l+m-l (1.5) 

T - T l+m 1 

hold 1f T ~ T (wh1ch we assume for our data). Presently we are l+IlI 1 

1nterested 1n the case m=2. Equatlon (1.5) becomes (w1th T = t ) 1 1 
g(tH2 )-g(tHI ) 

tI+2 - tHl 

Wh1Ch lS computable for 1=1,2, .•. ,n-2. 

5 

get l)-g(t) 1+ 1 

0.6) 

Not1ce that [T ,T 1, ... ,T Jp(-) = 0 1f P lS a polynom1al of 1 1+ l+m 

order m or less (degree m-l or less). (From equat10n (1.6) we see 

that (t - t )[t t t ]g(o) measures a d1fference 1n slopes; the L+2 1 l' 1+1' 1+2 

d1fference 1n slopes be1ng zero 1f g lS llnear.) 

Now we apply the (scaled) second-d1v1ded d1fference operator 

(t 2 - t )[t ,t l,t 2] to (1.4) and 1nterchange the order of the 1+ 1 1 1+ 1+ 

1ntegral and d1v1ded d1fference operators to obta1n 

where 

and g = f(2). 

b 

d = (g(s)N (s)ds 1,2 ~ 1 

d 1,2 

= Y1+2 - YHI 

t1+2 - tHl 

1=1,2, ••• ,n-2 (1. 7) 

(1.8) 

(t1+2 - s)+ - (t1+l - s)+ 

tH2 - tHl 

(t1+l - s)+ - (t1 - s)+ 

t 1+l - t1 

We call N 2 the (normal1zed) 11near B-spl1ne (or 1, 

(1. 9) 

B-spl1ne of order 2) w1th knots t , t 1 and t 2' 1 1+ 1+ The graph of N 1,2 

lS d1splayed 1n f1gure (1.2). 
6 

Flgure (1.2): The Norma11zed Llnear B-spllne. 
7 

We have shown that 1f f lS an 1nterpolant 1n the Sobolev space 

(f E A), then g = f(2) sat1sf1es (1.7). Let the set B cons1st of all 

funct10ns Wh1Ch are 1n L

2

[a,b] and wh1ch sat1sfy (1.7). 

Now cons1der the problem 

F1nd g~~ E B such that II g~~112 ~ II gl12 for all g E B (1.10) 

A un1que Solut1on eX1sts Slnce (1.10) lS a m1n1ma1 norm problem over 

a nonempty closed convex set 1n a H1lbert space. Furthermore, the 

(2) Solut1ons of problems (1.1) and (1.10) are related V1a gi~ = f* • 

Hence, to compute fi~ we can f1rst calculate g* and then 1ntegrate g* 

tW1ce. Slnce much of our emphas1s w1ll be on g*, rather than f*, we 

shall call g* the 1nterpolant of m1n1mal norm. 

For brev1ty we denote the 1ndex m n-2, the B-sp11ne N = N 2' 

1 1, 

and the d1v1ded d1fference d 

1 

d 2' We also def1ne the vector-valued 

1, 

b 

(Tx) = J x(t)N (t)dt 1 1 

a 

1=1,2, ••• ,m. 

To solve problem (1.10) we shall show that g*, the 1nterpo1ant of 

m1n1ma1 norm, 1S the 1ntersect1on of two spec1f1c sets--one an orthogona1 complement of a subspace and the other a translate of a subspace 

--In L2[a,b] V1a a var1at1on of the PrO]ect1on Theorem. If W lS a 

closed subspace of a H11bert space Hand 1f x 1S an arb1trary element 

of H, then the PrO]ect1on Theorem states that there eX1sts a un1que 

element w 1n W sat1sfY1ng o 

II x - \Y II ~ II x - \Y II 

o 

and character1zed by 

for all \Y E W (1.11) 
(x - w ,w) = 0 

o 

for all w e:: W. 

8 

Hence x - w LS 1n W ~. the orthogonal complement of W. The proof of 

o 

the PrO]ect1on Theorem can be found 1n any book deal1ng w1th H1lbert 

spaces (for example, [L, page 517]). The next propos1t1on w1ll serve 

as the actual form of the PrO]ect1on theorem Wh1Ch we shall usc. 

Propos1t1on ([L, page 64]): Let W be a closed subspace 1n a H1lbert 

space H. For a hxed element x 1n H dehne V: = x + W. Then there 

eX1sts a un1gue element x 1n V of m1n1mal norw. Furthermore, x 1S - 0- - 0-

.L 

1n W 

(The translate V 1S called an aff1ne set or llnear var1ety.) Not1ce 

that x 1S the 1ntersect1on of the orthogonal complement of Wand the 

o 

translate V of W. In fact, (1.11) reveals that x = x - w . 

o 0 

Def1ne 

W: = {z E: L2 [a, b] : Tz = 8} 

Wh1Ch 1S a closed subspace 1n L2[a,b]. Let ge:: L2[a,b] be any element 

such that Tg = 1. (Equlvalently, let g be any element of B.) Then 

B = g + Wand B corresponds to the llnear varlety 1n the propos1t1on. 

Hence g* 1S the unlque element 1n W ~ sat1sfY1ng Tg* = d. 

We conSlder the contents of W ~. Any element wh1ch 1S orthogonal 

to each N 1S also orthogonal to dny llnear comb1natlon of the B1 

~ splJ.nes. Hence S: = span(Nl'N2

, ••• ,Nm)lSa subset of W • We now 

show that W ~ 1S a subset of S (and hence S = W ~) by contrad1ctlon. 

~ Assume that there eX1sts an element y wh1ch 1S In W but not III S. 

Slnce S 1S a closed subspace there eX1sts (by the PrO]ect1on Theorem) 
an element s 1n S such that 

o 

II y - s II ~ II y - s II 

o 

for all s E: S 

w1th y - s 1n the orthogonal comp]ement of S. Th1S 1mpl1es o 

T(y - s ) = 8 or (y - s ) E: W. However y - s lS also 1n W.L Slnce 

o 0 0 

both y and s are 1n W.L. Therefore (y - s ) = 8 and S = W.L. 

o 0 

In summary, g* lS character1zed by 

m 

LaN 

1=1 1 1 

(Slnce gi~ lS 1n the span of the B-spl1nes) where the coeff1c1ents 

m 

9 

( LaN ,N ) 

J=l J J 1 

d 

1 1=1,2, •.• ,m (1.12) 

(slnce Tg* = i). Equat10n (1.13), a system of m llnear equat10ns 1n m 

unknowns, can be wr1tten 1n matr1x notat1on as 

A a = d (1.13) 

where the symmetr1c matr1x A has entr1es A = (N ,N ). 1J 1 J 

Because the B-spl1nes are llnearly 1ndependent, the matr1x A, a 

Grahm matr1x, lS nons1ngular and hence a un1que Solut1on eX1sts for 

any glven .Q.. Furthermore, Slnce N has support [t,t 2]' the matr1x 1 1 1+ 

A 1S tr1d1agonal. m For any ~ E R we have 
m T 

x Ax = r x (Ax) 

= 

1 - 1 1=1 

m m 

r X (N , r X N ) 

1=1 1 1 J=l J J 

m m 

= ( r X N , r X N ) 

1=1 1 1 J=l J J 

m 2 = II r X N II 2 

1=1 1 1 

~ 0 

w1th equa11ty hold1ng 1f and only 1f X = 8. The matr1x A lS hence 

10 

pos1t1ve def1n1te and (1.13) can be solved by Gauss e11m1nat1on w1thout p1vot1ng, or, better st11l, by Cho1esky decompos1t1on. 

We note also that 

The entry A , the 1ntegra1 of the product of two p1eceW1se 11near 1J 

polynom1als, can be computed exactly by Slmpson's rule app11ed on each 

sub1nterval [tk ,tk+1 ]. Denot1ng 6tk : = t k+l - tk and zk the m1dpo1nt 

of the 1nterval [tk ,tk

+1 ] we have for 1=1,2, ... ,m 

A 11 

= (6t 1/6)[N (t )2 + 4N (z )2 + N (t 1)2] 1+ 1 1 1 1 1 1+ 

+ (6 t 16)[ N (t ) 2 + 4N (z ) 2 + N (t ) 2 ] 1+2 1 1+1 L 1+1 1 1+2 

(t 2 - t )/3. 1+ 1 
11 

We also compute for 1=1,2, •.. ,m-l 

A = A 1,1+1 1+1,1 

I

t 1+2 

N N l(t)dt t 1 1+ 

1+1 

The solutIon gi~' beIng a lInear combInatIon of lInear B-sp11nes, 

IS pIeceWIse lInear (and contInuous) wIth knots t. After IntegratIng 1 

g* tWIce and applYIng the InterpolatIon condItIons, we obtaIn f* whIch 

IS piecewIse CUbIC (wIth knots t ) wIth two contInuous derIvatIves. 

1 

1=1 

1=2,3, .•. , n-1 

l=n 

and 68 = Bl+1 - 8

1 • On [t1 ,t1+1 J f~ IS defIned by a unIque CUbIC 

polynomIal p~~ and hence f* can be determIned by specIfYIng the numbers 

Pu (J)(t ) for 1=1,2, •.. ,n-1 and J=0,1,2,3. Then 

><1 1 

+ (1.14) 

2' 3' 

for t£ [t t 1]. Of course, (1.14) can be more effIcIently evaluated 

1, 1+ 

by uSIng nested mu1t1p11catlon. 
The po1ynomla1 p* solves the dlfferentla1 equatlon 1 

61 + (~6 /~t )(t-t ) 1 1 1 

12 

(1.15 ) 

on the lnterva1 [t ,t 1] wlth boundary condltlons p* (t ) = y and 1 1+ 1 1 1 

P• (t ) = y . Therefore '<1 1+1 1+1 

0.16) 

where the constants c and e are evaluated as e = y and 1 1 1 1 

c = ~ Y 1 _ (61+ I + ~s 1) ~ t 

1 M 2 6 1 

1 

(1.17) 

P (O)(t) = y *1 1 1 

p. (I)(t ) = c 

*1 1 1 (1.18 ) 

~B /~t 1 ] 

where c lS glven by (1.17). A complete FORTRAN program for computlng 1 

the natural CUblC sp11ne lnterpo1ant lS glven In Appendlx A. 

Flgure (1.3) dlsplays the natural CUblC spilne lnterpolant that 

lS 10 cootrast to the po1yoomla1 10terpo1ant of flgure (1.1). 

We complete thlS chapter by poslng (and solvlng) a generallzdtlon 

of problem (1.1). For k flxed satlsfYlng 2 ~ k ~ n, let A(k) be the 

(k) be the set of 10terpo1ants (to the data) WhlCh are 10 L2 [d,b]. We 
13 

F1gure (1.3): The Natural CUb1C Sp11ne Tnteroolant. 
14 

cons1der the problem 

F1nd f~~ E A(k) such that II f{~(k) "2 ~ "f(k)" 2 for all f E A(k) 

(1.19) 

Let f be an element of A(k). S1nce (1.3) 1S va11d for f, we can 

1ntegrate by parts aga1n (assum1ng k > 2) to obta1n 

where 

f(t) 

(2) 

= f(a) + f(l)(a)(t-a) + f (a) (t_a)2. 2' 

In genera~ after 1ntegrat1ng by parts k-1 t1mes we obta1n 

l

b _ (t_s)k-1 (k) 

f(t) - qk-1 (t) + (k-1)' f (s)ds 

a 

or 

}

\ t_s):-l (k) 

f(t) = qk-l (t) + (k-l)' f (s)ds. 

a 

(1. 20) 

(1. 21) 

(1.22) 

Now we take the (scaled) k-th d1v1ded d1fference of (1.22) to 

obta1n 

b 

d, •

k = I g(s)N, •

k (S)ds 

a 

1=1,2, ... ,n-l (1. 23) 
15 

d l,k Ck-l)'(t k-t )[t t k]fCo) , 1+ 1 1, ... , 1+ (1. 24) 

N k(s) 1. , 

k-l (t -t )[t t ](o-s) l+k 1 1, ... , l+k + (1. 25) 

(the normallzed B-spl1ne of order k), and g = f(2). 

Let B(k) denote the set of elements (In L2[a,b) wlch satlsfy 

(1.23). Then the Solutlon fi~ to (1.20) lS related to the Solutlon 

to the problem 

Flnd g-l~ E B(k) such that II g-l~ (k) II 2 ~ II g (k) 112 for all g E B(k) (1.26) 

Vla gi:- fu(k). Furthermore, for some a E Rn- k 

we have "i\ 

n-k 

L: aN k' 

J=l J J, 

The coefflclents a l,a2

, ... ,an_

k are chosen to solve the llnear system 

of n-k equatlons In n-k unknowns represented by the matrlx equatlon 

A a = d where A lS symmetrlc and posltlve deflnlte wlth entrles 

A (N k,N k)' lJ 1, J, 

Slnce g-l~ lS a llnear comblnatlon of plecewlse polynomlals of order 

k, fA wlll be a plecewlse polynomlal of order 2k. We call fA the 

natural spllne lnterpolant of order 2k. 
2. A M1n1mal Norm Interpolatlon Problem 

1n the L [a,b] Spaces p 

For p such that I < p ~ 00 we def1ne the set 

G(p): = [g £ Lp[a.bJ: rg(t)~,(t)dt 

for ,:1.2 ..... n} 

= fb g (t) <l> (t) d t 

o 1 

a 

(2.1) 

where {<l>1}:=1 1S a set of elements 1n Lq[a,bJ, q 1S conjugate to p 

(p+q = pq 1f p~ 00 and q=l 1f P = 00), and g 1S a flxed element of 

o 

L [a,b]. Cons1der the problem p 

F1nd g* £ G(p) such that /lg* /I p ~ /I g /I p for all g £ G(p). (2.2) 

In chapter I we solved (2.2) for the spec1al case p=2; f1nd1ng 

from a 11near var1ety 1n a H1lbert space the element of m1n1ma1 norm. 

The ProJect1on Theorem came 1n handy to character1ze g* as well as to 

guarantee un1queness. However, for p ~2 L [a,b] does not have the 

p 

orthogonal1ty propert1es of a H1lbert space and hence, we cannot use 

the ProJect10n Theorem to solve (2.2). Instead we solve (2.2) 1n th1S 

chapter by ut1l1z1ng the Hahn-Banach theorem to character1ze gi< 

Unlqueness follows In the case I < p < 00 by the strlct convexlty of the 

norm. ThlS chapter, modeled after [deB(2)], motlvates the use of the 

Hahn-Banach theorem In chapter 5. 

Let A be the 11near functlonal deflned on the subspace 

16 
17 

v1a 

r ,il"1·1) (t )go (t)d t. 

a 

(2.3) 

Any element of G(p) (lnclud1ng g ) w1ll serve as an extens10n of ~ to 

o 

a bounded 11near funct10nal def1ned on all of L [a,b]. Hence, q 

III. lis ~ IIglip for all gEG(p). (2.4) 

Conversely, any extens10n of ~ to a bounded Ilnear funct10nal def1ned 

on all of L [a,b), be1ng 1dent1cal to A on S, 1S represented by an 

q 

element of G(p). 

The Hahn-Banach theorem guarantees the eX1stence of an element 

g E G(p) such that 

b I f(t)g(t)dt ~ a 

II A II • II f I / for all f E L [a, b J • 

s q q 

Th1s 1mplles thalli g II ~ II), II Wh1Ch, taken along W1 th (2.4), glVes us 

s 

118 II = II A /I and, therefore, a Solut1on to (2.2). Now we charactcrve 

s 

g. 

n i~ 

Let [a ¢ be an element such that 

1=1 1 1 

n n 

II [ a ¢ II = 1 

1=] 1 1 q 

and I.e [ a ¢ ) = II AI' . 1=1 1 1 S 

(Th1S ele~ent lS un1que 1f 1 < p < 00 Slnce the norm lS str1ctly 

convex.) Then 
n 

= A( Lac!> ) 1=1 1 1 

j

b n 

( r a c!> )(t)g(t)dt 1=1 1 1 

a 

n 

~ 1/ Lac!> I/·I/sl/ n=l 1 1 P 

= IIsll . p 

Therefore, equa11ty holds throughout and we have 

b 

I 

( ~ a c!> )(t)g(t)dt 1=1 1 1 

a 

n A 

Slnce g and Lac!> are a11gned, we must have 

1=1 1 1 

n n 

= II All • 1/ Lac!>" q-1 slgnum ( Lac!> )( t). 

s 1=1 1 1 1=1 1 1 

18 

We close th1S chapter by stat1ng the 1nterpo1at10n problem that 

goes along w1th solv1ng (2.2). Let p be a number such that 1 < p < 00 , 

let k be an 1nteger such that k ~ 2, and let f E L (k)[a,b]. Def1ne 

o p 

the sets 

F: {f £ L (k)[a,b]: f(t ) = f (t ) P 1. 0 1 

1==1,2, ... ,n} 
and 

b 

G: = {g E Lp[a,b]: 1 

a 

Then the prob1em5 

and 

g(t)N k(t)dt = 

1, 

d 1,k 1=1,2, .•. ,n-k} 

~ II [(k) II for all f E F 

p 

F1nd f~~ E G such that II g~~ II p ~ II gil p for all g E G 

are equ1valent and 

19 
3. The Convex Spllne Interpolant 

n The data {(tl,yl )} 1=1 are called convex If the pOlnt (t1 'Y1 ) 

2 2 

lles on or beneath the llne JOlnlng the pOlnts (t ,y ) and 

11 11 

(tl 'Yl ) whenever 1 ~ 11 < 12 < 13 ~ n. Equlvalently, 

3 3 

(where f lS any lnterpolant to the data) or 

d 

1 

for 1 = 1,2, .•• ,m(= n-2). 

Yl+2 - YHI 

tl+2 tl+l 

In thlS chapter we address the problem of flndlng, for convex 

data, the smoothest convex lnterpolant; that lS, the convex lnterpolant 

havlng second derlvatlve of mlnlmal norm over all smooth convex lnterpolants. The natural CUblC sp11ne lnterpolant. the smoothest of all 

lnterpolants, regrettably does not always preserve the convexlty of 

the data. In chapter 1 we showed that f*, the natural CUblC spllne 

lnterpolant, has second derlvatlvP 

m 

LaN 

J=l J J 

where the coefflclents ~,a2, .•. ,am satlsfy (1.13). If any a lS 

1 

negatlve, then f* lS actually concave on a subset of [a,b]. 

Let {(tl'Yl)}~=l denote convex data and let A denote the set 

(2) of convex lnterpolant~ ln L2 [a,b]. We assume that A lS nonempty. 

20 
21 

(There are convex data for WhICh A IS empty. For example, let 

Yl = Itll and tl = -2, t2 = -1, t3 = 0, t4 

convex Interpolant IS f(x) = lxi, WhICh IS 

= 1, and ts = 2. The only 

(2) not In L2 [-2,2].) 

USIng the Peano kernel theorem as we dId In chapter 1 we can show 

that If f £ A then T(f(2» = d where T:L2[a,b] +Rm IS gIven by 

(Tg ): = (g, N ). Hence 1 f 

1 1 

B { g E L2 [ a , b ] : g;;; 0 and T g = QJ, 

then problems 

FInd f* £ A such that IIf{~(2) 112 ;;; II f(2) 112 for all f £ A (3.1) 

and 

(3.2) 

are equIvalent and the SolutIons are related VIa gi~ = f*(2). SInce B 

IS a nonempty closed convex set, we consIder (3.2) as fIndIng the 

dIstance from a pOInt to a closed convex set In a HIlbert space. 

ProposItIon ([L, page 69]): Let x be an element of a HIlbert space H 

and let K be a nonempty closed convex subset of H. Then there eXIsts 

a umgue element k E K such that 

o 

II x - k II ;;; II x - k II for all k E K 

o 

Furthermore, k IS characterIzed by o 

(x - k k - k ) ::; 0 for all k E K. 

0' 0 

SInce we wIsh to fInd the element of mInImal norm In B, X corrcsponds to the zero functIon and hence g* IS characterIzed by 

(3.3) 
22 

Propost1on ([MSSW, propos1t1on 2.1]): If there eX1st coeff1c1ents 

b 

J ( ~ <l N ) N (t)dt = d 1=1,2, .•. ,m, (3.4) 

a 

J=1 J J + 1 1 

/II 

then g~~ = ( l. (l N ) • 

J=1 J J + 

Furthermore, such coeff1c1ents eX1st If there 

m 

eX1sts g e: B such that {N1}1=1 are 11near1y 1ndependent over the support 

of g. 

Proof: Assume <ll,<l2, .•• ,<lm sat1sfy (3.4). 

assume g e: B. Dehne (h) = (-h) + so that 

h= (h) - (h) . + -

Then 

«s)+' g-(s)+) 

(s + (s) , g - (s)+) 

m 

Denote s = L (l Nand 

J=1 J J 

= (s,g - (s)+) + «s) ,g - (s)+) 

= «s) ,g) - «s) , (s)+) 

= «s) ,g) 

~ O. 

The last 1nequal1ty lS val1d Slnce both (s) and g are nonnegat1ve 

funct10ns. Hence (s)+ sat1sf1es (3.3). 

We now show that we can f1nd coeff1C1ents (l1,(l2, ... ,(lm so that 

(3.4) holds by f011ow1ng the procedure employed 1n [MSSW]. 
23 

We beglo by cooslderlng the problem 

f

b m 

lnf (L:a N )2 (t)dt : L:a d 

J=lJ J + 1 1 

a 

1 (3.5) 

and showlng that lf the lnflmum lS attalned at some~, then for some 

posltlve constant C the coefflclcnts Cal' Ca2, ... , Cam satlsfy (3.4). 

-3~ -3" 

If the lnflmum of (3.5) lS attalned at Q, then a lS a crltlcal 

pOlnt of the Largranglan 

b 

L~, A) = f ( ~ aN) (t) d t 

J=l J J + 

a 

At a mlnl~um of L we must have 

b 

o = 2 I ( ~ aN) N (t)dt - J=l J J + 1 

a 

and a • d = 1 for some A. 

m 

+ A( 1- L: ad). 

J=l J J 

Ad 

1 1=1,2, .•• ,m 

Now mu1t1p1y (3.7) by a and sum over 1=1,2, ..• ,m to obtaln 

1 

J

b m m 

2 (L: aN) ( L: aN) (t)dt - J J + 1 1 J=l 1 =1 

a 

or 

m 

A L: a d 

1=1 1 1. 

b 

A = 2J ( ~ a N );Ct)dt ~ o. 

J=l J J 

a 

If A>O, then (3.7) reveals that 

o 

(3.6) 

(3.7) 

(3.8) 
24 

J

b m * 

( LaN) N (t)dt = d 

J=l J J + 1 1 1=1,2, ... ,m (3.9) 

a 

where a * = 2 a IA. If A = 0, then (3.8) reveals that 

J J 

where ~ • d = 1. 

g E B we have 

J

b ( E aN) + ( t) d t : 0 

J=l J J 

a 

m 

Th1S lmp11es that ( L a.N ) ~ o. 

J=l J J 

J

b ( ~ a N )g(t)dt 

J=l J J 

a 

m 

L a (N ,g) 

J=l J J 

m 

Lad 

J=l J J 

= 1 

However, for any 

Wh1Ch lS 1mposs1b1e because g 1S nonnegat1ve on [a,b]. We conclude 

that A lS str1ct1y pos1t1ve and, 1f the 1nf1mum 1n (3.5) lS atta1ned 

by some~, that (3.4) lS solvable. We now show that the 1nf1mum lS 

atta1ned. 

Let {ik)}~=l be a m101m1z1ng sequence. If {IIQ:(k) 1/ }~=l lS 

unbounded, then d1V1de the obJect1ve functlOn of (3.6) by II~ 1/: and 

the constra1nt by I/a 1/ • There then eX1sts a such that - 00 
/I ~ II <Xl = 1, 

a·£!. = 0, and 

m 

We conclude that ~ a N 1S nonpos1t1ve, but not 1dent1cally 

J=l J J 

zero. Slnce we have assumed there eX1sts gEB such that the 

B-spl1nes are llnearly 1ndependent on the support of g, 

m m 

o = ~ a d 

J=l J J 

L: a (g,N ) 

J=l J J 

m 

= (g, L: aN) 

J=l J J 

< 0 

25 

Wh1Ch 1S a contrad1ct1on. Hence a m1n1m1z1ng sequence must be bounded 

and the 1nf1mum 1S atta1ned V1a a convergent subsequence. Th1S 

completes the proof of the propos1t1on. 

We note that the eX1stence of g E B, such that {N } mare llnearly 

1 1=1 

" 1ndependent over the support of g, 1n the prev10us propos1t1on 1S 

guaranteed 1f d > 0 for each 1. Then each g E B must be pos1t1ve on 

1 

some sub1nterval of [t,t 2]' the support of N , for each 1. 

1 1+ 1 

Now we cons1der the 1mpl1cat1on of allow1ng d

k = 0 for some k. 

As a spec1f1c example let t 

1 

T (1-1) for 1=1,2,3,4 and let £!. = (1,0) • 

If g* 1S the pos1t1ve part of a llnear comb1nat1on of B-spl1nes, then 

there must eX1st numbers a l and 0:2 sat1sfY1ng 
26 

1=1 

(3.10) 

Wh1Ch 1.mplles that a 2 = - 00 Th1S lS equ1va1ent to the Solut10n 

be1.ng 1dent1.cally zero on [1,3]. In fact, any g € B must be of the 

form g = g X[O,l]' It lS shown 1n [MSSW, theorem 3.1] that the 

So1ut10n to (3.2) lS 

r: 

m 

g = ( raN) Xr * J=1 J J + 

m 

= [a,b]/{U 

J=1 

(t ,t 2) : d J J+ J O} • 

Hence the So1ut10n to (3.2) w1th t = (1-1) for 1=1,2,3,4 and 

1 

d = (I,O)T lS 

Unless otherWlse stated we assume d > ° for each 1 for the remalnder 

1 

of th1S chapter. 

Before we cons1der how Lo compute the coeff1clents a r a 2 , ... ,am 

Wh1Ch sat1sfy (3.4), we glve a procedure for 1ntegrat1ng g*. 

Detlne S , 6S , 6t , and 6y as 1n chapter 1. We 1ntegrate g~~ on each 1 1 1 1 

sub1nterval ft,t 1] separately, form1ng a p1.eCeW1Se polynom1al, by 1 1+ 

solv1ng the d1fferent1al equat10n 

p (2)(t) = (S 

*1 1 

6B 

+ _l(t_t » 

lit 1. + 

1 

for t 5 t ~ t w1th boundary cond1t10ns p*(t1 ) = Y1 and 1 1+1 

(3.11) 
27 

Two 1ntegrat10ns glves us 

(3.12) 

and 

llt 

P ( t ) _---::;.1 ~ (6 + 66 (t - t » \ c ( t - t ) + e 

*1 = 6(66)2 1 1 1 + 1 1 1 

1 

(3.13) 

for constants c and e. We proceed by cases. 

1 1 

Case] occurs when both 61 and 61+1 are nonnegat1ve. The nonnegat1V1ty constra1nt 1S not act1ve 1n th1S case and so (3.13) 1S equ1valent to (1.16), although w1th ffiod1f1ed constants c and e. The values 

1 1 

P*1(J)Ct1) for J = 0,1,2,3, are glven by (1.18). 

Case 2 occurs when 6 < 0 and S 1 > O. In th1S case p~ can be 1 1+ ~1 

def1ned by two po1ynoffi1a1s: a 11near po1ynoffi1a1 q11 def1ned on 

[t ,T ] - where the nonnegat1v1ty constra1nt 1S act1ve and hence the 

1 1 

second der1vat1ve 1S zero - and a CUb1C po1ynoffi1a1 def1ned on 

[T ,t 1] where 1 1+ 

T 

1 

t - 111 6 6t /66 1 

App1Y1ng the boundary cond1t10n P*1(t1 ) = Y1 we obta1n e1 

App1Y1ng P*1(t1+1 ) = Y1+1 we get an equat10n for c1 : 

Solv1ng for c we have 

1 

+ c 6t 

1 1 

+ y 1 

(3.14) 

y • 1 
From (3.11), (3.12), and (3.13) we obta1n 

C (1 -t ) + y 1 1 1 1 

C 

1 

where 1 and care glven by (3.14) and (3.15) respect1vely. 1 1 

28 

(3.15) 

(3.16) 

Case 3 occurs when 81 > 0 and 81+1 < O. In th1S case P*l 1S 

def1ned by a CUb1C po1ynom1a1 q 1 on [t ,1 ] and by a 11near po1y- 1 1 1 

nom1a1 q12 on [1,t 1] w1th 1 def1ned by (3.14). These po1ynom1a1s 1 1+ 1 

are determ1ned by the values 
29 

q (2)(t ) = B 11 l 1 (3.17) 

C (1 -t ) + e 111 1 

C. 

1 

q(2)(1 ) = 0 

12 1 

where c and e are glven by 1 1 

and 

C 

1 

e 

1 y - 1 

«(3 )3Llt 

1 1. 

2(Ll(3 )2 • 

1 

6(Ll13 ) 2 

1 

Case 4 occurs when (3 and (3 1 are both nonpos1t1ve. In th1S 1 1+ 

case we obta1n a 1Jnear po1ynomla1 def1ned on [t ,t 1] and determlned by 1 1+ 
30 

P''!-l (t1) = Y1 

p(l)(t ) *1 1 6y /6t 1 1 

(3.18) 

P;~)(t1) = 0 

~3)(t ) = O. ><1 1 

Slnce g* lS p1eceW1se 11near and cont1nuous (w1th knots at the 

t 's and T IS), f* w111 be p1ecew1se CUb1C w1th two cont1nuous 

1 1 

der1vatlves (If d > 0 for each 1). We call f{~ the convex CUblC 

1 

sp11ne 1nterpo1ant. 

Now we turn our attent10n to the task of numer1ca11y ca1cu1at1ng 

the coeff1clents Ci1 ,Ci2 , ••• ,Cim WhlCh sat1sfy (3.4). We contlnue to 

assume that d > 0 for each 1. 

1 

where 

F (a) 

1 -

b 

= I( ~ aN) N (t)dt J=l J J + 1 

a 

We w1sh to solve F(~) = d. 

1=1,2, ..• ,me (3.19) 

One method lS to use Jacob1 1terat10n. An 1nlt1a1 guess 

_ (0) (0) (0) (0) T (k) 00 

x = (xl ,x2 , ... , xm ) lS chosen and a sequence {x } - k=o 

lS generated by ca1cu1at1ng ~(k+l), once ~(k) 1S known, by solvlng 

(k) (k) (k+1) (k) (k) F (x 1 , ••. ,x l' x , x +1 ' •.• ,x ) d 1 1- 11m 1 
31 

(k+1) for Xl for each 1. A mod1flcat1on, the Gauss-Se1del 1terat1on, 

(k+l) (k) 1nvolves ca1cu1at1ng ~ , once ~ lS known, by solvlng 

(k+1) for X for 1=1,2, ... ,m 1n succeSS1on. Both Jacob1 and Gauss1 

Sledel 1terat10ns converge globally a5 proved 1n [IMS] 

Now we cons1der Newton's method to solve G(~) = F(~) - d = 8. 

(0) (k \ <Xl We p1ck a sU1table 1n1t1al guess 2£. and form a sequence {~ '} k=o 

by solv1ng 

(3.20) 

[or X(k+l) once x(k) lS \<nown. S 1nce nG nF (3 20) v = v , we can express . 

alternately a5 

(3.21) 

The entr1es of the Jaboc1an matr1x 'VF are 

('V F) (a) = J\ ~ akNk f N N (t) d t 

1J - k=l + J 1 

a 

(3.22) 

m 

where ( L: \:Nk)o 1S the characterlstlc funct10n for the support of 

k=1 + 

m 

( L: aN) . We see that 'VF 1S symmetr1c and tr1d1agonal at each~. 

J=1 J J + 

We now character1ze those a for whLch ('VF)~) lS pos1t1ve def1n1te. 

Lemma (3.1): The Jacob1an ('Vf)(~) 18 pos1t1ve def1n1te 1f and only 1f 
m 

( r akNk) does not van1sL ldent1ca11y on any of the sub1nterva1s 

k=l + 

[t,t 2] for 1=1,2, •.• ,m. 1 1+ 

Proof: m For any ~ E K we have 

T ~ (VF)(~)~ = 

m 

L: 

1=1 

:1: 0 

m 

x L: (VF) (a)x 1 J=l 1J - J 

If ( L: aN) does not van1sh 1dent1cally on [t,t 2] for each 1, J=l J J + 1 1+ 

32 

then equal1ty holds 1f and only 1f x = 0 for each 1. If there eX1sts 

1 

m 

some k such that ( L: aN) 1S 1dent1cally zero on [tk,tk+2], then 

J=l J ] + 

equal1ty does hold for the nonzero vector ~ def1ned by Xl elk 

for each 1. Th1S completes the proof of the lemma. 

FrolT, (3.20) we see that 

m 

F ~) = L: a r ( m L: ~~ )~NJN1 (t)dt 1. J=l J k=1 

a 

m 

r a (yF) (a) 

J=1J 1.1-
33 

so thaL F(cd (Vf)(a) a. Newton's method - equatIon (3.22) - takes 

the form 

(3.23) 

Th eorem (3 ._ ?) : If ( II F')( ~ (k» IS posItIve d f h e lnlte, t en (nF)(~(k+l) v 

IS posItIve defInIte for each k and, hence, Newton's method - equatIon 

3.23) - IS always well-defIned. 

Proof: HavIng the known values (k) we wIsh to determIne the values Xl ' 

(k+l) X satIsfYIng 1 

J ( m x (k+l)N )N (t)dt 

S(k) J!l J J 1 

d 

1 1=1,2, ... ,m 

where S(k) IS the support of ( ~ x (k)N ). Slnce (VF)(~(k» 1S 

J=l J J + 

(3.24) 

posltlve def1nlte, then S(k)U[t,t 2] conta1ns an lnterval for each 1. 

1 1+ 

Slnce d > 0, then 

1 ( ~ x (k+1)N )+ b 1 f L. 1S POSI t1 ve on some su 1nterva 0 

y=l J J 

mF)(~(k+1» d f Th 1 \v lS posltlve e 1n1te. lS comp etes 

the proof of the Theorem. 

Note that lf x(O) has all posLtlve components (for example, lf 

x (0) = 1 for each 1, then S(O) = [a,bl and ~ x (l)N IS the second 

1 J=l J J 

derlvatlve of the naLural CUblC spllne Interpolant. 

Now we assume that dk = 0 for some k. In thIS case specIal care 

must be exercIsed SInce {xk (J)}~=O may dIverge to _co, preventlng any 

numerlcal convergence. We already know that dk = 0 ImplIes that the 
34 

dala pOlnts (tk,yk ), (tk+l,yk+l ), and (lk+2'Yk+2) are colllnear and, 

hence, any convex lntcrpolant must be llnear on [tk ,lk+2]' Equlvalently, the second derlvatlve of any convex Lnterpolant must be zero 

m 

( L x N) {X + X } J=l J J + [a,tk ] [tk+2 ,b] 

Slnce the value of x

k lS lmmaterlal and the k-th equatlon lS 

automatlcally satlsfled, the number of equatlons and unknowns each 

reduce by one. For computatlonal convenience (3.23) can stlll be used 

wlth the followlng modlflcatlons: (VF)kk = 1, (VF)k,k+l = 0, and 

(VF\,k_l = O. 

If d

k = 0, then the solut1on lS d1scontlnuous at tk lf x

k_l > 0 

and 18 d1scont1nuous at t k+2 1f Xk+l > 0. If the Solutlon lS dlscontlnuous, then f* wlil have only one contlnuous derlvatlve. 

A further problem lS encountered when dk_l and dk+l are both 

zero, but dk lS nonzero for some k. Any nonnegatlve functlon g WhlCh 

sntlsf1es the (k-l)-st and (k+l)-st equat10ns can not sat18fy the k-th 

We conclude that there does not eXlst any convex 1nterpolant In 

L

2

(2)[a,b] (and no Solutlon to the problem as posed). However, we can 

flnd a convex lnterpolant whose second derlvatlve 15 of the form 
3S 

m 

( L x N ) { x + x } J=l J J + [a,tk_l ] [tk+3,b] 

sat1sfY1ng all but the k-th equat1on. We already know that th1S convex 

hence, p1ecew1se llne8r on [tk_l ,tk

+3 ]. If d

k 1S nonzero, then there 

w1ll be a d1scont1nu1ty 1n slope at t

k

+l . For the conven1ence of 

utll1zeng (3.23) we can set d

k to be 7ero to sat1sfy the k-th equat1on. 

The dlscont]nu]ty 1n slope wlll show up after we lntegrate the Solutlon 

Lo obtaln the Interpolant. 

Flgure (3.1) dIsplays the natural CUblC spllne lnterpolant to the 

functlon 

f(t) 1 

(O.OS+t)(l.OS-t) 

at the knots tl = 0, t2 = 0.1, t3 = 0.4, t4 = 0.7, ts = 0.8, and 

t6 = 1.0. Flgure (3.2) dlsplays the convex spllne lnterpolant to thlS 

functlon. Table (3.1) shows the convergence results for Jacobe, GaussSeldel, and Newton's method lteratlons taken from [IMS]. Note the 

quadratlc convergence character1st1C of Newton's method. These convergence results are tYPlcal. 
Flgure (3.1): The Natural CUblC Spllne Inter~olant. 

'l", .uw'iF""':l:.u.w~"""'1 .. ~~·.u".w!l'li'~'.w"·.w .. '-'l;'~.w,, '.LU" '~!i~, ..... IS 

F1 :;:Jre (3.2): The Convex (uhle SphIW l'1terpola'1t 
37 

TABLE 3.1 

IIF(x(n») - d 112 

Iteratlon Gauss 

Number Jacobl Selde1 Newton 

1 .46 x 102 .27 x 102 .19 x 102 

2 .28 x 102 .11 x 102 .85 x 601 

3 .75 x 101 .42 x 101 ,29 x 101 

4 .12 x 102 .18 x 101 .49 x 100 

5 .26 x 101 .75 x 100 .14 x 10-1 

6 .49 x 101 .31 x 100 .11 x 10-4 

7 .10 x 101 .13 x 100 .71 x 10-11 

8 .21 x 101 .55 x 10-1 .49 x 10-12 

9 .43 x 100 .23 x 10-1 

10 .86 x 100 .96 x 10-2 

20 .11 x 10-1 .16 x 10-5 

30 .14 x 10-3 .26 x 10-9 

40 .18 x 10-5 .58 x 10-l3 

50 .24 x 10-7 

60 .30 x 10-9 

70 .39 x 10-11 
4. The Shape-PreserYlng Spl1ne Interpolant 

We addres5ed ln chdpter 3 the problem of flndlng, for convex data, 

the smoothest convex lnterpo]ant. We begln thlS chapter by conslderlng 

the problem of flndlng, for (oncave data, the smoothest concave lnterpolant. Then we contlnue the chapter by examlnlng the problem of 

flndlng, for general data, the smoothest lnterpolant Whlch IS locally 

convex where the data are locally convex and IS locally concave where 

the data are locally concave. 

Let {(t ,y )}n I denote concave data and let A denote the set of 1 1 1= 

(2) all concave Interpolants In L2 [a,b]. Assume A IS nonempty. USIng 

the Peano kernel theorem as we dId In chapter ], we see that, If f £ A, 

then 

b 

J f(2)(t)N,(t)dt = d, 

a 

1=1,2, ... ,m(=n-2) 

Equlvalently, we have T(f(2~ = d. 

DefInIng 

g S 0 and Tg = &} 

we conclude that the problems 

FInd f.::_ £ A such that IIr/ 2)" 2 s /If(2) /I 2 for all f £ A (4.1) 

(the problem of fIndIng the smoothest concave lnterpolant) and 

FInd g-l:- £ B such that II g,:- II 2 ~ II g /I 2 for all g , B 

38 
(2) 

are equIvalent and the solutIons are related VIa g* = f* . 

Of course, the smoothest concave Interpolant to the concave data 

{(tl'Yl)}~=l IS the negatIve of the smoothest convex Interpolant to 

the convex data {(t ,-y »n]. We hIghlIght thIS WIth the followIng 1 1 1= 

proposlt1.0n. 

PropOSItIon [MSSW]: If there eXIst coeffICIents ul ,U2 'U3 ' ... U

m 

sat Isfnng 

f-a 

m 

m 

( LaN) N (t)Jt 

J=l J J - 1 

d 

1 1=1,2, .•. ,m (4.3) 

then g~ = - ( L uN). Furthermore, such coeff1clents eXIst If there -- -" ]=1 J ] - -

eXlf3ts g E B such that {N1}:=1 are lInearly Independent over the 

support of g. 

We note that the eX1stence of g E B, such that {N1}:=1 are 

llnearly Independent over the support of g, In the preVIOUS propOSItIon 

IS guaranteed If d < 0 for each 1. Then each g E B 1S negat1ve on 

1 

some sub1nterval of [t ,t 2 J, the support for N , for each 1. 

1 1+ 1 

Now we conSIder the probJem of f1nd1ng, for general data, a 

smooth shape-preserv1ng 1nterpolant - a smooth Interpolant WhICh IS 

locally convex where the data are locally convex and IS locally concave 

where the datd are locally concave. AssumIng for the moment that d 

1 

IS nonzero for each 1, we defIne the sets 
and 

T 1 

Now we def1ne the sets 

d > O} , 1 

d > O} , 1 

and f(t) = y 1=1,2, •.. ,n} 1 1 

(whlch we assume l~ nonempty) and 

B: gX n = , > 0 gX n < 0 and Tg = _d} = , 

1 2 

We conclude that the problems 

40 

F1nd f.;:_ S A such that IIf/ 2) II 2 ~ II f(2) 112 for all f E A (4.4) 

and 
41 

(2) 

arc equIvalent dnd g* : i* . 

The followlng proposltlon glves the Solutlon to (4.5). We see 

that f,cX~ has the character of the convex sp11ne lnterpolant, f~cX~ 

1 2 

has the character of the concave spllne lnterpolant, and f*X~ has the 

3 

character of the natural &pllne lnterpo1ant. 

Proposltlon (MSSW]: If there eXlsts coefhcents ~,a2, ... ,o.m 

satlsfYlng 

m 

+ ( L: Ct N )X~} N (t) d t 

J:1 J J 3 1 

then 

d 1:1,2, ••. ,m 1 

(4.6) 

A Furthermore, such coefflclents eXlst lf there eXlsts g E B such that 

m A { N~ 1=1 are Ilnearly lndependent over the support of g. 

We note that thc eXlstence of ~ E B, such that {N}:=l are 

IJnearly lndependent over the support of g, In the preVlOUS proposltlon 

lS guaranteed 1f d l~ non7cro for each 1. Then each gE B lS non7cro on I 
42 

some sub1nterval of [t,t 2]' the support of N , for each 1. 1 1+ 1 

We now 501ve (4.6). Def1ne F : Rm + R

m where F = (F ,F2, ... ,F )T 1 m 

F (x) 

1 -

m 

( L x N ) N (t)dt J=l J J + 1 

m 

( L x N ) N (t)dt J=l J J - 1 

+ 1 (~x N )N (t)dt S""2 J=l J J 1 

3 

1=1,2, ... ,m (4.7) 

We use Newton's method to solve F~) = £. P1ck1ng a sU1table 

1n1t1al guess ~(O) we produce a sequence {~(O), ~(l), ...• } by solv1ng 

f (k+1) (k) or x once x lS known. The Jacob1an matr1x has entr1es glven 

by 

b 

(VF) (a) = f P(a)N (t)N (t)dt 1J - - J 1 

a 

(4.9) 

where 
43 

m m 

P(~) o 0 ( LaN ) X~ + ( LaN ) X~ + X~ . J=l J J + 1 J=l J J - 2 3 

(4.10) 

From (4.9) we see thatVF 1S symmetnc and tr1d1agona1 at eacha. 

We also note that 

m 

P(x)( L x N ) 

- ]'=1 J J 

m 

( L x N ) X~ 

J=l J J + 1 

m 

- ( L x N ) X~ 

J=l J J - 2 

m 

+ ( L x N )X~ 

J=l J J 3 

so that F(~) = (VF)(~)~ and, hence, (4.8) reduces to 

(4.11) 

The fo11owlng lemma (wlth proof Slml1ar to ltS counterpart ln 

chapter ~ characterlzes those a for WhlCh (VF)(£) lS posltlve deflnlte. 

Lemma(4.1): The Jacoblan (VF)(£) lS posltlve deflnlte If and only 1f 

P(£) does not van1sh ldentlca11y on any of the sublnterva1s 
44 

The fo11owlng theorem 18 modeled after theorem (3.2). 

Theorem (4.2): 1f (VF)(~(O» lS pos1t1ve def1n1te, then Newton's 

method - equat10n (4.10) - 18 always wel1-def1ned. 

Note that 1f x(O) 1S glven by x (0) = slgnum (d ) for €ach 1, - 1 1 

then P(~(O» 1S the character1st1C funct10n for the 1nterva1 [a,b] and 

~ x (l)N lS the second derlvat1ve of the natural CUb1C sp11ne 

J=l J J 

1nterpo1ant. 

If dk = 0 for some k, then we already know that any shape-preserv1ng 1nterpo1ant must be Ilnear on [tk ,tk

+2 ]. In fact any g E B must 

sat1sfy 

The Solutlon 10 thlS case lS of the form 

g =h{X +X } * [a,tk ] [t b] k+2, 

where 

m m m 

h (r aN) Xn - ( raN ) Xn + ( raN )Xn J=l J J + 1 J=l J J - 2 J=l J J 3 
45 

Slnce the value of a k cS ImmaterIal - the k-lh equatIon Fk(~) = dk of 

(4.11) beIng automatIcally satIsfIed - the number of equatIons and 

unknowns reduce by one each. For computatIonal convenIence we can 

stIll use (4.10) by settIng (VF)kk = 1, (VF)k,k+1 = 0, and (VF)k,k_1=0. 

Once we solve F(~) = £ we proceed to Integrate g* WhIch IS 

pIecewIse lInear (but not necessarIly contInuous, even If dk IS nonzero for each k) to obtaIn f* Wh1ch IS p1ecewlse CUblC. On the 

equatIon 

p (2)(t) 

1 

13 III + ~s //:,t )(t-t ) 1 (4.12) 

for t ~ t ~ t 1 If [t ,t 1] C ~3' 1 1+ 1 1+ 

p (2)(t) = (13 + (/:,13 //:,t )(t-t » 1 1 1 1 1 + (4.13) 

-(13 +(/:,13 //:,t )(t-t » 1 1 1 1- (4.14) 

for t ~ t ~ t 1 If [t ,t 1] C ~2 wIth boundary condl tlons 1 1+ 1 1+ 

and Yl+l' 
The functlon p lS elther a CUblC po1ynomla1 or pleceWlse CUblC 

1 

46 

glven hy two po1ynomla1s q~] and q12 deflned on separate sublnterva1s 

of [t ,t 1]' The Solutlon p to (4.11) lS glven by (l.18). The 

1 1+ 1 

Solutlon to (4.12) lS, dependlng on slgnum (6 ) and slgnum (6 1)' 1 1+ 

glven by (1.18), (3.16), (3.17), and (3.18). The Solutlon to (4.13) 

lS determlned by (l.18) If (3 ~ 0 and (3 1 ~ 0, by (3.16) If (3 > 0 

1 1+ 1 

and (3 1 < 0, by (3.17) If (3 < 0 and (3 1 > 0, and by (3.18) If 

1+ 1 1+ 

(3 ~ 0 and (3 1 ~ O. 

1 1+ 

F1gures (4.1), (4.3), (4.5} and (4.7) dlsp1ay the natural CUb1C 

sp11ne 1nterpo1ants to the glven data. F1gures (4.2), (4.4), (4.6), 

and (4.8) d1sp1ay the correspond1ng shape-preserv1ng 1nterpo1ants. 

Tables (4.1), (4.2), (4.3), and (4.4) glve convergence results for 

Newton's method. Note the quadrat1c convergence character1stlc of 

Newton's method. 

Appendlx B llsts a FORTRAN program for computlng the shape-preserv1ng CUb1C spl1ne 1nterpolant. 
· 1..Q15 

\ 

Flgure (4.1): The Natural CUblC Spllne Interpolant. 

Flgure (4.2): The Shape-Pre~ervlnci CUblC Spllne Interpolant. 

47 
. f\ /1 j \ 

,·1 \ J 

\ / 

\JI 

.I"""II,I"II'II"! • l'5..a l..flJlS 

Flgure:(4.3): The Natural CUblC Spllne Interpolant • 

. I-AlII5 

Figure (4.4): The Shape-Preservlng Cubic Spllne Interpolant. 

48 
49 

"'[" 

\ 

Flgure (4.5): The Natural CUblC Spllne Interpolant. 

U '!II,I"""",I,.! II.i,,!!!!!p! .. ,,'w.l''''''II'\ D • • .. '!II • I'- J~IS 

Flgure (4.6): The Sha~e-Prcscrvlng CUblC S~llnc Interpolant. 
Flgure (4.7): The Natural Cubic Spline Interpolant. 

ob""" •• I" •••• ",I •• "".,.I""""'/II,."".I •. ", ••• ,/ •• ",., •• 1, •• , ••• ,,/ o " • • • a ,. • ,. _-4U:IS 

Flgure (4.8): The Shape-Preservlng CUblC Spllne Interpolant. 

50 
Table 4.1 

Iterat10n Number 

1 

2 

3 

4 

5 

6 

7 

8 

II F (~ ( n )) - £11 2 

0.13 x 101 

0.67 x 100 

0.25 x 100 

0.42 x 10-1 

0.12 x 10-2 

0.88 x 10-6 

0.58 x 10-12 

0.64 x 10-13 

':>1 
Table 4.2 

Iterallon Number 

1 

2 

3 

4 

5 

6 

7 

8 

" F (~ ( n » - ~ II 2 

0.12 x 101 

0.56 x 100 

0121 x 100 

0.36 x 10-1 

0.11 x 10-2 

0.85 x 10-6 

0.54 x 10-12 

0.70 x 10-13 

52 
53 

Table 4.3 

Iteratlon Number /I F (~ ( n » - i II 2 

1 0.24 x 101 

2 0.16 x 101 

3 0.12 x 101 

4 0.90 x 100 

5 0.53 x 100 

6 0.20 x 100 

7 0.26 x 10-1 

8 0.42 x 10-3 

9 0.97 x 10-7 

10 0.37 x 10-12 

11 0.21 x 10-12 
Table 4.4 

Iteratlon Number 

1 

2 

3 

4 

5 

6 

7 

8 

54 

IIF(~(n) - £ 112 

0.29 x 101 

0.13 x 101 

0.50 x 100 

0.11 x 100 

0.56 x 10-2 

0.16 x 10-4 

0.13 x 10-9 

0.26 x 10-l2 
). COII..,trdlIH'd MIIlIIIII/dtIUI1 In d DUell Sp.)(C 

LeL C be d convex cone] n d normed dual space' X w~th predual Y. 

n Assume Yl'Y2""'Yn are elements of Y and deflne T: X~R by 

Tx 

Let B: = {x E C Tx = £} for a glven vector d. Conslder the 

problem 

Flnd x~~ E B such that II x-::_II ~ II x II for all x E B (5.1) 

of WhlCh (1.10), (3.2), and (4.5) are speclal cases. In thlS chapter 

we study eXlstence and characterlzatlon of Solutlons to (5.1). The 

follo~lng lemma glves sufflclent conctltlons for eX1stence of a Solut1on. 

Lemma(5.l): Ii B * lS nonempty, 1f C 1S weak closed, and 1f Y lS 

separable, then there eX1sts a solutlon to problem (5.1). 

Proof: Let y: = lnf {II x II 

sequence 1n C such that 

and 

x E C and Tx = d } • 

Tx = d 

n 

5') 

Let { x} be a 

n 

(5.2) 

(5.3) 
56 

for each n. SInce Y IS separable, by Alaoglu's theorem there eXIsls 

a weak* convergent 5ubsequence of {x } wIth weak* lImIt x. 

n 

IS weak'~ closed we have x E C, from (5.2) we have Tx = i, and from 

(5.3) we have Ilxll ~ y (and hence IIxll = y). ThIs completes the 

proof of the lemma. 

Throughout thIs chapter we assume that B IS nonempty, C IS weakoc 

closed, and Y IS separable. SInce x,~ = e If i = e, we assume also 

that d ~ e. The followIng proposItIon gIves us suffIcIent condItIons 

for C beIng weak* closed. 

ProposItIon (5.2): If C IS normed closed and If Y IS a reflexIve 

space, then C IS weak· closed. 

Proof: Assume {x} IS a sequence In C WIth weak* lImIt x. We want 

n 

to show that x IS In C. We do thIs by contradIctIon. If x IS not an 

element of C, then there eXIsts an element y (an element of both the 

dual and predualof X) whIch serves to separate x from C In the sense 

that 

for each nand 

x(y» K 

n 

X(y) < K 

for some constant K. ThIs ImplIes that 

11m x (y) ~ x(y) n 

n -+ co 
57 

Wh1Ch 1~ a contrad1ct1on. Therefore x E C and C 1S weak* closed. 

Th1s completes the proof of the propos1t1on. 

For y > 0 we deflne the convex set G(y) C Rn by 

G(y): {Tx X E C and II x II ~ y}. 

We now show that G(y) YG(l) and G(y) 1S closed. 

Propos1t1on (5.3): For each y > 0 we have G(y) = YG(l). 

Proof: By def1n1t1on 

G(y) {Tx x E C and II x II ~ y} 

{ Tx ~ E C and II x/y II ~ l} 

{T(x/y) ~ EC and Ilx/yll ~ I} 

y{ Tw w E C and II w II ~ l} 

yG(l) . 

Propos1t1on (5.4): The set G(l) 1S closed. 

Proof: Assume {z} 1S a sequence 1n G(l) wh1ch converges to~. We 

-n 

want to show that z 1S an element of G(l). Equ1valently, we want to 

show that x E C eXlsts such that II x II ~ I and Tx = z. 
For each n there eXlsts x E C such that "x II ~ 1 and Tx 

n n n 

By Alaoglu' s theorem there eXlsts a subsequence of II x II WhlCh 

n 

converges weaki~ to some x E C. Hence II x II ~ 1 and Tx = z. ThlS 

completes the proof of the proposltlon. 

We define 

i~ 

Y : = lnf{ y d E G(y)} • 

Equlvalently, 

y lnf{y There eXlsts x E C such thdt 

Tx = d and Ilx II ~ y} 

= lnf { II x II x E C and Tx = d} 

By lemma (5.1) we know that there eXlsts x_::_ E C such that 

58 

'I • 

11 

(5.4) 

(5.5) 

"i~4 II x* II = y and TXi~ = d. We call x~. an lnterpolant of mlOlma] norm. 

We now attempt to characterlze xh Vla the Hahn-Banach theorem. 

We begln by deflnlng a functlonal p : Y -+ R by 

p(y) sup {x(y) x E C and II x /I ~ 1 } • 

Notlce that lf C = X (the unconstralned problem), then p lS the norm 

on Y. In general, Slnce we are taklng the supremum over a subset of 

the closed UOlt ball U ln X, we have p(y) ~ "y" for all y E Y. 

Slnce e lS an element of C, we have p ~ O. In convex analysls p 15 

called the support functl0nal of the convex set {x E C: II x" ~ l}. 
59 

Slnce C 1S weak* closed, the supremum 1S atta1ned at some element of 

{x E C : Ilx II ::;; l}; that 1S, for any Y E Y lhere eXlsts an x (a 

funct10n of y) such that x E C, Ilx /I ::;; 1, and p(y) = x(y). In fact 

we have /Ix II = 1 un1ess x = O. The follow1ng two propos1t1ons reveal 

that p 1S cont1nuous, subadd1tlve, and pos1t1ve homogeneous. 

Lemma(S.S): The funct10nal p 1S cont1nuous. 

Proof: Assume Yl and Y2 are elements of Y and def1ne Y = Yl - Y2' 

Let x be the element 1n {x E C : /Ix /I ::;; I} such that p(Y2) = x(Y2)' 

Slnce Ix(Y)1 ::;; I/y", we have 

or 

Therefore, 

The elements Yl and Y2 can be 1nterchanged to obta1n 
60 

and hence 

Lemma (~: The funct10nal p 1S subadd1t1ve and pos1t1ve homogeneous 

(hence convex). 

Proof: Assume Yl and Y2 are 1n Y. To show that p 1S subadd1t1ve we 

must show that 

By deflm t10n 

x £ C and "x" ~ l} 

x £ C and /I x /I ;;; I} 

x E C and "x,,;;; I} 

Now assume a > 0 and Y E Y. To show that p 1S pos1t1ve homogeneous we must show that 

p(ay) ap(y). 
61 

By deflnl tl0n 

p(o.y) sup {x(o.y) x E C and II x II ;;l l} 

0. • sup {x(y) x E C and II x II ;;l l} 

o.p(y). 

ThlS completes the proof of the lemma. 

As an example we compute p for the case C = {x E L [a,b]: x ~ O} p 

where 1 < P < 00. For an arbltrary element g 1n L [a,b], the predual q 

of L [a,b] where p + q = pq, we have for any fEe w1th Ilf II ;;l 1 by p p 

the Mlnkowskl lnequallty 

r f(t)g(t)dt s r f(t)g/t)dt 
Assum1ng g+ ~ 0, let 

Then we have f £ C, II f II :::: 1, and 

P 

b I f(t)g(t)dt :::: Ilg)1 

a 

Hence 

Peg) sup (r f(t)g(t)dt f £ C and II f II 

a 

IIg+ II q' 

If g+ 0, then peg) :::: o. 

Lemma(5.7): For all a £ Rn we have 

n " n 

E a d 

1::::1 1 1 ~ y"p( Lay ). 

1::::1 1 1 

~ 1} 

P 

Proof: Sl.nce /:::: 1nf{y: £ £G(y)}, we have ££G(/~ + £) for any 

£ > O. Hence for every pos1t1ve 1nteger n there eX1sts x £ C 

m 

6'2 

(5.6) 
63 

~} 

such that TXm = d and "x

m 

1/ ;:; y + 11m. Therefore, for any ~ E R

n 

n n 

L: a d = 

1=1 1 1 

L: a x (y ) 

1=1 1 m 1 

n 

x ( L: a y ) 

fTl 1=1 1 1 

n 

;:i Ilxm II p (1:la1Y ) 

n 

;:i (y* + l/m)p( L: a y ). 

1=1 1 1 

Now let m -+00 to obta1n (5.6). Th1s completes the proof of the lemma. 

S1nce we know that G(Y>~) 1S closed from propos1t1on (5.4), we 

could have used x* 1n place of x 1n the proof of lemma (5.7). The 

m 

next lemma states that there PX1StS a nonzero vector BERn such that 

equall.ty holds 1n (5.6) . 

Propos1t10n (5.8): There CXlStS a vector ~ E R

n such that "B" 1 

and 

n 

B • d = Y ( L: By). 

1=1 1 1 

(5.7) 

Proof: The vector d 1S an element of G(y'<), but not 8n element of 

G(y'~ E) for any E > O. Hence the closed convex set G(y - E) and the 
64 

vector.Q. can be strIct ly sepdratcd by a hyperplune. ThIS Imp1 u~s 

the eXIstence of a nonzero vector ~(E) such that 

~ (c) • y < S(c) '.Q. 

for all y c G(y - c) and wlthout loss of generalIty we may assume 

that II S(c)" = 1. EqUlva1cnt1y, we have 

~ (e) • Tx <- ~ (c) • d 

and by the 11ncarlty of T 

n 

x( ~ S (e)y ) < §(c) • d 

1=1 1 1 

* for all x c C such that "x" ;:; y - c. Hence we obtaln 

n 1'-

(y - c)p( [ R (c)y ) < ~(c) • d 

1=] 1 1 

We can take the 11m t as F -+ 0 to 0 btaJ n a vec tor .§. such tha t "~,, 1 

and 

" n 

y "'oC ~ By) ~ ~ • d. 

1=1 1 1 
We have the reverse 1nequalLty from lemma (5.7) and therefore 

f? • {~ n 

d = Y pC r By). 1=1 1 1 

ThLs completes the proof of the lemma. 

Let ~ be a llnear funct10nal def1ned on the subspace 

by 

n 

A( r ex y ) 

1=1 1 1 

so that (5.6) can now be wr1tten 

n 

L: ex d 

1=1 1 1 

for all y c: S. 

65 

The Hahn-Banach theorem scates that there eXlsts an element w 

1n X such that 

w(y) A( y) for all YES (5.8) 

and 

for all y E Y. (5.9) 
61J 

Theorem (5.9): The Hahn-Banach extens10n w lS an 1nterpolant of 

m1n1mal norm. 

Proof: From (5.8) we see that Tw = d so that w 1nterpolates the data. 

To complete the proof we show that w E C and Ilwll = y*. 

We show that w 1S 1n C by contrad1ct10n. Assume w 1S not an 

element of C. Slnce C 1S weak* closed, there eX1sts an element Yo 

1n Y Wh1Ch str1ct1y separates w from C 1n the sense that 

w(y ) > x(y ) o 0 

for all x E C. (5.10) 

Slnce C lS a cone we have AX E C whenever A > 0 and X E C. Hence 

(5.10) 1mphes 

(or p(y ) a 0) and 

o ~ x(y ) o 

for all x E C 

w(y ) > O. 

o 

However, from (5.9) and (5.12) we have 

wh1ch 1S a contrad1ct10n. Hence w must be an element of C. 

(5.11) 

(5.12) 
67 

Lastly, we show that II w II y. We already know that 

;; 

y ~ II w /I (5.13) 

slnce wEB (w E C and Tw = ~). Because p lS bounded above by the 

norm on Y, (5.9) Y1e1ds 

-~~ 

w(y)~y I/yl/ for all y E Y 

and hence 

I/w II ~ y * • (5.14) 

Taken together, (5.13) and (5.14) 1mply that I/wl/ y. * ThlS completes 

the proof of the theorem. 

Recall that for a glven element y 1n Y there eX1sts an element 

o 

elther II Xo II 

x(y). Furthermore, o 

1 or x lS the zero element. The follow1ng lemma 

o 

w1ll lead us to the conclus1on that, 1f p lS d1fferent1able at y , 

o 

~emma (~: Let f be a funct10nal def1ned on a normed llnear space 

Z. If f lS d1fferent1able at x E Z and 1f there eX1sts a llnear 

o 

funct1ona1 A SUcll that 
68 

f (L. ) + A( ?-Z ) :$ f (z) a a (5.15) 

for all z 1n some ne1ghborhood of zo' then A (Vf)(z ). 

o 

Proof: Let z 

pelds 

z ~ tu where t > 0 and U E Z. Inequal1ty (5.lS) o 

fez + tu) - fez ) 

A(U) :$ o 0 

t (5.16) 

Slnce (5.16) holds for all t > 0 (and suff1c1ently small) and for all 

u E Z, we have A~(Vf)(z). Subst1tut1ng -u for u 1n (5.16) Y1e1ds 

o 

f(7 - tu) - f(L. ) 

A (u) ~ o 0 

t (5.17) 

for all t > 0 (and suff1c1ent1y 5ma1l) and for all u E Z. Taken 

together, (5.16) and (5.17) 1mply A = (Vf)(z ). 

o 

Corollary (5.11): 

Proof: Slnce p(y ) o 

If P 1S d1fferent1able at y E Y, then p'(y ) o ---- 0 

p (y ) + x (y - y );;; p(y) o 0 0 

x • 

o 
for all y £: Y. By the prcvl.ous lemma we have p'(y ) o 

completes the proof of the corollary. 

InequalIty (5.6) motIvates the problem 

n 

x. Thl.s 

o 

69 

Inf{ P ( Lay ) : a • d = I} . a 1=1 1 1 

(5.18) 

NotIce that If a IS any vector satIsfyIng £ • d = 1 and If x IS any 

element of B, then 

and hence 

n 

1 = Lad 

1=1 1 1 

n 

= xC Lay ) 1=1 1 1 

n 

~ II x II pC Lay ) 

1=1 1 1 

n 

pC L Ct Y ) 

1=1 1 1 

_1_ 

~ II xII 

ThIs l~plles that the InfImum IS posItIve (and, In fact, IS bounded 

{~ -1 * n below by (y) • If the InfImum IS attaIned at some £ £: R and If p 

n {: 

IS dIfferentIable al L Ct Y , then we are led to a SolutIon to (5.1) as 

1=1 1 1 
70 

the next theorem reveals. 

~4 .:~ 

Theorem (5.12): If there eX1sts a E Rn such that a • d = 1 and 

n 

p( l: a*y ) 

1=1 1 1 

and 1f P 1S d1fferent1able at 

n 

1nf { p ( l: a y ) : a· d = l} 

a 1=1 1 1 

then 

* n 

y pl( l: a y ) 

1=1 1 1 

lS an 1nterpolant of m1n1mal norm. 

Proof: Problem (5.18) has Lagrang1an 

L(~, A) 

·k 

If there eX1sts a Solut10n a 

* * 

n 

A( L: a d - 1). 

1=1 1 J 

~~ 

to (5.18), then there eX1sts A 

that (~ ,A ) lS a stat10nary p01nt of (5.19). Hence 

o 1=1,2, ... ,n 

(5.19) 

so 

(5.20) 
n " 

where x = p'C ~ a:Y1)' 

1=1 

x £ C, II xii 1, and a d = 1. 

vJe fust show that A * > O. MUlt1p1y C 5.20) by a * and 1 

sum over 1 to obtaln 

I C n * ) Slnce x = p ~ a

1 Y1 ' we have 

1=1 

so that 

n ';'_ 

xC ~ a y ) 

1:1 1 1 

A • 

n r- pC z: Ct y ) 

1=1 1 1 

71 
72 

Actually, we know that Slnee the 1nf1mum lS pos1L1ve, we have A-:: > O. 

We Cdn also show th1s by contrad1ct10n. If A = 0, then 

n ~~ 

xC L: a y ) :;; 0 

1=1 1 1 

for all x E.. C. (5.21) 

Let S be any 1nterpo1ant 1n C. (We know that there eX1sts an 1nterpo1ant 1n C Slnce B lS nonempty.) Then 

n ~~ 

s( Lay ) 

1=1 1 1 

n ~r 

L: a d 

1=1 1 1 

Wh1Ch contrad1cts (5.21). Therefore, A * > O. 

1 

Now we show that 1. From (5.70) we see 

lnterpolant 1n C. Hence 

II x" 

or 

15 an 

Let w be an ]nterpo1ant of m1n1ma1 norm sat1sfY1ng (5.9). Then 

(5.22) 
Equ1valently, we have 

wh1ch leads to 

n ~~ 

w( L 0. Y ) ;;; 

1=1 1 1 

n 

w( L 0."- Y ) 5 

1 1 1=1 

~~ n ,. 

Y x( L 0. 'y ) 

1=1 1 J 

Taken together, (5.22) and (5.23) 1mply 

-'~ i~ 

1 Y A • 

ThJS concludes the proof of the theorem. 

73 

(5.23) 

We cons1der now the problem of determ1n1ng when the 1nf1mum 1S 

atta1ned 1n (5.18). From proposlt1on (5.8) we know that there eX1st 

a nonzero vector ~ such that 

{:_ n 

o ;;; ~ • i = Y p( L By). 

1-=1 1 1 
74 

~~ 

If B °d > 0, then the 1nf1mum lS atta1ned 1n (5.18) at a = ~/(~ 0 i). 

Propos1t10n (5.13): If d 1S 1n the relat1ve 1nter10r of 

s: { r r E G(y) for some y }, 

then there eX1sts a vector ~ such that 

>c n 

1 = ~ °d ~ y' p( E By). 

1=1 1 1 

Proof: We prove by contradlct10n. Assume that every vector ~ wh1ch 

,. n 

Ii . d = y' p( E By) - 1=1 1 1 

also satlsfles Ii 0 i = O. W1thout loss of general1ty 1t can be 

assumed that there eX1sts a nonzero vector ~ such that 

and 

,. n 

o = Ii 0 d = yAp( E B

1 Y1) 

1=1 

for all y E G (/~). 
75 

In any relat1ve ne1ghborhood of £ there 1S a vector z such that 

~ • ~ < O. If z were an element of S, then there would be an 

element ~ 1n G(y ) such that z = ar for some a > O. However, we 

would then have 

~ • z 

Wh1Ch lS a contrad1ct1on. Therefore Z lS not an element of Sand d 

lS not 1n the relatlve lnter10r of S. Th1s completes the proof of 

the propobltlon. 
References 

[de B( 1) ] C. de Boor (1978): "A Prac tlca1 GUlde to Sphnes," 

Sprlnger-Ver1ag, New York. 

[de B(2)] C. de Boor (1976): On "best" lnterpoJotlon, J. Approx. 

Theory, 16:28-42. 

[IMS] L.D. lrvlne, P.W. Smlth, S.P. Marln (1985): Con~tralned 

76 

enterpo1atlon and smoothlng, General Motors Research 

Report No. MA-312, Warren, MJchlgan. 

f L] D.C. Luenberger (1973): "Optlmlzatlon by Vector Space 

lHSSW] 

Methods," John We1ey ;:JOd Sons, Inc., New York. 

C.A. Mlcc~e111, P.W. Smlth, J. Swetlts, J.D. Ward (1985): 

Constralned L approxlmatlon, J. Constructlve 

p 

Approxlmatlon, 1:93-102. 
AppendIX A 

A Program for ConstructIng the Natural CUbIC SplIne Interpolant 

to GIven Data. 

77 
00001 

OOOO~C 

0OOO3C 

00004C 

OOOO~ 

00006 

00007 

OOOOBC 

0OO09C 

0OO10C 

oOOtlC 

0001:2C 

0OO13C 

OOO14C 

0OO15C 

000lbC 

00017C 

00018C 

00019C 

OOO::OC 

00021C 

0002~C 

000:)3 

000:24 

00025 

0002t,C 

00027C 

00028C 

000~9C 

00030 

0003t 

0003~ 

00033 100 

00034C 

00035C 

0003bC 

00037C 

00038C 

0OO39C 

()0040C 

00041 

00042 

00043 

00044 

00045 

0004b 

00047 

OOOL:8 200 

00049 

000::;0 

F'F:DGRAI1 UNCON (INPUl • OUTPUT, 1 AF'[5=HJf-'lIT ,l APE6=DllTPUT) 

WE FORM THE NATU~AL CUBIC SPLINE INTERPOLANT. 

INTE.GER N.M.I 

REAL T(50),F(50),D(50),X(SO),A(SO),F'P(4,50) 

REAL AA(SO),BB(SO),CC(SO) 

78 

lHE ARRAYS (T) ANn (F) - EACH OF SIZE M, THE NUMBER 

OF [lInf,. POINlS - C(lNTAIN lHE COMPONENTS OF TH[ IIA1A. 

THE DATA FILE IS OF THE FOLLOWING FORM 

11 

T(1),F(J) 

T<:2).F(~) 

T(M),F(M) 

WHERE W[ ASSUME (T) HAS STRICTLY INCREASING COMPONENTS. 

REA[I(3,*' M 

READ(3,l) (T(I),F(I), I=l,M) 

N= ~1-:2 

THE ARRAY (D) CONSISTS OF THE SCALED 

SECOND DIVIDED DIFFERENCES. 

DO 100 I=l,N 

D(I)= ( F(I+2)-F(I+l) )/( T(I+2)-T(I+l) ) 

C - ( F(I+l)-F(I) )/( T(I+l)-T(I) ) 

CONTtNUE 

THE SECOND DERIVATIVE OF THE NATURAL CUBIC SPLINE 

INTERPOLANT IS A LINEAR COMBINATION OF LINEAR B-SPLINES. 

WE CALCULATE THE COEFFICIENTS. 

M(l)= 0.0 

8B(I)= (T(3)-T(1»/3.0 

CC(I)= (T(~)-T(~»/6.0 

no 200 I=2,N-l 

AA(I)= (T(ttl)-T(t»/6.0 

B8(1)= (T(I+2)-T(I»/3.0 

CC(I)= (T(I+:2)-T(I+l»/6.0 

CONTINUE 

Mon= (T(Nt! )-Ton )/6.0 

8BeN)= (T(Nt2)-T(N»J3.0 
00051 

OOO~2 

00053C 

0nD~4C 

00055C 

OOO~6 

00057 

000~8 

00059 

00060 300 

D006lC 

OOD6~C 

00063C 

00064C 

00065C 

00066C 

OOD67C 

00068 

00069 

00070 

00071 

00072 

00073 

00074 

00075 

00076 400 

00077 

00078 

00079 

00080 

0008lC 

0008~C 

00083C 

00084 

00085 

00086 450 

00087 500 

00088e 

00089C 

00090C 

00091C 

00092C 

00093C 

00094C 

00095 

00096 

00097 

00098C 

00099 

00100 

CC(N)= 0.0 

[ALL TRJD(AA.BB,CC,D,N) 

A(l)= 0.0 

A(M)= 0.0 

DO 300 J=2,N11 

A(I)= D(I-I) 

CONTINUE 

NOW WE COMPIITE THE NUMBERS PP(J,I) - THE VALUE 

OF THE (J-l)ST DERIVATIVE OF THE NATURAL CUBIC 

SPLINE INTERPOLANT EVALUA1ED AT T(I). 

no 400 K=I,Ntl 

DF= F(Ktl)-F(h) 

DT= l(Ktl)-T(K) 

DA= A(Ktl)-A(K) 

PP(4,K)= DA/DT 

PP(3,K)= A(h) 

PP(~,K)= DF/DT - (A(K)/~. + DA/6.)*DT 

PP(I,K)= F(K) 

CON1JNUE 

PP(4,M)= 0.0 

PP(3,M)= 0.0 PP(~,M)= 0.0 

PP(l,M)= F(M) 

DO 500 h=l,M 

WRITEC6,450) K,T(K),(PP(I,K), 1=1,4) 

FORMAT(5X,I5,5FI4.6) 

CONT[NUE 

WE CREATE A DATA FILE FOR PLOTTING THE (JDER)-TH 

DERIVATIVE OF THE NATURAL CUBIC SPLINE INTERPOLANT 

DY EVALUATING IT AT (MM) EQUALLY SPACED POINTS, 

INCLUDING THE ENDPOINTS. WE ASSUME THAT (JDER) 

HAS VALUE 0, 1, 2, OR 3. 

JDER= 0 

MM= 201 

CALL DATAFL(l,PP,M,MM,JDER) 

STOP 

END 

79 
00001 

00OO2C 

00003C 

00004C 

00005C 

00006C 

00007 

00008 

OOOO? 

00010 

00011 

00012 

00013 

00014 

00015C 

00016C 

00017C 

00018 

00019 

00020 

00021 200 

000:'>2 300 

00023 

00024 

00025C 

OOO:'>bC 

00027C 

000~8C 

OOO~9 

00030 

00031 

00032 

00033 

00034 

00035 400 

00036 

00037 450 

00038 500 

00039 

00040 

Slll-lF;OUTI NE=" !1M T MFU TX, F'P, LI ,MM, .I!1EFn 

WE CREA1E A DATA FILE FOR PL01T1NG THE (JDER)-TH 

DER[VMTIVE OF THE PIECEWISE CUBIC POLYNOMIAL. WE 

ASSUME (JUER) HAS VALUE 0, 1, 2, OR 3. 

INTEGER LI,MM,JDER 

REAL TX(100),PP(4,100) 

LEFT= 1 

liMON!:> MM - 1 

WR 11 E ': 4 , * ) MM 

XE= ( TX(LI)-TX(l) )/FLOAT(MMONE) 

[10 500 IP=l,MN 

XT= TX(l) + XE*FLOAT(IP-l) 

WE FIND THE INTERVAL IN WHICH THE POINT (XT) LIES. 

IF ( LEFT .NE. LI ) THEN 

[10 200 IS~LEFT,LI-l 

IF ( XT .LT. TX(IStl) ) GO TO 300 

CONTINUE 

CONTINUE 

END IF 

LEFT= IS 

WE NOW COMPUTE THE VALUE OF THE POLYNOMIAL AT 

THE POINT (XT) BY USING MESTED MULTIPLICATION. 

H= XT - TX (LEFT> 

FAC= 4.0 - FLOAT(JDER) 

YT= 0.0 

!IO 400 M=4,JDERtl,-1 

YT= CYT/FAC)*H + PP(M,LEFT) 

FAC= FAC - 1.0 

CONTINUE 

WRITE(4,450) XT,YT 

FORMAT(F8.4,E18.9) 

CONTINUE 

RETURN 

END 

80 
00001 

00002 

OOOO:~ 

00004 

00005 

00006 

00007 

00008 

00009 

00010 

00011 

0001~ 111 

00013 

00014 

00015 

00016 2~2 

00017 

00018 

~ilIFeHOU1 JNE 1 F:JlI (SUIl ,DIAG, SUP, Fe, N) 

INTEGFf.: N, I 

F;EAl }-l( N) ,I11flG (N) , SlllH N) , SUP (N) 

IF (N.LE.t) THEN 

B(l)= B(I)/DIAG(l) 

RETUF;N 

ENII IF 

DO 111 I=:?,N 

SUB(I)= SUFe(I)/DIAG(I-l) 

DIAG(I)= DIAG(I) - SUB(I)*SUP(I-l) 

B(Il= B(I) - SUB(I)*B(I-l) 

CONTINUE 

B(N)= B(N)/DIAG(N) 

DO 2~2 I=N-l,l,-l 

BCI)= (B(I)-SUP(J)*B(Itl»/DIAG(I) 

CONTINUE 

RETUF;N 

END 

81 
~pendlX B 

A Program for Constructlng the Shape-Preservlng CUblC 

Spllne Interpolant to Glven Data 

82 
0000] 

00002C 

00003C 

00004C 

OOOO~C 

00006C 

00007C 

OOOOSC 

00009C 

OOOtOC 

00011C 

00012C 

00013C 

00014C 

00015C 

00016 

00017 

00018 

00019 

00020C 

00021C 

00022C 

00023C 

00024C 

00025C 

0002bC 

00027C 

0002SC 

00029C 

00030C 

('0031C 

COOJ2C 

00033C 

0003~C 

OO",:,~C 

C~~3~C 

NHiJ7 

l..,OOJf 

000:;9 

Ov040C 

00041C 

0(·042: 

O('043C 

OO(144i:' 

()004~,C 

O·:>Oi!t, 

\,0047 

C")04;iC 

COO~":'; 

00050C 

['r:OGF;AM MIN! ItWUT ,rJlITr'UT, fflF'E"5=-l/WUT, 1hPEb=OUH'lIl) 

WE. CDMf'Ul E A SHM'E -F f.:E.SERV J NG I NTE:.RF'OLANT 

TO GI ..... EN DMA. 

NOTE ON THE SIZE OF THE ARRAYS: 

THE ARRAYS (T), (F), AND (ti, MUST BE OF LENGTH 

AT LEAST H, THE NUMBER OF DATA POINTS. THE 

ARRAY (TX) AND THE SECOND COMPONENT OF THE 

ARRAY (rp) SHOULD BE or LENGTH 2M. THE ARRAYS 

(x), (n. AND (D) MUSf BE OF LENGTH AT LEAST M-2. 

THE ARRA) (ID) HUST BE OF LENGTH A1 LEAST M-l. 

REAL T(50),F(50),X(50),Y(SO),A(50) 

REAL TX(100),PP(4,100) ,TL, TF.,c.PS 

INTEGER M,N,ITMAX,l,J,IFLAG,MH 

COMMON D(50),rp(~o) 

83 

THE ARRAYS (T) ANn (F) - EACH OF SIZE M, THE NUMBER 

OF DATA POINTS - CONTAIN THE COMPONENTS OF THE DATA. 

HIE DATA FILE IS OF THE FOLLOWING FORM 

M 

T(U,F(l) 

T(2),F(2) 

; 0'1, ,F(M) 

wHO;i ·,Jf ASSU'1E (1) HAS S~R1CTL,( INCREASING COMPONE.tns. 

' • .=..Art, -.~} M 

~.Et'\tl(:,.*) ~i(!),f(I>, "=:L;l) 

N~ M-2 

, ; f'~" r r (I ':-- (, ._ t- f1fJ ITT \!E NUMBER USED TO TEST FOR 

COr~\Ln:'["'U H~ HEWTOWS METHOD - SUBROUTINE (ZERO). 

(ITMAX' IS THE MAXIMUM NUMBER OF I1ERATIONS 

WHICH wr PERK[T ~OR NEWTONS METHOD TO CONVfRGE. 

THL (Ihl-, '.l ~ X) IS 1 HE I,NDT SC::QUENCE (T) WITH THE 

ElWr:'O;:H3 lL r,j~D TR [tE'_eTED. 
OOO~IC 

0005:2 

()OO~,3 

00054 

C)005~, 

00056 120 

00057C 

00058C 

00059C 

(l0060C 

00061C 

00062C 

00063C 

00064C 

00065C 

00066 

00067 

00068 

00069 

00070 

00071 

0007:2 130 

00073 135 

00074 

00075C 

00076 

00077 

00078 

00079 

00080 140 

oooalC 

00082C 

00083C 

00084C 

00085C 

00086C 

()OO87C 

00088 

00089C 

00090C 

ooor,> 1 

0009:2 

000S)3 

00094 

0009:=i 

00096C 

00097C 

00098 145 

00099C 

OOlOOC 

0010t 

TL= T(1) 

TR= T':M) 

DO 120 I=l,N 

XCI)= T(HU 

CONTINUE 

THE ARRAY CD) CONSISTS OF THE SCALED 

SECOND DIVIDED DIFFERENCES. 

84 

IT IS IMPORTANT THAT WE IDEN1IFY DIVIDED DIFFERENCES 

WHICH ARE ZERO. THIS MEANS THAT WE MUST COMPARE TWO 

FLOATING-POINT NUMBERS. 10 DO THIS WE ASSUME DCK) IS 

ZERO IF DCK) IS SMALL. 

XEPS= 1.0 

DO 130 J=1,20 

XE.PS= XEPS/I0. 

Z= 1.0 + XEPS 

IF C Z .EQ. 1.0 ) GO TO 135 

YEPS= XEPS 

CONTINUE 

CONTINUE 

YEPS= YEPS*1000. 

DO 140 K=l,N 

DCK)= C FCKt2)-F(Ktl) )/( T(Kt2)-T(Ktl) ) 

C - C FCKtl)-FCK) )/C T(~+l)-TCK) ) 

IF C ABSCDCK» .LE. YEPS ) D(K)= 0.0 

CON1INUE 

THE INITIAL GUESS (Y) FOR NEWTON'S METHOD 

WILL YIELD THE SECOND DERIVATIVE OF THE 

NATURAL SPLINE SOLUTION, EXCEPT POSSIBLY 

WHEN DCK)= 0.0 FOR SOME K. 

DO 145 ~=l,N 

IF C DCK) .GT. 0.0 ) THEN 

ELSE 

YCI\)= -1.0 

END IF 

CONTINUE 

WRITE (6,150) 
00102 150 

00103 

00104 160 

OOlO~J 

00106 170 

00107C 

0010eC 

OOl09C 

00110C 

OOllle 

00112C 

f)0113C 

00114C 

00115C 

OO1l6C 

00117C 

OOllSC 

OOl19C 

00120C 

001:>lC 

00122C 

001:',) 

00124 

0OI:.!~J 

00126 

00127 180 

00128 

00129 

00130 

00131 

0013? 

00133C 

00134 

00135 

001~t. 

00:13:' 

oU38 

00139C 

00140C 

00141C 

0014:2C 

00143C 

00144 

00145 

00146 185 

()0147C 

OO14GC 

(10149C 

0·:l150C 

00151C 

()O152 

FORMATe/,' DATA VALUES ',/) 

WFiIlE(6,160) (lIe 1), l=-l,N) 

FORMAT(5X,4E12.6) 

W~; 11 E ( 6 r 170 ) 

FOF:MA T (f f) 

ID(K)= 1 INDICATES THAT THE INTERPOLATING FUNC1ION 

IS CONSTRAINED TO BE CONVEX ON [TCK),T(Ktl)] 

AND, HENCE, ITS SECOND DERIVATIVE IS CONSTRAINED 

TO BE NONNEGATIVE ON THIS INTERVAL. 

ID(K)= -1 INDICATES THAT THE INTERPOLATING FUNCTION 

IS CONSTRAINED TO BE CONCAVE ON [T(K),T(Ktl)] 

AND, HENCE, 11S SECOND DERIVATIVE IS CONSTRAINED 

TO BE NONPOSITIVE ON THIS INTERVAL. 

IDCK)= 0 INDICATES lHAT THE INTERPOLATING FUNCTION 

IS UNCONSTRAINED ON [T(K),T(Ktl)J. 

DO 180 I=l,N-l 

ID(It1)= 0 

IF (D(I).GE.O.O .AND. DCIil).GE.O.O) ID(Itl)= 1 

IF (D(I).LE.O.O .AND. D(Itl).LE.O.O) ID(I+1)= -1 

C(lNTINUE 

IF ( D(1) .GE. 0.0 ) THEN 

IfI(1)= 1 

ELSE 

III( 1)= -1 

END IF 

IF ( [I(N) .GE. 0.0 ) 1HEN 

Ht(NtU= 1 

ELSE 

1[l(Ntl).:: -1 

ENit IF 

IF A NONZERO nATA VALllf D( I) LIES BETWEEN TWO 

ZERO DATA VALUES D(I-l) AND D(Itl), THEN D(I) 

IS lAKEN TO BE ZERO FOR COMPUTATIONAL PURPOSES. 

DO 185 I=::!,N-l 

IF ( D(I-1).EO.0.0 .AND. D(ltl).EQ.O.O ) D(I)= 0.0 

CONTINUE 

SUBROUTINE (ZERO) CALCULATES lhE PIECEWISE 

LINEAR SECOND DERIVATIVE OF THE SHAPE- PRES~RVING INTERPOLANT. 

CALL ZERO(Y,X,N,ITMAX,EPS,IFLAG,lL,TR) 

85 
00153C 

00154 

001~5 

0015t, 

00157 

()01~18 190 

00159C 

00160 

0016J ~OO 

00162 

00163 :210 

00164 

00165 :?:>O 

00166C 

0Olf>7C 

O(,168C 

00169C 

00170 

00l71e 

00172 

00173 230 

00174 

00175 

00176 240 

00177 250 

00178 

00179 260 

00180 

00181 ::>70 

00182C 

00183C 

00184C 

001El5C 

00186C 

00187C 

00188C 

00189C 

00190 

001 S'1 

00192 

00193C 

00194 

00195 

AO)= 0.0 

AI,M)= 0.0 

[1(1 190 I=2,Nt1 

A(l)= Y( [-1) 

CON1JNUE 

WF:ITE (6,200) 

FORMAT(/,' PIECEWISE LINEAR :>ND DERIVATIVE ',I) 

WRITE(6,210) (l,T(I),A(I), I=1,M) 

FORMAT(5X,I5,' ( ',F14.6,' , ·,F14.6,' )') 

WRITE (6,220) 

FORMAT (J /) 

SUBROUTINE (POLY) INTEGRATES THE RESULT 

FROM SUBROUTINE (ZERO). 

CALL POLY(A,T,PP,M,F,LI,TX) 

WRIlE(6,230) 

86 

FORMAT(/,' KNOTS AND COEFFICIENTS OF PIECEWISE CUBIC',!) 

DO 250 I=l,LI 

WRITE(6,240) I,TX(I),(PP(J,I), J=1,4) 

FORMAT(5X,I5,5F14.6) 

CONTINUE 

WRITE(6,260) IFLAG 

FORMAT(/,' ERROR CODE = ',15,/) 

WRITE(6,270) ITMnX 

FORMAT(/,' NUMBER OF ITERATIONS =',15,) 

SUBROUTINE (DATAFL) IS US~D TO CREA1E A 

DATA FILE FOR PLOTTING. WE EVALUATE THE 

UDER)-TH DERIVATIVE OF THE PIECEWISE CUBIC 

POLYNOMIAL AT MM EQUALLY SPACED POINTS, 

INCLUDING THE ENDPOINTS TL AND lR. WE 

ASSUME (JDER) HAS VALUE 0, 1, 2, OR 3. 

MM= 201 

JDER-= C 

CALL UAfAFL(TX,PP,LI,MM,JDER) 

STOP 

ENfl 
0000] 

00002C 

00003C 

(10004 

00005 

OOOOt. 

00007 

OOOOS 

00OO9C 

0OO10C 

OOOlle 

(lOO12C 

OClO13C 

0OO14C 

00Ol5C 

0OO16C 

0OO17C 

OOOlSC 

0OO19C 

0OO20C 

00021C 

0OO22C 

00023C 

OOO~4C 

0OO25C 

OOO;'6C 

OOO:>7C 

0OO23C 

VOO29C 

')OO30C 

0OO31C 

{IOO~2C 

00033( 

ooos>;c 

))03~,c.. 

... J Ij J:~6 

')00:::7 100 

vOO38 C 

00039 

(,OO40C 

0OO41C 

00047: 

(JOO43C 

0OO44C 

1;(IO<i:': 

('OO4tC 

·::004 i C 

0004Sr; 

'cI OO49C 

OOo~,(lr 

SUBROUTINE ZERO(A,X,N,JTMAX,EPS,IFLAG,TL,TR) 

INTEGrR N,ITMAX,~,J,LJ,L,IFLAG 

REAL A(N),X(N),FXI50),AL,XL,AR,XR,DT,DA,T,W 

REAL SUB(50l,DIAGI50),SUP(50),H(50),SUM1,SUM2 

REAL RATIO,GLEFT,GRIGH,EPS,FNORM1,TL,TR 

COMMON D(50),ID(50) 

INPUT PARAMETERS: 

A ••• INITIAL ESTIMATE FOR NEWlON'S METHOD. 

X ••• KNOT SEQUENCE WITH lHE ENDPOINTS DELETED. 

~7 

N ••• THE SIZE OF lHE ARRAY (A): l~E NUMBER OF UNKNOWNS. 

I1MAX ••• MAXIMUM NUMBER OF ITERATIONS FOR NEWTON'S METHOD. 

EPS ••• PARAMETER USED TO TEST FOR CONVERGENCE. 

TL,TR ••• LEF1- AND RIGHT-ENDPOINTS OF THE 

INTERVAL RESPECTIVELY. 

OUTPUT PARAMETERS! 

A ••• ThE CALCULATED ZERO IF CONVERGENCE OCCURRED. 

ITMAX ••• NUMBER OF ITERATIONS RE~UIRED FOR NEWTON'S 

METHOD TO CONVERGE. 

IFLAG ••• IFLAG= 1: CONVERGENCE INDICATED BY COMPARING 

THE Ll NORMS OF THE ITERATES 

IFLAG= 2: NUMBER OF ITERATIONS EXCEEDED ITMAX. 

F'F.IIJT 100 

FOr'MAl" (' ITERATION NUnBH: AND RESIDUAL! I,! 

I QUADRATIC LCNVERGENCE IS EXPECTED.',!) 

DO 350 LJ=l t ITMAX 

THE AR~AYS (SUB), (DIAG). AND (SUPl CONTAIN 

TH[ ELEMENTS OF THE TRIDIAGONAL POSITIVE-DEFINITE 

JACOBIAN MATRIX (J), EVALUATED AT THE VECTOR (A). 

IT SHOULD BE NOTED. THAT THE MATRIX EQUATlON SOLVER, 

THE SUP~OUTINE (TRln), DOES NOT TAKE ADVANTAGE OF 

THE SrhMETRY OF (J). HENCE (SUB) AND (SUP) ARE 

BOT~ NECESSARY. ALTHOUGH SUB(\)=SUP(K-]), EQUATIONS 

FO~ qO!H ARRAYS ARE WRIfTEN OUT IN FULL. 

l~ D(~)=O.O FOR SOME ~, lHEN THE NUMBER 
00051[' 

00052C 

OOO~J3C 

()0054C 

OOO~'~JC 

OOO~,6 

00057C 

000~J8 

00059 

00060 

00061 

00062 

00063 

000t..4 

00065C 

00066C 

00067 

00068 

00069 

00070 

OOO?] 

00072 

00073 

00074C 

O()075C 

00076 

00077 

0001'8 

00079 

OOOElOC 

00081 

00082 

00083 

00084 

00085C 

00086 

00087 

00088C 

00089 

00090C 

00091 

ooonc 

00093 

00094C 

00095 

00096 

0009?C 

00098 

00099C 

00100 

00101 

OF UN~NOWNS (AND EQUATIONS) REDUCE. IN ORDER 

TO PERMIT lHE COMPUTATION OF ONE JAfOBIAN 

MATRIX THE PROGF:AM SETS SUB(~\)=SlIP(I'\-l)=O.O 

AND DIAG(K)=1.0. 

DO 125 I\=1,N 

IF (K.EQ.1) lHEN 

AL= 0.0 

XL= lL 

ELSE 

Al= AO,-l) 

XL= X (\\-1 ) 

END IF 

IF (~.EO.N) THEN 

AR= 0.0 

XR= TR 

ELSE" 

AR= A(Kt1) 

XR= X (K+1 ) 

END IF 

IF AL.GE.O.O .AND. A(K).GE.O.O ) J1= 1 

IF AL.LT.O.O .AND. A(I\).GE.O.O ) Jl= ~ 

IF AL.GE.O.O .AND. A(K).LT.O.O ) Jl= 3 

IF AL.LE.O.O .AND. A(").LE.O.O Jl= 4 

IF ( A(I\).GE.O.O .AND. AFi.GE.O.O 

IF ( AO~) .LT .0.0 • Af~D. AR.GE.O .0 

IF ( A(K).GE.O.O .AND. AR.LT.O.O 

IF A ( ~: ) • LE .0 .0 .AND. AR.LE.O.O 

DT= X(K)-XL 

DA= A(I\)-AL 

IF ( I D it,) .EO. n THEN 

IF 0\ .NE.!) THEN 

IF (.11. EO. 1> THEN 

SUB("')= DT/6.0 

GLEFf= DT/3.0 

ELSE IF (Jl.EQ.2) THEN 

1= XL-CDT/DAJ*AL 

W= O.5~( X(KltT ) 

J2= 1 

1'1= '1 

J2= 3 

J2= 4 

88 
()010:> 

00103 

001()4 

001 O~~ 

(1010bC 

O()107 

OOlOaC 

00109 

00110 

OOll] 

(10112 

0011"3 

00114 

OOll ,)C 

00116 

00117e 

(j01l8 

00119 

00120C 

00121 

00122C 

001::3 

00124C 

001~5 

00126 

00:i:':7 

00128(' 

00129 

00l30( 

00131 

00132C 

,')133 

0\' 11.a 

!lO::: C 

001,FC 

f, "<i 3S 

G0139C 

00140 

00141C 

(.'C14~ 

00143 

O(/t44C 

00145 

00146C 

00147 

00148 

001" Y 

00 1 50 

001 ~.J 

00152 

SUB(K)= (x(t\)-1) io.O * ( «1-XI )/Dn*«x(tn-T) IrlT) 

c + 4.0*«W-XL)/DT)*«X(~)-W)/DT) ) 

GlEFT= (X(~)-T)/6.0 * ( (CT-XL)/DT)**2 

C t 4.0*«(W-~L)/DT)**2) + 1.0 ) 

ELSE IF (J1.EQ.3) THEN 

T::. XL-(DT/Dt'i)*M_ 

1rJ= 0.5* ( T+ XL ) 

SUB(K)= (T-XL)/6.0 * ( 4.0*«W-XL)/DT)*«XCK)-W)/DT) 

c + \ (T-XU/DT)*«X(~)-T)/[ln ) 

GLEFT= (1-XL)/6.0 * ( 4.0*«(W-XL)/DT)**2) 

c + «T-XL)/[lr)~*2 ) 

ELSE IF (Jl.En.4) THEN 

SUB(~:):: 0.0 

GLEFT= 0.0 

END IF 

ELSf IF (K.EQ.l) THEN 

SUB(1) = 0.0 

GLEFT= 0.0 

IF (Jt.EQ.l) GLEFT= DT/3.0 

nw IF 

ELSE !F ( ID(K) .ED. 0 THEN 

SUP(K)::. DT/6.0 

(?LEFT= £IT /3.0 

F'_SE ! F ( III ( K ) • EQ. -1 ) THEN 

IF (I\.NE.1) THEN 

IF (J1.EQ.4) :HEN 

5UB(I\)= lIT/6.;) 

GLEFT::. rll /3. \..) 

ELSE IF (J1.ED.3) THEN 

T=- >~L - (DT IDA) ).I(AL 

I~= 0.5* ( X O~) H ) 

SLlll(I\)::. (X(K)-T)/6.0 * ( «l-XL )/[11)*( O:(K)-T)/D1) 

c + 4.0*«W-XL)/[lT)*«X(~)-W)/DT) 

HF:Fl..:: \, '~\)-T)/6.0 * " 1,\T-XU/[lTl**2 

I:; + 4.0*«(~J-XU/[lT)**2) + 1.0 ) 

89 
001 ~,3C 

00154 

00t55C 

001 ~d~ 

00157 

001~B 

00159 

00160 

00161 

001b2C 

00163 

00164C 

00165 

00166 

001b7C 

0016B 

00t69C 

00170 

OOlllC 

0017:2 

00173 

00174 

00115C 

00176 

0017lC 

00178 

00119C 

00180 

00181 

0018:2 

00183 

00184 

00185 

00186C 

00187 

00188 

00189C 

00190 

00191C 

001 S'2 

00193C 

00194 

00195C 

00196 

00197 

001 S'8C 

00199 

00:200C 

00201 

00:20:2 

00203 

ELSE IF (Jl.EQ.:2) THEN 

T= XL - (liT IliA> * AL 

W= 0.5*( TtXL ) 

SUB(" ):. <T -xu /6.0 * ( 4.0* ( (W-XU lIlT) * ( (X (K )-W) /I1T) 

C + « T-XUlDT>*( (X(")-T)/[lT) ) 

GLEF1= (T-XL)/6.0 * ( 4.0*«(W-XL)/DT)**:2) 

C t «T-XL)/DT)**2 ) 

ELSE IF (J1.E(~.l) THEN 

SUB(")= 0.0 

GLEFT= 0.0 

END IF 

ELSE IF (K.EQ.l) THEN 

SUB( 1)= 0.0 

GLEFT= 0.0 

IF (J1.EO.4) GLEFT= Dl/3.0 

END IF 

END IF 

IF (" .NE,!) 

IF ~ D(K-l) 

SLlB(K) 

GLEFT 

END IF 

END IF 

[IT= XF:-X (K) 

DA= AR-ACK) 

THEN 

.En. 0.0 ) THEN 

= 0.0 

= 0.0 

IF ( I[I(Ktl) .EO. 1 ) THEN 

IF (K.NE.N) THEN 

IF (J2.EO.l) THEN 

SUP(,,)= DT/6.0 

GF: I GH= [IT /3. 0 

ELSE IF (J:2.EO.:2) THEN 

T= X(Kl-(DT/DA)*A(K) 

W= 0.5*( XRH ) 

SUP(,,)= (XR-ll/6.0 * ( «T-X(K»/DT)*«XR-T)/DT) 
nO.-'04 

(IO~O~I 

0020t. 

OO~07(, 

OO~08 

00209C 

00210 

00~11 

()O~l':! 

00213 

00~J4 

00215 

(lO~16C 

OO~17 

00218C 

()0~19 

00220 

002:11 C 

()0~22 

00223C 

002:)4 

00225C 

00.2:!6 

00227 

002:')8 

00229C 

()0230 

0\;:~31 C 

00.22;.2 

OO~33C 

00234 

00~35 

00236C 

00237 

00238C 

00.239 

00240C 

OO~41 

00~4~C 

00243 

00244 

()0~45C 

00.24t> 

00~47C 

00248 

00249 

()02~!0 

00251 

OO~52 

00253 

002~14C 

c 

br-:IGH:: (Xr:-l )/(,.0 * ( (OJi -1 )/1IT)**2 

C + 4.0*«(XR-W)/DT)r*2) ) 

ELSE IF (J':!.ED.3) THEN 

T= X(K)-(nT/nA)*A(~) 

w~ 0.5*( T+X(~) ) 

SUP(~)= (T-X(~))/6.0 * ( 4.0l«W-X(K»/DT)*«XR-W)/DT) 

C + «f-X(K»/Df)l«XR-T)/Uf) ) 

GRIGH= (T-X(~)/6.0 * ( 1.0 + 4.0*«(XR-W)lDT)**2) 

C + «XR-T)!DT)l*2 ) 

ELSE IF (J2.EO.4) THEN 

SUP(~):. 0.0 

GF:IGH= 0.0 

END IF 

ELSE IF (K.EO.N) THEN 

SllF' ( N ) = 0.0 

GRIGH= 0.0 

IF (J2.EQ.1) GRIGH= DT/3.0 

ENII IF 

ELSE IF ( ID(Kt1) .EO. 0 ) THEN 

SUP(K)= DT/6.0 

GRIGH= DT/3.0 

ELSE IF ( IIt(K+1) .EQ. -1 ) THEN 

IF (')2.E(~.4) THEN 

SUP(K)= DT/6.0 

GF:IGH= DT/3.0 

ELSE rr (J~.EQ.3) THEN 

T= X(K)-(DT/DA)*A(K) 

W= 0.5* ( XF:tT ) SUP(~)= (XR-T)/6.0 * ( (CT-X(K»/DT)*«XR-T)/DT) 

C + 4.0*«(W-X(~»/DT)*«(XR-W)!DT) 

GRIGH= (XR-T)/6.0 * ( «XR-T)/UT1**2 

C + 4.0*((XR-W)/UT)**2) ) 

91 
OO~5~, 

00256C 

OO~~,·? 

00258 

002~}9 

00~60 

00261 

00262 

002t)3C 

00::'64 

OO~65C 

0026f, 

00267 

00268C 

00269 

OO~70C 

00~71 

OO~72C 

00273 

00274 

00~75 

00276C 

00277 

00278C 

00279 

00~80C 

OO~Ol 

0028~ 

002E.l3 

00284 

00285 

00286 

00287C 

00288 

00289C 

00290C 

00291 

00292 

00293 

00294 

00295 

00296C 

002S'7 125 

00298C 

0025'9 

00300 

00301 150 

00302C 

()0303C 

()0304C 

ELSE IF (J2.ED.2) THEN 

T= X(K)-(D1/DA'*A(K) 

W= 0.5*( T+X(~) ) 

SUP(~'~ (T-X(~')/6.0 * ( 4.0*«W-X(K»/DT'*«XR-W)!DT) 

C + «T-X(K»/DT'*«XR-T)/DT) ) 

GRIGH= (1-X(~)/6.0 * ( 1.0 + 4.0*«(XR-W)/D1)**2) 

C + (tXR-T)/DT)**2 ) 

ELSE IF (J2.ED.1) 1HEN 

SlIP(K'= 0.0 

GRIGH= 0.0 

END IF 

ELSE IF (K.EO.N) THEN 

SUP(N)= 0.0 

GRIGH= 0.0 

IF (J2.EQ.4) GRIGH= D1/3.0 

END IF 

END IF 

IF (K.NE.N) THEN 

IF ( D(~+l ) .EO. 0.0 

SUP(~)= 0.0 

GF:IGH = 0.0 

END IF 

END IF 

DJAG(~)= GLEFTtGRIGH 

IF ( tICK) .EO. 

DIAG(K)= 

SLJB(K)= 

SUPOd= 

END IF 

CONTINUE 

DO 150 L=l,N 

H(U= D(L) 

CONTINUE 

0.0 ) 

1.0 

0.0 

0.0 

) THEN 

THEN 

WE SOLVE THE MATRIX EOUATION JX=H, THE ARRAY (H) 

92 
O()30~jC 

00306C 

00307e 

00308C 

00309 

00310C 

00311 

00312 

00313 

00314 

00315 200 

00.316C 

00317C 

00318C 

00319C 

00320C 

003:! 1 C 

00322C 

o032~~ 

00324 

00325 

00326 

00327 250 

00328 

003:'9 

00330 300 

00331C 

00332C 

00333 

00334 

00335 

00336 

00337 

00338 

00339 

00340 350 

00341 

00342 400 

00343 

00344 

0034~ 

BEING IDENTICAL TO THE ARRAY (D). THE SOLUTION 

IS RlTUHNED IN THE ARRAY (H). 

CALL TRID(SUB,DIAG,SUP,H,N) 

SLIM1= 0.0 

DO 200 L=l,N 

tHL>= H(l) 

SUM1= SUM1 + ABS(A(L» 

CONTINUE 

THE FUNCTION EVALUATION SUBROUTINE COMPUT MAY 

BE DELETED. IN THIS CASE THE FOLLOWING EIGHT 

LINES ARE TO BE TlELETED AND T:1E ARRAY (FX) 

CAN BE TAKEN FROM THE REAL STATEMENT AT THE 

BEGINNING OF lHIS SUBROUTINE. 

tALL COMPUT(A,FX,N,X,TL,TR) 

FNORM1= 0.0 

DO 250 L=1,N 

FNORMl= FNORHl + FX(L)*FX(L) 

CONTINUE 

FNORM1= SQRT(FNORMl) 

WRITE(6,300) lJ,FNORM1 

FORMAT(I5,E15.6) 

IF (LJ.NE.l) THEN 

RATIO= ABS(SUHl-SUM2) 

AB= EPS*SUM2 

IFLAG= 1 

IF (RATIO .LE. AB) GO TO 400 

END IF 

SUM2= SlIM1 

CONTINUE 

IFLAG= 2 

CONTINUE 

ITMAX= LJ 

RE TURN 

END 

93 
00001 

00002C 

00003C 

OOOO'~C 

()OOO~C 

0()006C 

()0007 

00008 

00009 

00010 

00011 

00012C 

00013 

00014C 

000]5 

00016 

()00J7 

00018 

00019 

00020 

000:) 1 

00022e 

00023 

00024 

00025 

00026 

00027 

00028 

()0029 

00030C 

00031C 

00032C 

00033 

00034 

00035 

00036 

00037C 

00038 

00039 

00040 

00041 

00042C 

00043 

00044C 

00045 

00046C 

00047 

00048C 

00049 

00050C 

SUBROU1INE COMPUTCA.FX,N,X,TL,lR) 

SUBROUTINE CCOMPU1), THE FUNC1ION EVAULATJNG 

SUBROUTINE, IS OPTIONAL. 

~EAl ACN),FXCN),Fl,ALO,AHI,TLO,THI 

REAL GLEF,GRIG,TS,X(N) 

INTEGER N,K,Jl ,J2 

COMMON D(50),IDC50) 

[10 100 K=l,N 

IF ( D(~) .NE. 0.0 ) THEN 

IF (K.EQ.l> THEN 

ALO:::: 0.0 

llO= Tl 

ELSE 

END IF 

ALO=- ACK-l) 

TlO= X(K-l) 

IF (K.EO.N) THEN 

AHI= 0.0 

THI:::: TR 

ELSE 

END IF 

AHI:::: A(K+1) 

THI= X(KtU 

IF CALO.GE.O.O .AND. A(K).GE.O.O) Jl= 1 

IF (ALO.lT.O.O .AND. A(K).GE.O.O) Jl= :2 

IF CALO.GE.O.O .AND. A(K).LT.O.O) Jl:::: 3 

IF (ALO.LT.O.O .AND. ACK).LT.O.O) Jl= 4 

IF (A(K).GE.O.O .AND. 

IF (A(K).lT.O.O .AND. 

IF (ACK).GE.O.O .AND. 

IF (ACK).LT.O.O .AND. 

DT= X(~)-TLO 

AHI.GE.O.O) 

AH 1. GE • 0 .0) 

AHI.lT.O.O) 

AH 1. l T .0.0 ) 

IF ( IDCK) .EO. 1 ) THEN 

IF (Jl.EO.1) THEN 

GlEF= nT*( 2.0*A(K) t AlO )/6.0 

J2= 1 

J2= :2 

J2= 3 

. .12= 4 

94 
00051 

00052C 

000~,3 

00054 

OOO~J5 

00056C 

00057 

00058C 

00059 

00060 

OOOl)1 C 

00062 

00063C 

00064 

00065C 

00066 

00067C 

00068 

00069C 

00070 

00071C 

00072 

00073C 

00074 

00075C 

00076 

ooonc 

00078 

00079C 

00080 

00081 

00082 

00083C 

00084 

000B5C 

00086 

00087 

00088C 

00089 

00090C 

00091 

ooonc 

00093 

00094C 

00095 

00096C 

00097 

00098C 

00099 

OOlOOC 

00101 

ELSE IF (Jl.EQ.~) THEN 

TS= TLO - ALO*D1/( AC~)-ALO ) 

Fl= ( (TStX(K»lO.5 - TLO )/Df 

GLEF= (X(K)-TS)*A(K)*( ~.o*rl t 1.0 )/6.0 

ELSE IF (Jl.EQ.3) THEN 

TS= TLO - AlOtDT/( ACK)-ALO ) 

GLEF= ( (TS-TLO)**~ )*AlO/(6.0*DT) 

ELSE IF (Jl.Ea.4) THEN 

GLEF= 0.0 

END IF 

ELSE IF ( ID(K) .Ea. 0 ) THEN 

GLEF= DT*( ~.O*A(K) t ALO )/6.0 

ELSE IF ( ID(K) .EO. -1 ) THEN 

IF (Jl.Ea.4) THEN 

GLEF= DT*( ~.O*A(K) + ALO )/6.0 

ELSE IF (Jl.Ea.3) THEN 

1S= TLO - ALO*DT/( ACK)-ALO ) 

Fl= ( (TStX(K»*0.5 

- TLO )/DT 

GLEF= (XCK)-TS)*A(K)*( ~.O*Fl t 1.0 )/6.0 

ELSE IF (Jl.Ea.2) THEN 

TS= TLO - ALO*DT/( ACK)-ALO ) 

GLEF= ( (TS-TLO)**~ '*ALO/(6.0*DT) 

ELSE IF (Jl.EQ.l) THEN 

GLH= 0.0 

END IF 

END IF 

DT= THI-X(K) 

IF ( ID(Ktl) .Ea. 1 ) THEN 

IF (J2.Ea.l) THEN 

95 
00102C 

00103 

OOl04C 

00105 

00106C 

00107 

00108 

00109C 

OOll 0 

OOlllC 

0011::! 

00113 

00114 

00115C 

00116 

001l7C 

00118 

00119C 

00120 

00121C 

001:)::! 

00123C 

00124 

00125C 

001:?6 

00127C 

00128 

00129C 

00130 

00131C 

0013~ 

00133C 

00134 

0013::: 

00136C 

00137 

00138C 

00139 

00140 

00141 

00142C 

00143 

00144C 

00145 

00146C 

00147 

00148C 

00149 

00150C 

(1015]C 

0015~ 

GRIG= DT*( 2.0*A(K) t AHI )/6.0 

ELSE IF (J2.EO.~) THEN 

1S= X(K) - ACK)*DT/( AHI-A(K) ) 

GRIG= ( (THI-TS)**2 )*AHI/(6.0*DT) 

ELSE IF (J2.EO.3) THEN 

1S= X(K) - AC~)*DT/C AHI-A(K) ) 

F1= THI-0.5t( TStXCK) ) )/DT 

GRIG= ( TS-XCK) )*A(K)*( 1.0 t 2.0*F1 )/6.0 

ELSE IF (J2.EO.4) THEN 

GRIG= 0.0 

END IF 

ELSE IF ( IDCKt]) .EO. 0 ) THEN 

GRIG= DT*C 2.0*ACK) t AHI )/6.0 

ELSE IF C ID(Kt1) .EO. -1 ) THEN 

IF (J2.EQ.4) THEN 

GRIG= DT*( 2.0*ACK) t AHI )/6.0 

ELSE IF (J:?EO.3) THEN 

lS= X(K) - A(~)*DT/( AHI-ACK) ) 

GRIG= C CTHI-TS)**~ )*AHI/C6.0*DT) 

ELSE IF CJ2.EQ.2) THEN 

TS= X(K) - ACK)*DT/( AHI-A(K) ) 

F1= ( THI-0.5*( TStXOO ) )/IIT 

GRIG= ( TS-X(K) )*A(K)*( 1.0 t 2.0*F1 )/6.0 

ELSE IF (J2.EO.l) THEN 

GRIG= 0.0 

END IF 

END IF 

I F it,. HE • 1) 1 HEN 

96 
00153 

0()1~4 

0015se 

0015b 

00157 

00158 

00159C 

00160 

00161C 

00162 

00163 

001b4 

00165C 

00166 100 

00167 

00168 

IF ( [1(1\-1) .EeL 0.0 ) GLEF= 0.0 

END IF 

IF (t(.NE.N) THEN 

IF ( [1(1\+1) .EQ. 0.0 ) GRIG= 0.0 

EtHi IF 

FX(K)= GLEF + GRIG - D(K) 

~LSE IF ( DCh) .EQ. 0.0 ) THEN 

FX(K)= 0.0 

E.ND IF 

CONTItWE 

RETUF:N 

END 

97 
98 

00001 

00002C 

00003C 

00004C 

OOOO:5C 

00006C 

00007C 

ooooec 

00009C 

OOOlOC 

(,001lC 

000l2C 

OOOl3C 

OOOl4C 

00015C 

00016C 

00017C 

000l8C 

00019C 

00020C 

OOO:llC 

0002:?C 

00023C 

00024C 

00025C 

00026C 

00027C 

00028 

00029 

00030 

00031 

00032 

00033 

00034 

00035 

00036 

00037 

00038C 

00039 

00040 

00041 

00042 

00043 

00044 

00045 

00046 

00047 

00048C 

00049C 

OOO~JOC 

SU8ROU1INE POL)CA,T,PP,M,F,LI.1X) 

SUBROUTINE POLY INTEGRATES BAC~ TWICE THE 

POSITIVE PART OF THE PIECEWISE LINEAR SECOND 

DERHJATl VE WHERE THE DATA SliGGESTS THAT THE 

INTERPOLATING CURVE SHOULD BE CONVEX, THE 

NEGATIVE PART OF THE PIECEWISE LINEAR SECOND 

DERIVAfIVE WHERE THE DATA SUGGESTS THAT THE 

INTERPOLATING CURVE SHOULD BE CONCAVE, AND 

THE REMAINING PORfrON OF THE PIECEWISE LINEAR 

SECDND DERIVATIVE ON THE TRANSI1ION INTERVALS. 

lHE INTEGRATION YIELDS A PIECEWISE CUBIC 

POLYNOMIAL WIfH KNOTS GIVEN BY THE SEQUENCE 

(lX). THIS CUBIC POLYNOM1AL INTERPOLATES THE 

DATA AND ITS COEFFICIENTS ARE DENOTED BY THE 

NUMBERS PP(J,I) - THE VALUE OF THE (J-1)ST 

DERIVATIVE OF THE FUNCTION EVALUATED AT TX(I). 

FOR X SUCH THAT TX(I).GE.X.LT.TX(Itl) THE VALUE 

OF THE CUBIC POLYNOMIAL IS 

PP(!,I) 

+ PP( 2,1) * 

+ (112 ) PP ( 3 , I) * 

+ (116) PP ( 4. 1) * 

INTEGER M,J,L,LI 

X-TXCI) 

X-TX<I) )*'*2 

X - TX ( I ) ) **3 

REAL A(50),T(50),PP(4,100),F(50),TX(100),TAU 

REAL DF,Dl,DA,C,E 

COMMON D(50),ID(50) 

LI= ] 

MN1= M-l 

DO 100 L=l,MNl 

DF= F(Lt1)-F(U 

DT= T(U1)-T(U 

DA= A(Ltl)-A(U 

JP= 0 

IF CL.EQ.1> THEN 

IF ( D(1) .EQ. 0.0 ) JP= 1 

ELSE IF (L.EG.MN1) THEN 

IF ( [I(M-2) .EQ. 0.0 ) JP= 1 

ELSE 

Co:: D(L-i)iD(U 

IF ( C .EQ. 0.0 ) JP= 1 

END IF 
\ 

\ 

00051 

00052C 

00053 

00054 

00055 

00056 

00057 

00058 

00059C 

00060 

00061C 

OOOb2 

00063 

00064 

00065 

00066C 

00067 

00068C 

00069 

00070C 

00071 

00072 

00073 

00074 

00075 

00076 

00077 

00078C 

00079 

00080C 

00081 

00082 

00083 

00084 

00085 

00086 

00087 

00088 

00089 

00090 

00091 

00092 

00093 

00094C 

000S)5 

00096C 

00097 

00098 

00099 

00100 

00101 

IF (JP. EO. 1) r HEN 

F'P(4,LI)= 0.0 

F'P(3,LI)= 0.0 

F'F'(2,LI)= ItF/DT 

F'F' ( 1 , LI ) = F (L ) 

TX(U)= T(L) 

LI= LIt1 

ELSE IF (JP.EQ.O) THEN 

IF ( A (l) • GE • 0 • 0 .AND. 

IF (M(l).LT .0.0 .AND. 

IF ( A (U • 6T • 0 • 0 .MW. 

IF (M(L> .LLO.O .AND. 

A(Ltl> .GE.O.O) 

A(ltl) .GT .0.0) 

A(Ltl>.LT .0.0) 

A( Ltl> .LE .0. 0) 

IF ( II1(L> .Ea. 1) THEN 

IF U.EO.1 ) lHEN 

c= DF/DT - (DA/6.0 + A(L)/2.0)*DT 

F'F' ( 4 , LI )::: DM/D r 

F'F' (3, LI ) = A (U 

F'F' (2, LI ) = C 

PF' (1 , LI ) = F ( L> 

TXeLI)= T(U 

LI=LItl 

ELSE IF (J.Ea.2) THEN 

TAU= Tel) - AeL)*DT/DA 

C= DFIDT - (ACLt1)**3)*DT/(6.0*DM*DA) 

F'F'(4,LI)= 0.0 

PF'(3,LI)= 0.0 

PP(2,LI)= C 

F'F' ( 1 , LI ) = F (U 

F'P(4,Llt1)= DA/DT 

F'F'(3.LItl)= 0.0 

F'P(2,LItl)= C 

PF'(1,LIt1)= C*(TAU-T(L» t F(l) 

neLl)= T(U 

TX ( LI tl ) == TAU 

LI== LI+2 

ELSE IF (J.Ea.3) THEN 

TAU= T(l) - A(L)*DT/DA 

E= F(l) - (ACL)**3)*DT*DT/C6.0*DA*DA) 

c= DFIDT + (ACL)**3)*DT/C6.0*DA*DA) 

PP(4.LI)= [I,;/DT 

F'F' ( 3 , LI ) = A ( L> 

99 

J= 1 

J-- ~ ~ 

J= 3 

J= 4 
00102 

00103 

00104 

00105 

00106 

00107 

00108 

00109 

00110 

OOlllC 

00112 

00113C 

00114 

00115 

00116 

00117 

00118 

00119 

001~OC 

00121 

00122C 

001:>3 

00124C 

00125 

00126 

00127 

00128 

001:)9 

00130 

00131 

0013:!C 

00133 

00t34C 

0013~ 

00136C 

00137 

00138 

00139 

00140 

00141 

00142 

00143 

00144C 

00145 

00146C 

00147 

00148 

00149 

00150 

00151 

00152 

PP(2,LI)= C + ~(l)*A(l)*DT*0.5/DA 

F'P (1 , LI ) = F ( U 

F'F' ( 4, LI + 1) = 0.0 

PP ( 3 , LI t1 ) = 0.0 

F'F'(;:?,l1tl)= C 

PP(1,lIt1)= C*( TAU-T(l) ) + E 

TX(U)= T(U 

lXClI+l>= TAU 

LI= LIt2 

ELSE IF (J.EO.4) THEN 

PF'C4,LI)= 0.0 

PP(3,lD= 0.0 

F'F' C~, LJ ) = DF lIlT 

F'F'(l,LI)= F(L) 

lX(Ll)= T(U 

LI=LI+1 

EN[I IF 

ELSE IF ( ID(L) .EO. 0 ) THEN 

c= DF/DT - (DA/6.0 t A(l)/2.0)*I1T 

F'P ( 4 , l I ) = DIVDT 

PP ( 3 , U ) = PI ( U 

F'P (2, LI ) = C 

F'P (1 , LI ) = F ': U 

TX(LI)= T(U 

LI= LI+1 

ELSE IF ( ID(U .EO. -1 ) THEN 

IF (J.EO.4) THEN 

C= DF/DT - (DA/6.0 t A(l)/2.0)*DT 

H'( 4.LI) = DA/DT 

F'P (3, LI) = A (U 

F'P ( 2 , LI ) = C 

PP (1 , LI ) = F ( U 

TX(LI)= T(U 

LI= lIt1 

ELSE IF (J.EO.3) THEN 

TAU= T(l) - A(l)*DT/DA 

C= DF/DT - (A(Lt1)**3)*DT/(6.0*DA*DA) 

PP(4,LI)= 0.0 

PP(3,LI)= 0.0 

PP(2,LI)= C 

F'P ( 1 , LI ) = F ( U 

100 
001 ~)3 

00154 

(l01~::=; 

00156 

001~,7 

00158 

00159 

00160C 

(}0161 

00162C 

00163 

00164 

00165 

00166 

00167 

00168 

00169 

00170 

00171 

00172 

0017~ 

00174 

00175 

00176 

001i'7C 

00178 

00179C 

00180 

00181 

00182 

00183 

00184 

00185 

00186C 

00187 

001SSC 

00189 

00190C 

001S'1 

00192C 

00193 100 

00194 

00195 

00196 

00197 

00198 

00199 

00200 

PP(4,LJtl)= DA/DT 

PF' (3, LI + 1) = 0.0 

F'F' ( :> , LJ + 1) = C 

PP(l,lItl)= C'(TAU-T(l» t F(l) 

lX(U)::. T(U 

TX(Ut1>= TAU 

LI= lH2 

ELSE IF (J.EQ.2) THEN 

TAU= T(l) - A(lJtDT!DA 

E= F(l) - (A(l)**3'*DT*DT/(6.0*DA*DA) 

C= DF/DT t (A(l'*l3)lDT/(6.0*DAlDA) 

F'P ( 4 , LI ) = DA/D r 

- PP(3,lI)= A(l) 

PP(2,lI)= C t A(l'*A(l'*DT*0.5/DA 

PF-' ( 1 , LI ) = F ( l) 

PP(4,UtU= 0.0 

PP(3,LItl)= 0.0 

PP (2, l It 1) = C 

PP(1,lIt1)= C*< TAU-T(l) ) t E 

TX(U)= r(L) 

1 X < LI + 1) = TALI 

U= LIt:? 

ELSE IF (J.EQ.l) THEN 

PP(4,LI)= 0.0 

PP(3,LI)= 0.0 

PP(2,LI)= DF/DT 

PP ( 1 , LI ) = F (l ) 

TX ( LI ) = T ( U 

LI=LItl 

END IF 

ENT! IF 

END IF 

CONTINUE 

PP(4,LI)= 0.0 

PF'(3,U)= 0.0 

PP(2,LI)= 0.0 

F-'P ( 1 , LI j = F ( M ) 

TXeLl)= r<M) 

RETURN 

END 

101 
102 

Subroutlnes TRID and DATAFL are 11sted In Appendlx A. 
1 Report No I 2 Government Accession No 

NASA CR-178001 

4 Title and Subtitle 

Minimal Norm Constrained Interpolation 

7 Author(s) 

Larry Dean Irvine 

9 Performing Organization Name and Address 

Old Dominion University 

Norfolk, Virginia 

12 Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Washington, DC 20546 

15 Su pplementary Notes 

3 Recipient's Catalog No 

5 Report Date 

September 1985 

6 Performing Organization Code 

8 Performing Organization Report No 

10 Work Unit No 

11 Contract or Grant No 

NAG-1-363 

13 Type of Report and Penod Covered 

Contractor Report 

14 Sponsoring Agency Code 

505-31-83-02 

This report is a dissertation submitted to Old Dominion University for partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computational 

and Appl ied Mathematics. Langley Technical Monitor: Robert E. Smlth, Jr. 

16 Abstract 

In computational fluid dynamics and in CAD/CAM, a physical boundary is usually 

known only discreetly and most often must be approximated. An acceptable approximation preserves the sal ient features of the data such as convexity and concavity. 

In this dissertation, a smooth interpolant which is locally concave where the data 

are concave and IS locally convex where the data are convex IS described. The 

Interpolant is found by posing and solving a minimization problem whose solution 

is a piecewise cubic polynomial. The problem is solved indirectly by using the 

Pea no Kernal theorem to recast it into an equivalent minimization problem having 

the second derivative of the interpolant as the solution. This approach leads to 

the solution of a nonl inear system of equations. It is shown that Newtonls method 

IS an exceptionally attractive and efficient method for solving the nonlinear system of equations. Examples of shape-preserving interpolants as well as convergence 

results obtained by using Newtonls method are also shown. A FORTRAN program to compute these interpolants is listed. The problem of computing the interpolant of 

minimal norm from a convex cone in a normal dual space is also discussed. An 

extension of de Boorls work on minimal norm unconstrained interpolation is presented. 

17 Key Words (Suggested by Author(s)) 18 Dlstnbutlon Statement 

Surface definition, Interpolation, 

Spline functions 

19 Security Oasslf (of thiS report) 20 Security Classlf (of thiS page) 

UNCLASSIFIED UNCLASSIFIED 

Unclassified -- Unl imited 

Subject Category - 64 

21 No of Pages 22 Price 

108 A06 

N·30S For sale by the National Technical Information Service. Springfield. Virginia 22161 
End of Document