Algorithmen voor het vereffenen en aanpassen van experimentele data

Algorithmen voor het vereffenen en aanpassen van

experimentele data

Citation for published version (APA):

van Ginneken, C. J. J. M. (1980). Algorithmen voor het vereffenen en aanpassen van experimentele data. (Eindhoven University of Technology : Dept of Mathematics : memorandum; Vol. 8017). Technische Hogeschool Eindhoven.

Document status and date: Gepubliceerd: 01/01/1980 Document Version:

Uitgevers PDF, ook bekend als Version of Record Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

Onderafdeling der Wiskunde

Memorandum 1980-17

september 1980

Algorithmen voor het ~ereffenen en aanpassen van experimentele data door

Technische Hogeschool Onderafdeling der Wiskunde PO Box 513, Eindhoven Nederland

1. Inleiding

2. Definities en notatie

3. Methoden van Schoenberg en Reinsch

4. Een automatische keuze voor de parameter A

5.

Numerieke experimenten 6. Procedures Literatuur

...

blz. 1 1 2 4 11 22 30

1. Inleiding

Zij gegeven de data (x.,Y.), i =. O,l, ••• ,n, met Xo < xl < ••• < x en waarbij

~ ~ n

Y

i een, bijvoorbeeld door meting verkregen, benadering is voor een onbekende functiewaarde f(x

i). Gevraagd wordt om een functie 9 te bepalen die de on-nauwkeurigheden in Y

i vereffent.

In het geval dat bijvoorbe~ld bekend is dat de gemeten functie f behoort tot een door enkele parameters gekarakteriseerde klasse is het aan te b~velen de kleinste kwadraten-aanpassing te nemen als benadering voor f. In het geval dat deze informatie ontbreekt kan met een van de te beschrijven algorithmen een aanpassing met behulp van splines worden berekend.

In paragraaf 3 worden de methoden voorgesteld door Schoenberg en Reinsch beschouwd. Deze methoden bevatten parameters waarvoor bij gebruik geschikte waarden moeten worden gekozen, wat soms moeilijk is. In paragraaf 4 wordt

een methode beschreven die sterk verwant is aan die ui t paragraaf 3, maar welke -parametervrij is.

In paragraaf 5 worden de verschillende methoden voor een aantal voorbeelden onderling vergeleken. In paragraaf 6 worden een aantal procedures beschreven. 2. Definities en notatie

In het volqende zijn (Xi,Y

i ), i = O,l, ••• ,n, n ~ 1 gegeven data met Xo < Xl < ••• < x

n' De getallen di , i = O,l/ ••• ,n zijn gegeven zogenaamde gewichtsfactoren.

Definities

Ck[Xo'XnJ is de verzameling van de functies, die k-maal differentieerbaar zijn op [xO,xnJ en waarvan de k-de afgeleide continu is;

gk[xo'X

n] is de verzameling van de functies, die k - 1 maal differentieer-baar zijn en waarvan de k-de afgeleide continu is in (x

i ,xi+1), i = O,l, ••• ,n-l. en kwadratisch integreerbaar op [xO,x_n

];-de functionalen xn E (f) :=

J

[f(m) (x)]2dx , f E H(m)[X x] m 0' n Xo n 2 Kef) :=

2

_{dix-(f(xi ) - Y i) ;} i=O

8

2m-1 is de verzamelinq van de natural spline functies van de qraad 2m--- 1 met knooppunten x_O,x₁' ••• ,x

n•

Pro memorie qeven we de definitie vaneen natural spline.

Een natural spline functie van de qraad 2m - 1 met knooppunten x

O

,x

1, ••• ,xn is een op [xO,xnJ qedefinieerde functie s met de volgende eigenschappen

(i) (ii) (iii)

op het interval [xi,x, 1J is seen polynoom van de qraad 2m - 1

2m-2 ~+

s ~

c

[xO,x J

S (k) (x

O)

=

s(k?(Xn)

=

0, k

=

m,m +1 , ••• , 2m 2 - . 3. Methoden van Schoenberg en Reinsch

Om een functie te bepalen die de onnauwkeurigheden in Y

i verdisconteert stelde Schoenberq (1964), uitqaande van een idee van Whittaker (1923), het volqende v~~r.

Zij m ~ n+1 en A > 0 qeqeven. Bepaal gA die de functionaal

minimaliseert.

-(3.1)

De parameter A is een maat voor het compromis tussen de "qladheid" van gA' die wordt qemeten met Em{qA) en de afwijking van

qA

van de data Y

i ' die wordt qemeten met de som van kwadraten K(gA>.

Schoenberq toonde het volgende aan. Er bestaat een eenduidig :bepaalde gA en deze behoort tot de klasse 8

2m-i " V~~r een bewijs zie Greville (1969). Als A ~ 0 dan qaat

qA

over in het polynoom go met graad ~ m - 1 waarvoor K minimaal is en als A + = nadert

qA

tot de natural spline die de data Yi interpoleert. Er qeldt dat K(qA) een monotoon niet stijgende en E(gA) een monotoon niet dalende functie van A is.

In practische toepassingen is het kiezen van een geschikte waarde voor A een probleem. Een geschikt~ waarde kan dan bijvoorbeeld worden vast-gesteld door een aantal waarden te proberen. Daarom stelden 5choenberg

(1964) en Reinsch (1967) ook nog de volgende aanpak v~~r.

Zij (in plaats van A) een getal 5 0 < 5 < ~ gegeven. Bepaal h5 die m

de.functionaal Em(g), 9 € H [xO,xnJ minimaliseert onder de

nevenvoor-waarde K(g) S 5.

In het geval dat bijvoorbeeld bekend is dat de data Y

i onderling onaf-afhankelijke ruiscomponenten e. bevatten met variantie

cr~,

ligt het voor

1 -2

de hand de gewichtsfactoren d

i gelijk te nemen aan cri en 5 ongeveer ge-lijk te nemen aan n + 1.

Er bestaat het volgende verband tussen de problemen (3.1) en (3.2).

Zij K(gO) := lim K(gA) en zij 0 < 5 < K(gO)' dan heeft (3.2) precies een A+O

*

oplossing h

_{s •}

Zij A uit (3.1) zodanig dat K(g~.>

=

5, dan geldt h5

=

9 .•

A' A

*

Dit is eenvoudig in te zien. Immers K{gA) is een niet stijgende functie

*

van A met 0

=

K(g~) < 5 < K(gO)' dus er m

is een A met K(g *}

=

5. Voor

A

iedere 9 € H [xO,xnJ met K(g) S 5 geldt

.

*

*

*

*

Em(cj *) + A K{g *)

=

E

(q ,)

+ A 5 S E (g) + A K(g) S E (g) + ). 5,

). A m).* m m

waaruit volgt dat E (g ) S E (g) en dus is 9 een oplossing van (3.2).

. m).* m ).*

Ook geldt voor iedere oplossing h van (3.2)

5

*

*

*

Em(h

s) + ). K(hs) S E (g ) _{m ).*} + ). S

=

E (g *) _{m A} + ). K(g *). _A Omdat de oplossing van (3.1) eenduidig is volgt hieruit dat hs

=

9 •.

).*

Als S = 0, dan is hS de natural spline die de data Y

i interpoleert en als S ~ K(gO)' dan is bijvoorbeeld hs

=

go een oplossing.

Algorithmen voor het berekenen van hs zijn gegeven door Reinsch (1967) voor het gevalm

=

2, Woodford (1970), Lyche and Schumaker (1973,1974) voor willekeurige m ~ 1. De bepaling van hsgebeurt steeds, indien

o

< S < K(gO)' door het oplossen van de van A afhankelijke vergelijking K(gA) = S. De bepaling van gAt bij gegeven

A,

komt neer op het oplossen van een stelsel lineaire vergelijkingen, waarvan de matrix een bandmatrix is met bandbreedte 2m + 1.

4. Een automatische keuze voor de parameter A

Alhoewel het in een groot aantal gevallen mogelijk is om een geschikte waarde voor Suit (3.2) aan te geven, blijven er situaties waarin dit moeilijk is. Daarom willen we nu een strategie aangeven om een geschikte waarde voor

A

uit .(3.1) te bepalen. We beperken ons hierbij tot m =

2.

De uitbreiding naar willekeurige m is echter wel mogelijk.

We beschrijven eerst in het kort hoe de oplossing van (3.1) bij gegeven A numeriek berekend kan worden. Een natural cubic spline g(x) € S3 kan

op het interval [~i,xi+1J geschreven worden als

xi₊1 - x [ 1 2 2 ] g(x} = h ' g ( x_i)

-'6

g "(xl.,)(h i+L - (x i+1 -x» + . i+~ , ~ (4.1) waarin h, L := x, 1 - x .• l.+~ l.+ l. We stellen a

i := g(xi), ci := g"(xi), i = O,l, ••• ,n. Uit de continuiteit van de eerste afgeleide in de punten x" j

=

1,2, ••• ,n-1 volgen de

vergelij-J

kingen

= a_j+1 - a_j hj+~

Omdat

Co

=

_cn

=

0, kunnen deze relaties geschreven worden als (4.2) Tc

=

Q T a,

. T T

waarin c := (c

1'c2, ••• ,cn_1) I a :={aO,a1, ••• ,an} , Teen positief defi-niete tridiagonale (n - 1) x (n - 1) matrix is en Q een tridiagonale

(n + 1) x (n - 1) matrix 1 T==-6 2 _(h1/2 + h3/2 ) h3/2

...

Q== x-x. l. g"(X) == ~--hi+~

o

o

-

_-

--1 h3/2

--

-h n-3/2

-

-o

-1 -1 ... :-h3/2 - hS/2 ...-1 -1 - _ h_{n-3/ 2} hS/2

~1

-1 ... -h_n__{3/ 2 - h} n_1/ 2 ... -1 h_n_₁/ 2

waaruit (na enig rekenwerk) volgt dat

(4.3)

₌

c Tc, T

of met (4.2)

(4.4) _{E2 (g) == a} T _{Q T} -1 T _{Q a •}

o

Zij D de diagonaalmatrix met D .. := d

i , i = O,l , ••• ,n, dan geldt

l.l.

-1 T

Omdat QT Q symmetrisch is geldt dat als (4.5) minimaal is dat dan a voldoet aan het stelsel

(4.6) (QT Q + AD)a -1 T

=

'ADy • Uit (4.2) en (4.6)- volgt het stelsel

(4.7) (Q D Q + T -1 ~T)c

=

AQ y T en de relatie

(4.8) a = y -

TD

1 -1 Qc.

Door het-stelsel (4.7) op te lossen, waarvan de matrix positief definiet is met bandbreedte 5, wordt c. = g~(x.) verkregen en vervolgens uit (4.8)

~ 1\ ~

a_i - gA(X

i). De berekening van gA(X) voor willekeurige x gaat dan eenvoudig met (4.1).

We willen nu uit de verzameling functies g" (~), 0 < A <

00,

een geschikte kandidaat zoeken, die de ruis in de gegeven data y. verdisconteert. Deze

~

kiezen we op grondvan de volgende overwegingen. Als gA(X} geschikt is, dan moet de vector

een ruissignaal voorstellen. Zij 0A(x} de natural cubic spline die r" inter-poleert, dus 0A(x

i}

=

(rA)i. Voor de "energie" E2(OX) volgt uit (4.4) (4.9)

waarin Indien van de

E2 (oA) = rAB rA, , T B := QT-1QT.

r" een ruissignaal verwachtingswaarde

is met onder ling ongecorreleerde componenten, waar-nul is voor iedere component, dan geldt voor de

(4.10)

waarin D de diagonaalmati:ix is bestaande uit de diagonaalelementen van B. Uit (4.10) destilleren we nu het volgende criterium. Probeer A zo te kiezen dat

(4.11) F()d

Het is (mij) niet gelUkt om, eventueel onder bepaalde voorwaarden, precies aan te geven of de vergelijking (4.11) wortels heeft en hoeveel. We kunnen echter weI op basis van gevoelsmatige argumenten enig inzicht geven over het kwalitatieve gedrag van F(A) in de praktijk.

Veronderstel eenvoudigheidshalve dat de punten x. equidistant zijn, dus

~ ~

x_i+₁ - x.

=

h. De elementen D .. zijn positief en in de orde van grootte

-3 ~ ~~ T~

van h • Op grond hiervan zal de orde van grootte van rADr

A gelijk zijn

T 3

aan (rArA)/h , notatie:

(4.12)

Er geldt

(r,).

=

(g,(x.) - f(x.» + (f(x.) - y.),

I\J. 1\ ~ J. J. ~

(y. is een benadering voorf(x.». Als A klein is, dan is g, bijna lineair

~ J. 1\

en dan zal in het algemeen (mits f niet bijna lineair is) de term f(x.) - y.

~ J.

verwaarloosbaar zijn t.o~~. gA(X

i ) - f(xi ). De spline 0A(X) interpoleert .dus ongeveer de functie gA (x) - f(x), welke een gladde functie is. Op grond

hiervan verwachten we dat, als A klein is, de orde van grootte van E 2(OA) gelijk is aan f.xn [f"(x)]2dx , dus

Xo

T (Xn 2

E2 (OA)

=

rABr_A", [f" (x)] dx.

o

Uit (4.12) en (4.13) voIgt dat voor niet bijna lineair is dat dan

r~BrA

groot wordt dan gaat r

A en dus F(A)

kleine waarden van A en h en als f

or....

« rADr_A en dus dat F(A) < O. Als A naar O.

We doen, op grond van het voorafgaande, de volgende keuze voor de aan-passing gA(X)

- indien F(O) ~ 0 kiezen we gA(X)

indien FCO) < 0 kiezen we gACX), waarin A1 de kleinste wortel is van (4.11); (A_l kan eventueel oneindig zijn).

Een belangrijk onderdeel bij de numerieke berekening van een nulpunt van '(4.11) is de berekening van F (A) voor een gegeven waarde van A.

Zij

en

dan geldt dat

en hieruit voIgt met behulp van (4.2) en de definitie van B = QT-lQT

(4.14)

Het getal yTBY en de vector QTy zijn onafhankelijk van A. De vector c

vol-T -1 T

doet aan (4.7), waarin de matrices Q D Q, T en de vector Q y als constanten

T

voorkomen. De berekening van rABr

A komt dus in feite neer op het oplossen van (4.7). De berekening van

...,

-qaat bij gegeven c en D eenvoudiq via (4.8). De elementen b

kk kunnen als volgt berekend worden.

Er geldt

waarin e

k, k

=

O,l, ••• ,n, de k-de eenheidsvector is. De matrix T is

positief definiet en kan daarom geschreven worden als (Cholesky splitsing)

.T=LXLT,

'dus

T T

waarin Zk kan worden berekend uit LZ

k

=

Q eke Op ana loge wijze wordt y By uit (4.14) berekend. In het geval dat de punten x. equidistant zijn kunnen

1.

...,

we de volgende formules voor Dkk afleiden

(4.15)

13

(1 + n n-1 DOO = Dnn

=

_{h 3(1 _} a - a - a ), an) 1 3613 (1 + n k n-k D = - (-48 + a. - a - a

»,

kk h3 _{1 -} .n a waarin h := (x - x )/n en a

=

7 - 4/3. n

a

k = 1,2, ••• ,n-1, -1 T

Bewijs: Er geldt B

=

QT Q, waarin de (n + 1) x (n - 1) matrix Q en de (n - 1) x (n - 1) matrix T worden gegeven door

1 -2 _~

0

4

0

\ , 1 1 -2 "

"

"- 1, 4 1

"

"

,

"-1 1, "

,

_"

_"

h

_"

_"

_"

Q

= -

,

"- _T

=

,

,

,

6

_"

_"

_"

_",

6

..

, ,

"

, ,

,

,

"

,

"1

_"

,

,

"

"

_"

0

'

"

,

"

~2

0

"

,

_"

,

,

,

,

1

" "

,

'1 '1 4

Zij nu P :=

~T

I dan geldt

1

Q=-h p - 6T , waarin a-1 rasp. e n-1

de.eerste resp. (n - 1) de eenheidsvector is. Voor de matrix B vinden we

T -1 I ~ I ' - . / . e₁ P e l l / " " ' - . . I .,... _______ L _________________ L __________ _ ·(4.16) B = -_.. 6 ₃ h I -1 I

X

I

P - 12I + 36P

I

X

_______ ~ _________________ 4 __________ _ I I T -1

><

I

X

Ie P e I I n-l n-1 I t

waaruit volgt dat de diagonaalelementen D~ van de matrix B bekend zijn als de diagonaalelementen van de matrix P bekend zij.n. We berekenen deze op de volgende manier.

Zij e

k de k-de eenheidsvector en zek)

(k) (k)

:= (z1 , ••• ,zn_1) -de oplossing van

pz(k)

=

dan geldt (4.17) (p-1) (k) kk

=

zk • Definieer weglaten (k) . (k) z := z

o

n := 0, dan volgt uit (4.17), als we de bovenindex k

(4.18.1) zi_l + 4zi + zi+1

=

0 1 s i S k-l, (4.18.2) zk_l + 4zk + zk+l = 1,

Uit (4.18.1) en

Zo

=

0 volgt

(4.19) z.

=

a(A - A ) i -i

1.

o

S i S k,

waarin

A

=

-2

~

13

en a een constante is. Evenzo volgt u1t (4.18.3) en z

=

0

n

(4.20) k S i S n,

waarin ~ een constante is •

. Uit (4.19) en (4.20) voor i

=

k en uit (4.18.2) volgen nu twee relaties waaruit a en

a

kurtnen worden opgelost. Na enig rekenwerk vinden we

-1

h

n n-k k ( P )

=

z = --...;;~- (1 + a - a -a) , kk k 6(1 _ an) .... 2 waarin a = A = 7 - 4/3 •.

Substitutie in (4.16) levert de relaties (4.15).

s.

Numerieke eXferimenten

In deze paragraaf worden, aan de hand van een aantal voorbeelden,de methoden beschreven in §3 en §4 onder ling vergeleken.

Voorbeeld 1:

Gegeven zijn de data x.

=

i x n/23, i

=

0,1, ••• ,23 en y. de op 1 decimaal

1. 1.

afgeronde waarde van sin(x

i). We nemen de gewichten di

=

1 en omdat we de ruiscomponenten in y. opvatten als random trekkingen uit een homogene

ver-1.

deling op (-0.05,0.05), waarvan de variantie 0.01/12 is, nemen we bij de methode van Reinsch S

=

24 x 0.01/12.

In figuur 1 is de grafiek getekend van de bij deze gegevens berekende natural spline h

s

.

Bovendien zijn in deze figuur de data (xi'Yi) aangegeven

. 5

In figuur 2 is de grafiek getekend van hi en in figuur 3 de grafiek van

s

h" • De aanpassing berekend met de methode uit §4 levert, tot op

teken-s

nauwkeurigheid, identieke resultaten. Opvallend is het minder gladde karakter van hI! in figuur 3. Oit is mijns inziens inherent aan de methode

S

omdat uit formule (4.8) blijkt dat de correctie op y, welke de ruis in y zou moeten representeren, gelijk is aan

!

0-lQc , met c

i

=

h~(xi)

•

O~---~---~

o

n/2

fIG.

1 .5

o

-.5 -1

o

!

-.5 -1

o

o

Voorbeeld 2: 1\/2 FIG. 2 HI2 FI G. 3

Gegeven zijn de data x.

=

i, i

=

0,1, ••• ,23 en y. random trekkingen uit een

~ ~

homogene verdeling op (0,1), waarvan de variantie 1/12 is. We nemen d. _~ = 1 en s~= 2 bij de methode van Reinsch. De bij deze gegevens berekende natural

spline hs levert in figuur 4 getrokken rechte lijn OPe De data (xi'Yi) zijn H

aangegeven met een x. De aanpassing berekend met de methode uit §4 levert identieke resultaten. 1 x

_x

x

x

x

x x

x

_x

.5

l~~x~~

____

---~x~---;x~~-x;---x:-~x~-r

x

x

x

x

o

o

Voorbeeld 3:

x

x 10 FIG. 4

x

•

.

20 1

Gegeven zijn de_data

Xo

=

0, xi

=

(i -

2'

+ r

i) x 11'/23, i

=

1,2, ••• ,22, x₂₃ - ~ en Y

i

=

sin(xi) + 0.1 x si' waarin ri en si random trekkingen zijn uit een homogene verde ling op ~-0.5, 0.5). We nemen d

i

=

1 en S = 24 x 0.01/12. Infiguur 5 is de gra~iek getekend van de berekende spline hS en bovendien zijn de data (Xi,Y

i) gemarkeerd door een x. In figuur 6 resp. 7 is hs resp. h" getekend. De figuren 8,9 en 10 zijn analoog aan 5,6 en 7 maar nu voor

S

.5

o

.5 -1 •

x

x

• 5

OA---__________

~

________________________________

~

o

o

nn

FIG. 5 nl2 FIG. 6

n

o

~~---~ -.5 -1 .0

nl2

F JO. 7 .5 OX---~---~x

o

n/2

n

FlO. 8

1

o

-1

o

o

-.5 -1

o

tt/2 FIG. 9 R/2 FIG.tO • 1\

De punten Xi zijn als in voorbeeld 3 en Y

i is de op 1 decimaal afgeronde

waarde van sin(x.). We nemen d. en S weer als in voorbeeld 3.

l. l. .

In figuur 11 is de grafiek getekend van de berekende spline hs en de

data (x.,y.) zijn gemarkeerd door een x en in figuur 12 is de spline weergegeven

l. l.

berekend met de methode van §4. Opvallend is het minder gladde karakter van de spline uit figuur 12. Een mogelijke verklaring hiervoor zqu kunnen zijn dat door het afronden op 1 decimaal niet zo goed voldaan is aan de veronderstelling dat de ruiscomponenten in y. onderling onafhankelijk zijn.

1 .5 O~---~---~

o

·5 lt/2 fIO.ll o~---~---~

o

R/2 FJO.12 Voorbeeld 5:

Gegeven zijn de data x. _~

=

i _. x ~/n, i

=

O,l, ••• ,n en Yi de op d decimalen -2d afgeronde waarde van sin(x.). We nemen d. = 1 en S

=

(n + 1) x 10 /12. Zij -s - sl -s _{. 2} -s3 ~ ~

de natural spline die de functie sin (x) interpoleert in de punten x. ~ de natural spline die de data y. interpoleert

~

de spline berekend met de methode van Reinsch uit §4 de spline berekend met de methode ui t § 4 •

Zij r!P) de als volqt qedefinieerde maat voor de afwijking van de p-de af~eleide van s en s: ";J m ' r(p):= 1 JD n

+

1

Ir

i=O

In de volgende drie tabellen zijn voor n = 23, 45, 90 en d = 1,2, ••• ,6 de afwijkingen r(p)getabelleerd.

m

rr'abel 1 n.= 23.

Aarttal Orde van Interpolatie Reinsch r1ethode §4

!decimalen d de afqeleide p _{rep) r(p) r(p)}

-. • - - . * - . -

-

' -. -. '"" ... '" -"--- 1 2 3 1 0 2.5 10-2 2.310-2 2.310-2 1 2.5 10-1 6.110-2 6.210-2 2 8.1 9.4 10-2 9.510-2 2 0 1.9₁₀-3 2.7 10-3 . 1.610-3 1 2.1 10-2 1.110-2 5.610

-3

2 4.5 10-1 2.510-2 2.410-2 3 0 2.3 10-4 2.210-4 4.010-4 1 2.0 10-3 1. \0-3 2.610-3 2 1.1 10-1 1.210-2 1.210-2 4 0 2.4 10-5 3.110-5 3.110-5 1 1.9 10-4 2.010-4 2.010-4

.

" 2 7.0 10-3 3.010-3 3.010-3 5 0 2.9₁₀-6 3.9 10-6 2.910-6 1 2.9₁₀-5 3.2 10-5 2.910-5 # 2 6.5 10-4 6.310-4 6.510-4 6 0 1.6 10-7 3.210-7 1.610-7 1 1.6₁₀-6 2.1₁₀-6 1.4₁₀-6 2 4.7 10-5 4.910-5 4.6 10-5

Tabel 2 n ... 45

Aantal Orde van Interpolatie Reinsch Methode §4

decimalen d de afqeleide P _rep) rep) rep)

1 2 3 1 0 2.7 10-2 1.210-2 2.110-2 1 5.2 10 -1 4.9 10-2 6.310-2 2 4.2 10 1 1.0 10-1 1.110-1 2 0 2.6 10-3 2.210-3 1.410-3 1 3.8 10-2 9.010-3 6.410-3 2 2.7 3.0 10-2 4.710-2 3 0 2.4 10-4 2.110-4 2.110-4 1 4.8 10-3 1.610-3 1.610-3 2 3.0 10-1 1.110-2 1.110-2 4 0 2.7 10-5 3.310-5 4.810-5 1 4.0 10-4 3.210-4 4.810-4 2 3.1 10-2 3.410-3 3.410-3 5 0 2.9 10-6 3.710-6 4.310-6 1 5.7 10-5 5.710-5 6.710-5 2 2.8 10-3 1.410-3 1.310-3 6 0 3.0 10-7 3.910-7 3.010-7 1 5.6₁₀-6 9.0₁₀-6 5.8₁₀-6 2 3.2 10-4 3.210-4 3.310-4

.

Tabel 3 n = 90

Aantal Orde van Interpolatie Reinsch Methode §4

decimalen d de.afgeleide p r(p) rep) rep)

1 2 3 1 0 2.7 10-2 8.410-3 7.710-3 1 1.1 3.6 10-2 5.710-2 2 1.2 102 7.910-2 9.510-1 2 0 2.6 10-3 9.010-4 1.910-3 1 8.0 10-2 5.310-2 9.610-3 2 1.7 101 2.310-2 3.010-2 3 0 2.5 -4 10 2.110-4 1.710-4 .1 8.9 10-3 1.610-3 1.310-3

-2 1.2 9.9 10-3 1.010-2 4 0 3.0₁₀-5 2.8 10-5 3.910-5 1 8.5 10-4 3.610-4 5.110-4 2 1.4 10-1 4.010-3 4.410-3 5 0 2.7 10-6 1.710-6 1.610-6 1 _{9.4 10-:5} 5.9₁₀-5 8.4₁₀-5 2 1.3 10-2 1.610-3 -3 1.5 10 6 0 3.0 10-7 3.510-7 2.110-7 J 1 1.1 10-5 1.410-5 1.410-5 2 1.2₁₀-3 6.7 10-4 6.610-4

Uit de gegeven numerieke resultaten mege blijken dat de kwaliteit van de methode uit §4 vergelijkbaar is met de methode van Reinsch. Oak vele andere numerieke experimenten hebben dit bevestigd.

6. Procedures

In deze paragraaf worden procedures beschreven voor de berekening van de diverse cubic spline aanpassingen uit het voorafgaande. Er is een verzameling procedures voor het geval dat de steunpunten x. equidistant

l.

zijn. Deze procedures worden beschreven in §6.1. Een verzameling procedures voor het geval dat de steunpunten x, niet equidistant zijn wordt beschreven in

l.

96.2. In beide paragrafen worden achtereenvolgens procedures gegeven voor de be-rekening van de parameters van de cubic spline:

- gAo uit (3.1) (zie §6.1.1 of §6.2.1) - g S ui t (3. 2) (zie § 6. 1 .2 of § 6 • 2.2)

volgens de methode van §4 (zie §6 .. 1.3 of §6.2.3) - die de data Y

i interpoleert (zie §6.1.4 of §6.2.4).

Bovendien wordt een procedure beschreven o~, uitgaande van de parameters van de spline, de functiewaarde of een van de afgeleiden van de spline in een gegeven punt te berekenen (zie 96.1.5 of §6.,2.5).

6.1. Procedures .voor esuidistante punten

6.1.1. De procedure Reinscheqlambda

procedure Reinscheqlambda (n,xO,xn,y,d,lambda,a,da,d2a,term); value n,xO,xn,lambda1 integer n,term;

~ xO,xn,lambdaj ~ arra~ y,d,a,da,d2a;

Formele parameters n xO. xn Y lambda d

Bevat het rangnummer van het laatste datapunt. Voorwaarde n ~ 2. Bevat xO.

Bevat xn. Voorwaarde: xn > xO.

Array met grenzen [0 : nJ dat bij aanroep en na beeindiging de waarden y. bevat.

l.

Bevat de.waarde voor de gladstrijkparameter

A

(zie toelichting). Voorwaarde: A ~ O.

Array met grenzen [0 : nJ dat bij aanroep en na beeindiging de gewichtsfactoren d. bevat. Voorwaarde d. > O.

at da, d2a Array's met qrenzen [0 : nJ. Indien na beeindiqinq term

=

0 dan bevatten a[iJ, da[iJ, d2a[iJ de waarden qA(X_i), ql(x_i ), q~(xi) (zie toelichting).

term Bevat na beeindiging van de procedure informatie over de

af-loop van het proces.

Toelichting

Indien

term

=

0: de procedure heeft gA bepaald (zie toelichting) term = 3: gedurende de berekeningen is, ten gevolge van

afrondfouten, een theoretisch onmogelijk resultaat ontdekt

term = 4: de input voldoet niet aan de genoemde restricties n ~ 2 etc.

Zij qegeven de data (xi'Yi)' i

=

O,l, ••• ,n met equidistante Xi en posi-tieve qewichtsfactoren d

i en zij qeqeven een reeel getal A ~ O. Zij gA de lvoldoend gladde) functie die de functionaal

[g"(X)]2dx + A

minimaliseert. Dan qeldt gA is een natural cubic spline (zie §3). De procedure berekent, zo mogelijk, de waarden gA(X

i), qi(xi ) en g~(xi)' i = O,l, ••• ,n, (zie §4).

).1.2. De procedure ReinscheqS

procedure ReinscheqS (n,xO,xn,y,d,S,a,da,d2a,lambda,term); value n,xO,xn,S; integer n,term;

Formele parameters

Voor de betekenis van n,xO,xn,y en d zie 6.1.1.

s

Bevat de waarde voor de gladstrijkparameter S (zie toelichting).

Voorwaarde: S > O.

a,da,d2a Array's met grenzen [0 : nJ. Indien na beeindiging term

=

0 of term = 1 dan bevatten a[i], da[i], d2a[iJ de waarden hs(X

i), lambda

term

hg:(x_i), h~(xi) (zie toelichting).

Bevat na beeindiging indien term

=

0 de bij de waarde van S

*

behorende waarde A {zie toelicnting} en indien term

=

1 de waarde O.

Bevat na beeindiging van de procedure informatie over de afloop van het proces.

Indien

term

=

0: de procedure heeft hS bepaald (zie toelichting)

term

=

1: de waarde voor S is zo groot dat de berekende aanpassing h S lineair is (zie toelichting)

term

=

3: gedurende de berekeningen is, ten gevolge van afrond-fouten, een_theoretisch onmogelijk resultaat ontdekt term = 4: de input voldoet niet aan de genoemde restricties

n ~ 2 etc ••

Toelichtin9

Zij gegeven de data (Xi'Yi) i

=

O,l, ••• ,n met equidistante xi gewichtsfactoren d. en zij gegeven een reeel getal S > O. Zij

~

(voldoend gladde) functie die de functionaal

2

[g"(x}] dx

minimaliseert onder de voorwaarde

en positieve hs een

Uit §3 volgt

- hS is een natural cubic spline

- als S niet al te groot is (voorwaarde: S < K(gO» dan is hS eenduidig

*

en bestaat er precies ~en A waarvoor de functie 9

*,

zoals de proce-A

dure Reinscheqlambda uit 6.1.1. die aflevert, gelijk is aan hs - als S groot is (voorwaarde: S ~ K(gO» dan is de lineaire functie

hs := go een oplossing.

De procedure berekent, zo. moqelijk, de waarden hs(X_i), h~(xi)' _{hs(Xi ) en}

, *.

f

A Indien bekend is dat de data y. onderling ona hankelijke ruiscomponen-2 ~

ten &i bevatten met variantie

_.

a.

dan wordt aanbevolen d. gelijk te nemen

~ ~

aan a~2 en S ongeveer gelijk aan n + 1. 6.1.3. De procedure CvGeq

procedure CvGeq (n,xO,xn,y,d,a,da,d2a,lambda,term); value n,xO,xn; integer n,term;

~ xO,xn,lambda; ~ array y,d,a,da,d2ai

Formele parameters

Voor de betekenis van n,xO,xn,y en d zie 6.1.1.

a,da,d2a .Array's met grenzen [0 : nJ. Indien na beeindiging term

=

0,1 of 2

lambda

term.

--_

...

dan bevatten a[iJ, da[i], d2a[i] de waarden g, (x.), g~ (x.),

Al 1. (1.1 ~

g~ (Xi) (zie toelichting).

1

Bevat na beeindiging indien term

=

0,1 of 2 de waarde Al (zie toelichting).

Bevat na beeindiging van de procedure informatie over de afloop van het proces.

Indien

term

=

0,1 of 2: de procedure heeft gA .bepaald (zie verder

1

term

=

3

term = 4

toelichting)

gedurende de berekeningen is, ten gevolge van afrondfouten, een theoretisch onmogelijk resul-taat ontdekt.

de input voldoet niet aan de geeiste restricties n ~ 2, etc ••

Toelichting

Zij gegeven de data (x., y. ) I i

=

0,1, ... , n met equidistante xi en posi tieve

J. J.

gewichtsfactoren d .• De procedure berekent, zo mogelijk, (zoals

beschre-J.

ven in §4) een waarde Al en de daarbij behorende aanpassing gA zoals

1

de procedure Reinscheqlambda die aflevert. Indien na beeindiging geldt dat

- term

*

0, dan is Al de ]Ueinste wortel van (4.11) term

=

1, _{dan is Al}

=

a

en de aanpassing gA line air

. 1

- term = 2, dan is het niet gelukt de kleinste wortel van (4.11) te berekenen en is Al zodanig dat in het interval [O,A_l] geen wortel ligt.

6.1.4 De ~rocedure splineinte;E0leq

procedure splineinterpoleq (n,xO,xn,y,dy,d2y,term); value n,xO,xn; integer n,term;

~ xO,xn; ~ array y,dy,d2y

Formele parameters n

XO,~ .. xn y

_Bevat het rangnummer van het laatste datapunt. Voorwaarde: n ~ 2. Bevat xO.

Bevat x_n• Voorwaarde: xn >xo.

Array met grenzen [0 : nJ dat bij aanroep en na beeindiging de waarden y. bevat.

J.

dy,d2y Array's met grenzen [0 : nJ. Indien na beeindiging term =

a

dan

term

bevatten dy[i], d2y[i] de waarden s'(x.), s"(x.), waarin s de

. J. J.

natural cubic spline is die de data y. interpoleert in de

knoop-J.

punten xi.

Bevat na beeindiging van de procedure informatie over de afloop van het proces.

Indien

term

=

0: de procedure heeft de interpolerende natural cubic spline s bepaald

term = 3: qedurende de berekeninqen is, ten qevolqe van afrond-fouten,een theoretisch onmoqelijk resultaat ontdekt term

=

4: de input voldoet niet aan de qenoemde restricties

n ~ 2, etc ••

Toelichtinq

Zij qeqeven de data (xi,y

i), i = O,l",./n met equidistante xi' Zij s de natural cubic spline die de data Y

i interpoleert, ofwel s minimali-seert de functionaal

fX

n [q"(x)]2dx

Xo

onder de voorwaarde q(x

i) =

Y

i • De procedure berekent, zo moqelijk, de waarden s(x.), S'(X

i), s"(x.), i =O,l, ••• ,n.

1. 1.

6.1.5. De procedure splinederiveq

~ procedure splinederiveq (n,xO,xn,a,da,d2a,p,t);

value n,xO,xn,p,t; integer n,pi

~ t,xO,xni ~ array a,da,d2a;

Formele parameters n

xO

Bevat het rangnummer van het laatste knooppunt. Voorwaarde: n ;;:: 1.

Bevat xO.

xn

a,da,d2a

p

Bevat xn. Voorwaarde: xn > xO.

"

Array's met qrenzen [0 n]. Bij aanroep en na afloop bevatten a[i], da[i], d2a[iJ de waarden s(x

i), S'(Xi), slt{xiJ van de natural cubic spline s die de data a[i] interpoleert in de knooppunten xi'

Orde van de te berekenen afqeleide.

t Waarde van het argument.

Toelichting

Zij seen natural cubic spline met equidistante knooppunten xi' i

=

O,l, ••• ,n. De procedure berekent uitgaande van de waarden s(x

i), s'(x

i ), S"(Xi ), p en t de waarde s(p) (t). Opmerking

Op de intervallen (~,XO) en (x

n,=) is s lineair met afgeleide SI(XO) resp. s' (x

n) •

I 6.2. Procedures voor niet equidistante punten

De procedures voor het geval dat de steunpunten x. niet equidistant ~

zijn, zijn analoog aan die uit §6.1. In de parameterlijsten komt nu in plaats van xO,xn een array x voor dat bij aanroep en na beeindiging de punten xi bevat. Voorwaarde Xo < xl < ••• < X

n' De beschrijving van de overige formele parameters vindt men bij de overeenkomstige procedure in §6.1. We geven nog de heading van elke procedure.

5.2.1. De procedure Reinschneqlambda

procedure Reinschneqlambda (n,x,y,d,lambda,a,da,d2a,term); value n,lambda; integer n,term;

!!!h

lambda;

!!!h

array x,y,d,a,da,d2a;

5.2.2. De procedure ReinschnegS

procedure ReinschneqS (n,x,y,d,S,a,da,d2a,lambda,term)i value n,Si integer n,termj

real S,lambda;

!!!h

array x,y,~,a,da,d2a; 5.2.3. De procedure CvGneq

procedure CvGneq (n,x,y,d,a,da,d2a,lambda,term); yalue ni integer n,termj

>.2.4. De procedure splineinterpolneq

procedure splineinterpolneq (n,x,y,dy,d2y,term); value n; integer n,term; real array x,y,dy,d2y;

;.2.5.

De procedure splinederivneg

real procedure splinederivneq (n,x,a,da,d2a,p,t);

Literatuur

In de volgende lijst worden naast de referenties waarnaar in de vooraf-gaande tekst wordt verwezen, nog enkele andere referenties gegeven waarin data fitting wordt behandeld.

[1] Boor, C. de: A practical guide to splines. New York etc.,

Springer-Verlag, 1978 (Applied Mathematical Sciences, vol. 27).

[2] Dierckx, P.: An algorithm for smoothing, differentiation and inte-gration of experimental data using spline f~nctions. J. Comput. Appl. Math.

1

(1975), 165-184.

[3] Greville" T.N.E.: Introduction to spline functions, in: Theory and applications of spline functions, T.N.E. Greville ed., New York etc., Academic Press, 1969, pp. 1-35.

[4] Lyche, T. and Schumaker, L.: Computation of smoothing and inter-polating natura~ splines via local bases. SIAM J. of analysis

1Q

(1973), 1027-1038.

[5 ] Lyche, T. and Schumaker, L.: Procedures for computing smoothing and interpolating natural splines. Communications of the A.C.M.

12

(1974), 463-467.

[6] Powell, M.J.D.: Curve fitting by splines in one variable, in: Nume-rical approximations to functions and data, J.G. Hayes ed., The Athlone Press, 1970, pp. 65-83.

[7] Reinsch, C.H.: Smoothing by spline functions. Num. Math. 10 (1967), 177-183.

[8] Reinsch, C.H.: Smoothing by spline functions II. Num. Math. 16 (1971), 451-454.

[9] Schoenberg, I.J.: Spline functions and the problem of graduation. P~oc. Nat. Acad. Sci. ~ (1964), 947-950.

[10] Whittaker, E.T.: On a new method of graduation. Proc. Edinburgh

Math. Soc. ~ (1923),63-75.

[11] Woodford, Ch.: An algorithm for data smoothing using spline functions • . BIT 10 (1970), 501-510.