Algorithmen voor het vereffenen en aanpassen van experimentele data

van Ginneken, C. J. J. M. (1983). Algorithmen voor het vereffenen en aanpassen van experimentele data. (Computing centre note; Vol. 19). Technische Hogeschool Eindhoven.

Document status and date: Gepubliceerd: 01/01/1983 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.

--

II. . .

_Bibliotheek

_{TechnischeHogpc,.!,,,,\,, "}

·...,.,pqlde datum

maart 1983

Eindhoven University of Technology Computing Centre Note 19

ALGORITHMEN VOOR HET VEREFFENEN EN AANPASSEN VAN EXPERIMENTELE DATA

Onderafdeling der Wiskunde

Memorandum 1980-17 september 1980

Algorithmen voor het vereffenen en aanpassen van experimentele data door

C.J.J.M van Ginneken

Technische Hogeschool Onderafdeling der Wiskunde PO Box 513, Eindhoven Nederland

Inhoudsopgave

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 bIz. 4 4 5 7 14 25 33

i

o

n

Y_i een, bijvoorbeeld door meting verkregen, benadering is voor een onbekende functiewaarde f(xi)' Gevraagd wordt om een functie 9 te bepalen die de on-nauwkeurigheden in y. vereffent.

~

In het geval dat bijvoorbe~ldbekend is dat de gemeten funct1e f behoort tot een door enkele parameters gekarakteriseerde klasse is het aan te bevelen 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 gebru1k geschikte waarden moeten worden gekozen, wat soms moeilijk is. In paragraaf 4 wordt

een methode beschreven die sterk verwant is aan die uit 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 volgende zijn (xi,Y

i), i

=

O,l, .••,n, n ~ 1 gegeven data met X

o < xl < •••< Xn• De getallen di , i

=

O,l, ••• ,n zijn gegeven zogenaamde gewichtsfactoren.

Definities Ck(xO'X

n] is de verzameling van de functies, die k-maal differentieerbaar zijn op [x IX ] en waarvan de k-de afgeleide continu is;

o

n

ak[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 (xi,xi+_{l ), i}

=

O,l, ••• ,n-l. en kwadratisch integreerbaar op [XO,x_n];

de functionalen X_n E (f) :=

f

[f(m) (x)]2dx , f € H(m)(X x] m _X 0' n

o

n ₂ K(f) :=

r

d. x'(f(x_i )

-

_{Yi) ;} 1=0 1.

-

5-S2m_l is de verzameling van de natural spline functies van de graad 2m..- 1 met knooppunten x_O,x₁" " , x_n'

Pro memorie geven we de definitie van een natural spline.

Een natural spline functie van de graad 2m - 1 met knooppunten xO,x1"",x n is een op [xO,xnJ gedefinieerde functie s met de volgende eigenschappen

(i) (ii) (iii)

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

2m-2 ~ 1+

S € C [xo,x ]

S(k) (xO) = s(kV(X

n) = 0, k = m,m+ , ••• , m- •1 2 2 3. Methoden van Schoenberg en Reinsch

Om een functie te bepalen die de onnauwkeurigheden in y. verdisconteert

~

stelde Schoenberg (1964), uitgaande van een idee van Whittaker (1923), het volgende voor.

Zij m S n+1 en A > 0 gegeven. Bepaal gx die de functionaal

(3.1)

minimaliseert.

De parameter A is een maat voor het compromis tussen de "gl adheid" van gA' die wordt gemeten met Em(gX) en de afwijking van gA van de data Yi' die wordt gemeten met de som van kwadraten K(gX)'

Schoenberg toonde het volgende aan. Er bestaat een eenduidig,bepaalde gA en deze behoort tot de klasse S2m_1' Voor een bewijs zie Greville (1969). Als X t 0 dan gaat gx over in het polynoom go met graad S m - 1 waarvoor K minimaal is en als A+ ~ nadert gA tot de natural spline die de data Y

i interpoleert. Er geldt dat K(gA) een monotoon niet stijgende en E(gX) een monotoon niet dalende functie van A is.

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

(1964) en Reinsch (1967) ook nog de volgende aanpak voor.

Zij (in plaats van A) een getal S 0 < S <

~

gegeven. Bepaal h

_s

die ) de functionaal Em(g), g E Hm[Xo'X

n] minimaliseert onder de nevenvoor- . (3.2) waarde K(g) $ S.

In het geval dat bijvoorbeeld bekend is dat de data y.

~

afhankelijke ruiscomponenten E. bevatten met variantie

~

de hand de gewichtsfactoren d. gelijk te nemen aan

cr~2

~ ~

lijk te nemen aan n + 1.

onderling

onaf-cr~,

ligt het voor

~

en S ongeveer

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

oplossing

:= lim K(gA) en zij 0 < S < K(gO)' dan heeft (3.2) precies een

HO

*

h

s

.

Zij A uit (3.1) zodanig dat K(g*,)_). = S, dan geldt h

s

= g .•_).*

Dit is eenvoudig in te zien. Immers K(g).) is een niet stijgende functie

*

van A met 0

=

K(g~) < S < K(gO)' dus er is een A met K(g )

=

S. Voor

).* m iedere g E H [XO,x n] met K(g) $ S geldt E ( .. ) A* ( ) mg

*

+ K g

*

A A * * * = Em(g ~) +A S $ E (g) + ). K(g) $ E (g) + A

s,

A m m

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

m A* m A*

OOk geldt voor iedere oplossing h van (3.2)

S

*

*

*

E (h ) + A K(h

s

)

$ E (g ) + A S = Em(g *) + A K(g *).

m S m A* A A

omdat de oplossing van (3.1) eenduidig is volgt hieruit dat h

_s

=

g

*..

A

Als S

=

0, dan is h

S de natural spline die de data Yi interpoleert en als S ~ K(gO)' dan is bijvoorbeeld h

_s

= go een oplossing.

-

7-Algorithmen voor het berekenen van h

_s

zijn gegeven door Reinsch (1967) voor het geval m

=

2, Woodford (1970), Lyche and Schumaker (1973,1974)

•

voor willekeurige m ~ 1. De bepaling van h

_s

gebeurt steeds, indien 0<5 < K(gO)' door het oplossen van de van A afhankelijke vergelijking K(gA)

=

5. De bepaling van gA' 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 5 uit (3.2) aan te geven, bIijven 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 weI 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) € 53 kan

op het interval [x

i,xi+1] geschreven worden als

(4.1 )

waarin hi+~ := _{xi +1 - xi·}

We stellen a

i := g(xi), ci := qlt(xi), i = O,l, ••• ,n. Uit de continuiteit van de eerste afgeleide in de punten Xj' j = 1,2, ••• ,n-l volgen de vergelij-kingen

=

c_n

=

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

=

QTa,

-

8-T T

waarin c := (c

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

(n + I) x (n - 1) matrix 1 T=-6 Q

=

2 (h1

/.2

+ _{h 3/ 2 )} h 3/ 2

-o

o

-

_-

_-

-...

o

-

- - _ hn-3/ 2 h_{n -3/ 2} 2(h_{n _3/ 2} + h_{n -1/ 2)}

o

... ... -1 ...

_....

hn-3/ 2 -1 -1 -h_{n -3/ 2 -.hn-1/ 2} ... _ -1 hn-1/ 2 g"(x)

=

x. 1 - x + J.+ C. 1 h J.+ i+1.2

waaruit (na enig rekenwerk) volgt dat (4.3) of met (4.2) T

=

c 're, (4.4) E T -1 T 2(g)

=

a QT Q a •

Zij D de diagonaalmatrix met D.. := d., i

=

O,l, •••,n, dan geldt

.1..1. .1.

(4.5) E T -1 T T

-

9--'1 T

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

(4.6)

Uit (4.2) en (4.6) voIgt het stelsel (4.7) (QT -10 Q + AT)c = AQTy

en de relatie

(4.8) a

=

y - - 01 -1QC

A

Door het stelsel (4.7) op te lossen, waarvan de matrix positief definiet is met bandbreedte 5, wordt c

i

=

g~(xi) verkregen en vervolgens uit (4.8)

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 gA(~), 0 < X <

'7J,

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

~

kiezen we op grond van de volgende overwegingen. Als gX(x) geschikt is, dan moet de vector

een ru1ssignaal voorstellen. Zij 0X(x) de natural cubic spline die r_X inter-poleert, dus

<\

(xi)

=

(r

A)i ' Voor de "energie" E2(OX) voIgt uit (4.4)

(4.9)

Indien r

Aeen ruissignaal is met onderling ongecorreleerde componenten, waar-van de verwachtingswaarde nul is voor iedere component, dan geldt voor de

~ 10 ...

(4.10)

waarin 0 de diagonaalmatrix is bestaande uit de diagonaalelementen van B. Uit (4.10) destilleren we nu het volgende criterium. Probeer ). zo te kiezen dat

(4.11) F()')

Het is (mij) niet gelUkt om, eventueel onder bepaaide 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().) in de praktijk.

Veronderstel eenvoudigheidshalve dat de punten x. equidistant zijn, dus

~

....

x_i

+

_{1 - x.}

=

h. De elementen D'i zijn positief en in de prde van grootte

-3 1. J. ~

van h • Op grond hiervan zal de orde van grootte van r).Dr). gelijk zijn

T 3

aan (r).r).)/h , notatie:

(4.12)

(4.13)

Er geldt

(Y_i is een benaderlng voorf{x

i

».

Als ). klein Is, dan is g). bijna lineair en dan zal in het algemeen (mits f niet bijna lineair is) de term f{x_i) - Y_i verwaarloosbaar zijn t.o.v. g).(x

1) - f(xi ). De spline o).(x) interpoleert

dus ongeveer de functle g).(x) - f(x), welke een gladde functie is. Op grond hlervan verwachten we dat, als ). klein is, de orde van grootte van E₂(O).)

gelijk is aan

~n

[f"(x)]2dx , dus

o

T

rXn

2

E2

«\)

= r).Br)...,

1x

[fll (x)] dx.

- 11

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

r~BrA

groot wordt dan gaat r

Aen dus F(A)

Kleine waarden van A en h en als f T....

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

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

- indien F(O) ~ 0 kiezen we gX(x)

indien F(O) < 0 kiezen we gA(x), waarin A

1 de kleinste wortel is van (4.11); (AI 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

y :=

en

dan geldt dat

T T T T

rXBr

_A

=

y By + a Ba - 2y Ba

. -1 T

en hieruit voIgt met behulp van (4.2) en de definitie v~~ B

=

QT Q (4.14) rXBrT

_X

=

y ByT + c Tc - 2cT TQTy.

T T

Bet getal y By en de vector Q y zijn onafhankelijk van X.

T -1

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

T

voorkomen. De berekening van rABr

_X

komt dus in feite neer van (4.7). De berekening van

De vector c

vol-T

Q y als constanten op het oplossen

12

-gaat bij gegeven c en D eenvoudig via (4.8). De elementen D

kk kunnen als voIgt 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)

dus

D

kk

. T T

waarin zk kan worden berekend U1t LZ

k

=

Q eke Op analoge wijze wordt y By uit (4.14) berekend. In het geval dat de punten xi equidistant zijn kunnen we de volgende formules voor D

kk afleiden

13

n n-l D_OO = D

₌

(1 + a - a - a ), nn h3(l

-

an) (4.15) D 1

3613

_{(l +} n k n-k _k l,2, ••• ,n-l, kk

=

- (-48 + a. a - a

»,

=

h3 1

-

an waarin h := (x_n - x )/n

_a

en a = 7 - 4/3. -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 1}

0

"

1 -2,

,,

1, 4 1

,

,

,

_,

_,

1 1,

"

"

h

_"

_"

_"

"

"

Q

= -

_{" "}

T

=

_"

_"

_"

6

_"

,

_"

6

...

_"

_"

"

_,

_"

_"

,

_"1

,

_"

_"

_"

,

,

"

,

"

,

"

"

_"

0

"

...

"

0

"

,

~2

,

,

"

,

1

"

"

,

'1 '1 4

i3

-Zij nu P := : T, dan geldt

1

Q=-h P - 6T , waarin ~1 resp. en_1

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

(4.16) B 6

=

h3 eTp-1e I

X

I

~

1 11 I "... _______ L L _ I -1 I

X

lP - 121 + 36P

I

X

_______ L L _ I I T -1

><

I

X

Ie P e

I

I

n-l n-l

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

Zij e k de k-de eenheidsvector en z(k) := pz(k) = e • k dan geldt (k) (k)

(zl , ••• ,zn_1) .de oplossing van

(4.17) (P-1) (k)

kk

=

zk • Definieer zcik)

weglaten

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

n

(4.18.2) _{zk_l + 4zk + zk+l}

=

1, (4.18.3) zi_l + 4z_i + zi+l = 0

1 s i S k-l,

- 14

-Uit (4.18.1) en

Zo =

0 volgt (4.19) z.

=

a(Ai - A )-i

1.

o

s; i s; k,

waarin A

=

-2

+

13

en a een constante is. Evenzo voIgt uit (4.18.3) en z

=

0

n

(4.20) k s; i s; n,

waarin

e

een constante is.

Uit (4.19) en (4.20) voor i

=

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

B

kunnen worden opgelost. Na enig rekenwerk vinden we

-1 (P )kk = zk = waarin a

=

).2

=

7 - 4/3. n n-k k (1 + a - a -a), 6(1 - an)

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

5. Numer1eke experimenten

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

Voorbeeid 1:

Gegeven zijn de data Xi

=

i x w/23, i

=

0,1, ••• ,23 en Y_i de op 1 decimaal afgeronde waarde van sin(x

i). We nemen de gewichten di

=

1 en omdat we de ruiscomponenten in Y

i opvatten als random trekkingen uit een homogene

ver-de ling 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 door )C.

_15

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-nauwkeurigheid, identieke resultaten. Opvallend is het minder gladde karakter van h" in figuur 3. Dit is mijns inziens inherent aan de methode

S

omdat uit formule (4.8) blijkt dat de correctie op y, welke de ruis

1 -1

in y zou moeten representeren, gelijk is aan

I

D QC, met c

i

=

h~(xi)• . 5

n

1112 FIG.

o l l E - - - -

. l - --L. ~

o

·5

o

-.5 -1

o

fin FlO. 2

n

o

-.5 -1

o

voorbeeld :2. 11/2 flO. 3

n

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_i = 1 en S.= 2 bij de methode van Reinsch. De bij deze gegevens berekende natural spline h

17

-aangegeven met een x. De aanpassing berekend met de methode uit §4 levert

•

identieke resultaten.

x

x x

x

x x

x

x x

.

5lr--~X:'--,-

__

---~X:""_--~X;----;X;---x~--;;x-l

x

o

o

x

x

x

x

x

10 fIG. 4

x

20 Voorbeeld 3:

Gegeven zijn de .data X

o

=

0, xi

=

(i -

t

+ ri) x

~/23,

i

=

1,2, ••• ,22, x

23

=

wen y,1.

=

sin(x.)1. + 0.1 x s., waarin r. en s1.' random trekkingen zijn1. 1. uit een homogene verdeling op ~-0.5, 0.5). We nemen d

i = 1 en S

=

24 x 0.01/12. In figuur 5 is de grafiek getekend van de berekende spline h

_s

en bovendien zijn de data (xi'Yi) gemarkeerd door een X. In figuur 6 resp. 7 is h

S

resp.

h; getekend. De figuren 8,9 en 10 zijn analoog aan 5,6 en 7 maar nu voor de spline berekend met de methode uit §4.

x

x

.5 O---..L---~_X

o

R/2 R FIG. I) .5

o

-.5 -1

o

R/2 fIG. 6 R

Or-=:::::::::=---0>1

-.5 -1

o

R/2 FIG. 1

n

.5

o

:L---..I---~X

o

R/2

n

FIG. 8

o

-1

o

o

-.5 -1

o

Voorbeeld 4: 11/2 fIG. 9 1112 FIG.IO 11 11

De punten xi zijn als in voorbeeld 3 en Y

i is de op 1 decimaal afqeronde waarde van sin(x.). We nemen d. en S weer als in voorbeeld 3.

J. J.

In fiquur 11 is de qrafiek getekend van de berekende spline h

_s

en de data (xi'Y

i ) zijn gemarkeerd door een x en in fiquur 12 is de spline weerqeqeven,

berekend met de methode van §4. Opvallend is het minder qladde karakter van de spline uit fiquur 12. Een moqelijke verklarinq 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.

- 2J-5

n

nl2 FIG.II

O'E---L---.:j(

o

5

n

lV2 FIG.12 O'f---l.---...:~

o

Voorbeeld 5:

Gegeven zijn de data xi

=

ix ~/n, i

=

O,l, •••,n en Y

i de op d decimalen_-2d afgeronde waarde van sin(x

i). We nemen di

=

1 en S

=

(n + 1) x 10 /12.

Zij

-s de natural spline die de functie sin (x) interpoleert in de punten xi 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 •

- 22

-Zij r(p) de als voIgt gedefinieerde maat voor de afwijking van de p-de

m

afgeleide van s_m en s: r(p):= 1

m n + 1

I

r

1=0

In de volgende drie tabellen zijn voor n

=

23, 45, 90 en d

=

1,2, ••• ,6 de afwijkingen r(p)getabelleerd.

m

trabel 1 n =: 23

Aantal Orde van Interpolatie Reinsch t4ethode §4

~ecimalen d de afgeleide p _{r(p) r(p)} r(p) 1 2 3 1 0 2.5 10-2 2.310-2 2.310-2 1 2.5₁₀-1 6.1 10-2 6.210-2 2 8.1 9.4₁₀-2

.

9.51O~2 2 0 1.9 10-3 2.710-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.\0-4 2.2 10-4 4.010-4 1 2.0 10-3 1.710-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 10-6 3.910-6 2.910-6 1 2.9 10-5 3.210-5 2.910-5 2 6.5₁₀-4 6.3 10-4 6.510-4 6 0 1.6₁₀-7 3.2 10-7 1.610-7 1 1.6₁₀-6 2.1 10-6 1.410-6 2 4.7 10-5 4.910-5 4.6 10-5

- 23-Tabel 2 n= 45

Aantal Orde van Interpolatie Reinsch Methode §4 decimalen d de afgeleide 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.910-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₁₀-2 3.4 10-3 3.410-3 5 0 2.9 10-6 3.710-6 4.310-6 1 5.7₁₀-5 5.7 10-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~610-6 9.0} 10-6 5.810-6 2 3.2 10-4 3.210-4 3.310-4

Tabel 3 n

=

90

Aant,al Orde van Interpolatie Reinsch Methode 54 •

declmalen d de afgeleide p _rep) rep) rep)

1 2 3 1 0 2. \0-2 8.4₁₀-3 7.7₁₀-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₁₀-3 9.0 10-4 1.910-3 1 8.0 10-2 5.310-2 9.610-3 2 _{.1.7 101} 2.3₁₀-2 3.0 10-2 3 0 -4 2.1₁₀-4 1.7 10-4 2.5 10 .1 8.9 10-3 1.610-3 1.310-3 2 1.2 9.9₁₀-3 1.0 10-2 4 0 3.0₁₀-5 2.8₁₀-5 3.9 10-5 1 8.5₁₀-4 3.6 10-4 5.110-4 2 1.4₁₀-1 4.0 10-3 4.410-3 5 0 2.7₁₀-6 1.7₁₀-6 1.6₁₀-6 1 _{9.4 10-:5} 5.9₁₀-5 8.4₁₀-5 2 1.3 10-2 1.610-3 -3 1.5₁₀ 6 0 3.0 10-7 3.510-7 2.110-7 1 1.1 10-5 1.410-5 1.410-5 2 1.2₁₀-3 6.7 10-4 6.610-4

Utt de gegeven numerieke resultaten moge blijken dat de kwaliteit van de methode uit §4 vergelijkbaar is met de methode van Reinsch. Ook vele andere numerleke experimenten hebben dit bevestigd.

- 25 _

6. Procedures

In dezeparagraaf 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

J.

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

J.

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

- gl uit (3.1) (zie §6.1.1 of §6.2.1) - gs uit (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. interpoleert (zie §6.1.4 of §6.2.4).

J.

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 §6. 1.5 of §6 •.2.5) •

6.1. Procedures voor equidistante punten 6.1.1. De erocedure Reinscheqlambda

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

~ xO,xn,lambda; ~ arra~ y,d,a,da,d2a; Formele earameters n xO xn y lambda d

Bevat het rangnummer van het laatste datapunt. Voorwaarde n ~ 2. Bevat x_O'

Bevat x • Voorwaarde: x > x

O'

n n

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

Bevat de,waarde voor de gladstrijkparameter A (zie toelichting). Voorwaarde: A ~ O.

Array met grenzen [0 : n] dat bij aanroep en na beiindiging de gewichtsfactoren d_i bevat. Voorwaarde d_i > O.

-

26-a, d26-a, d2a Array's met 9renzen [0 : nJ. Indien na be~indi9in9term

=

0 dan bevatten a[i], da[i], d2a[i] de waarden gA(X_i), gl(x_i), 9

_A

(x_i) (zie toelichting).

term Bevat na be~indigingvan de procedure informatie over de af-loop van het proces.

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.

Toelichting

Zij gegeven de data (xi,y i), tieve gewichtsfactoren d

i en de Xvoldoend gladde) functie

i

=

O,l, ••• ,n met equidistante x. en

posi-l.

zij gegeven een reeel getal

A

~ O. Zij gl die de functionaal

[gfl(x)]2dx + A

[n

X

_o

m1nimaliseert. Dan geldt 9

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

i), gi(xi) en g~(xi)' i

=

O,l, ••• ,n, (zie §4).

6.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;

s

21

-/

Formele parameters

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

Bevat de waarde voor de gladstrijkparameter S (zie toelichting). Voorwaarde: S > O.

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

= a

of term

=

1 dan bevatten a[i], da[i], d2a[i] de waarden hS(x_i), hs{x_i), hs(X

i) (zie toelichting).

lambda Bevat na beeindiging indien term

=

a

de bij de waarde van S

*

behorende waarde

A

(zie toelicnting) en indien term

=

1 de term

waarde O.

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

Indien

term

=

0: de procedure heeft h

s

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, eentheoretisch onmogelijk resultaat ontdekt term

=

4: de input voldoet niet aan de genOemde restricties

n ~ 2 etc • • Toelichting

Zij gegeven de data (xi'Yi) i

=

O,l, ••• ,n met equidistante Xi en positieve gewichtsfactoren d

i en zij gegeven een reeel getal S > O. Zij h

s

een

(voldoend gladde) functie die de functionaal

_ 28 _ Uit 13 volgt

- h

S is een natural cubic spline

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

_s

eenduidlg en bestaat er precies ~~n ~* waarvoor de functie 9

*,

zoals de

proce-~

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

h

_s

:c go een oplossing.

de waarden hs(x_i), h~(xi)' _{hs(xi ) en}

onderling onafhankelijke ruiscomponen-aanhevolen d. gelijk te nemen

1.

De procedure berekent, zo mogelijk,

*'

A • Ind~en bekend is dat de data Y_i

i 2

ten £i. bevatten met variantie O"i dan wordt aan

0"1

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)i value n,xO,xn; integer n,term;

real xO,xn,lambda; real array y,d,a,da,d2aj Formele parameters

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

a,da,d2aArray's met grenzen [0 : n]. Indien na be~indigingterm

=

0,1 of 2 dan bevatten a[i], da[i], d2a[i] de waarden g, (Xi)' gi (Xi)'

,

At

1

gl (Xi) (zie toelichting).

1

lambda Bevat na be~indiging indien term

=

0,1 of 2 de waarde ~1 (zie toelichting) •

term. Bevat na beiindiging van de procedure informatie over de afloop van het proces.

Indien

toelichting)

term

=

0,1 of 2: de procedure heeft g~

1

bepaald (zie verder term ... 3

term .,. 4

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

de input voldoet niet aan de geeiste restricties n Oi!: 2, etc ••

_ 2!1 _

Toelichting

Zijqeqevende data (xi,y

i ), i;:;O,l, ••• ,nmetequidistantexi en positieve qewichtsfactoren die De procedure berekent, zo mogelijk, (zoals beschre-ven in §4) een waarde A

1 en de daarbij behorende aanpassinq

qA

zoals

1

de procedure Reinscheqlambda die aflevert. Indien na be~indiginqqeldt dat

- term

*

0, dan is 1.

1 de kleinste wortel van (4.11)

- term ;:; 1, dan is 1.₁

;:; a

en de aanpassinq gA lineair 1

- term ;:; 2, dan is het niet qelukt de kleinste wortel van (4.11) te berekenen en is 1.

1 zodaniq dat in het interval [0,1.1] qeen wortel ligt.

6.1.4 De procedure splineinterpoleg

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

real xO,xn; ~ array y,dy,d2y Formele parameters n xO, __ xn y dy,d2y term

__ Bevat het ranqnwnmer van het laatste datapunt. Voorwaarde: n ~ 2. Bevat x

O• Bevat x

n• Voorwaarde: xn > xO•

Array met qrenzen [0 : n] dat bij aanroep en na be~indiqinqde waarden Yi bevat.

Array's met qrenzen [0 : n]. Indien na be~indiqinqterm ;:;

a

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

~ J..

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

knoop-~

punten xi.

Bevat na be~indiqinqvan de procedure informatie over de afloop van het proces.

Indien

term

=

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

_ 30 _

term = 3: gedurende de berekeningen is, ten gevolge vaq afrond-fouten,een theoretisch onmogelijk resultaat ontdekt term

=

4: de input voldoet niet aan de genoemde restricties

n e:: 2, etc••

Toelichting

Zij gegeven de data (x.,y.), i = O,l, ••• ,n met equidistante x{. Zij s

1. J. ...

de natural cubic spline die de data Yi interpoleert, ofwe1 s minimali-seert de functionaal

onder de voorwaarde g(x.)

=

y .• De procedure berekent, zo mogelijk, de

J. J.

waarden s(x.), s'(x

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

J. J.

6.1.5. De erocedure splinederiveg

~eroceduresplinederiveq (n,xO,xn,a,da,d2a,p,t)i

value n,xO,xn,p,t; integer n,pi real t,xO,xn; real arrax a,da,d2a; Formele earameters

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

Bevat x O•

Bevat xn. Voorwaarde: xn > xO.

n

xO xn

a,da,d2a Array's met grenzen [0

~

n]. Bij aanroep en na afloop bevatten a[i], da[i], d2a[i] de waarden s(x_i), s'{x_i), S"{X_i } van de natural cubic spline s die de data a[i] interpoleert in de knooppunten xi.

P Orde van de te berekenen afgeleide. t Waarde van het argument.

- 31

-Toel ichting

Zij seen natural cubic spline met equidistante knooppunten xi'

i

=

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

~

sl(x

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

Op de intervallen (-~,XO) en (xn'~) is s lineair met afgeleide Sl(X_O) resp. Sl(X

n).

6.2. Procedures voor niet equidistante punten

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

~

zijn, zijn ana100g 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 X

_o

< Xl < ••• < x_n' De beschrijving van de overige forme1e parameters vindt men bij de overeenkomstige procedure in 16.1. We geven nog de heading van elke procedure.

6.2.1. De procedure Reinschneq1ambda

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

real lambdai real array x,y,d,a,da,d2ai 6.2.2. De procedure ReinschneqS

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

real S,lambdai real array x,y,d,a,da,d2ai 6.2.3. De procedure CvGneg

procedure CvGneq (n,x,y,d,a,da,~2a,lambda,term);

value n; integer n,term;

_ 32_

6.2.4. De procedure splineinterpolneq

procedure splineinterpolneq (n,x,y,dy,d2y,term); value n; integer n,term; real array x,y,dy,d2Yi 6.2.5. De procedure s2linederivneq

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

33

-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.

[1J 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 •

..!.

(1975), 165-184.

[3J 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.

[4J Lyche, T. and Schumaker, L.: computation of smoothing and inter-polating natural splines via local bases. SIAM J. of analysis 10

(1973), 1027-1038.

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

..!.l

(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.

[9J Schoenberg, I.J.: Spline functions and the problem of graduation. Proc. 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 algoritlun for data smoothing using spline functions. BIT 10 (1970), 501-510.