Tillämpad matematik: Splineapproximation

Niels Möller 24 maj 1995

Sammanfattning

Jag löser problemet att anpassa en jämn kurva till punkterna (x_{k}, y_{k}) genom att minimera

∫_{a}^{b}(u′′(x))^{2}dx + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})^{2}

Det visar sig att optimalt u är en kubisk spline på [x_{1}, x_{n}], med naturliga randvillkor. Då λ → 0 närmar sig kurvan rät linje, som i minstakvadratmening approximerar punkterna. Då λ → ∞ närmar sig kurvan en interpolerande spline.

Innehåll

1 Problemställning 2 2 Existens 2 3 Variationsräkning 4 4 Vektorrumsstruktur 5 5 Skalärprodukt 6 6 Projektion 7 7 En lämplig bas 8 8 Parameterberoende 9

1 Problemställning

Problemet är att finna en reellvärd funktion u som minimerar

J(u) = ∫_{a}^{b}(u′′(x))^{2}dx + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})^{2} (1)

där (x_{k}, y_{k}) är n givna punkter, a < x_{1} < x_{2} < ⋯ < x_{n} < b, och λ > 0. Man måste förstås också specificera över vilket rum av u:n man vill minimera J. Jag betraktar U = mängden av alla u sådana att

• u är kontinuerligt deriverbar på [a, b],

• u′ är absolutkontinuerlig (vilket medför att u′′ existerar nästan överallt),

• u′′ ∈ L^{2}.

U är förstås ett vektorrum.

Den första uppgiften är att visa att J(u) verkligen har ett minimum i U. Därefter visar jag att den minimerande funktionen är en spline och att den är entydigt bestämd. Givet punkterna (x_{k}, y_{k}) kan den minimerande funktionen för varje λ bestämmas ur ett linjärt ekvationssystem, vilket används för att undersöka vad som händer i gränsfallen λ → 0 och λ → ∞.

2 Existens

Jag använder i detta avsnitt beteckningarna

(u, v) = ∫_{a}^{b} u′′v′′ + λ ∑_{k=1}^{n} u(x_{k})v(x_{k})
∥ u ∥ = \sqrt{(u, u)}

för u, v ∈ U. Definiera också en linjär avbildning σ : C[x_{1}, x_{n}] → ℝ^{n} genom

σ(u) = (u(x_{1}), …, u(x_{n})) (2)

Då är

∥ u ∥^{2}=∥ u′′ ∥_{2}^{2} +λ ∥ σ(u) ∥^{2} (3)

där normerna i högerledet är vanlig 2-norm på L^{2}[a, b] och ℝ^{n}.

Vi har, för alla u, v ∈ U, att

J(u + v) = J(u) + J(v) + 2(u, v) − λ ∑_{k=1}^{n} y_{k}^{2} (4)

Lemma 1 Det finns en funktion u^{∗} ∈ U sådan att J(v) ≥ J(u^{∗}), för alla v ∈ U.

Bevis Låt J^{∗} = inf_{v∈U} J(v). Eftersom J(v) ≥ 0 för alla v ∈ U så är J^{∗} ≥ 0. Låt {u_{i}} ⊂ U vara en följd sådan att

lim_{i→∞} J(u_{i}) = J^{∗} (5)

Jag visar att denna följd konvergerar mot en funktion u^{∗} ∈ U, och att detta u^{∗} uppfyller J(u^{∗}) = J^{∗}. Ekvation (4) ger, med u = (u_{i} + u_{j})/2 och v = ±(u_{i} − u_{j})/2, att

\frac{J(u_{i}) + J(u_{j})}{2} = J (\frac{u_{i} + u_{j}}{2} ) + ∥ u_{i} − u_{j} ∥^{2}
≥ J^{∗}+ ∥ u_{i} − u_{j} ∥^{2}

där vi i den sista olikheten har utnyttjat att J(u) ≥ J^{∗} för alla u ∈ U, och speciellt för (u_{i} + u_{j})/2. Eftersom J(u_{i}), J(u_{j}) → J^{∗} får vi

∥ u_{i} − u_{j} ∥→ 0 då i, j → ∞ (6)

Detta betyder enligt (3) att u′′_{i} är en cauchyföljd i L^{2}[a, b], och att σ(u_{i}) är en cauchyföljd i ℝ^{n}. Eftersom L^{2}[a, b] och ℝ^{n} är fullständiga finns det därför f ∈ L^{2}[a, b] och η ∈ ℝ^{n} så att

(7)  u′′_{i} → f i L^{2}
(8)  u_{i}(x_{k}) → η_{k}

Eftersom intervallet [a, b] är begränsat medför (7) att u′′_{i} → f i L^{1} så att vi kan definiera

g_{i}(x) = ∫_{x_{1}}^{x} u′′_{i}(ξ)dξ
g(x) = ∫_{x_{1}}^{x} f(ξ)dξ

Genom denna definition blir g absolutkontinuerlig, och

|g(x) − g_{i}(x)| ≤ ∫_{a}^{b} |f(ξ) − u′′_{i}(ξ)| → 0

så att g_{i} → g likformigt. Låt C_{1}, C_{2} vara två konstanter och sätt

u^{∗}(x) = C_{1} + C_{2}(x − x_{1}) + ∫_{x_{1}}^{x} g(ξ)dξ

Det är klart u^{∗} ∈ U och att ∥ u^{∗}′′ − u′′_{i} ∥_{2}→ 0, eftersom u^{∗}′′ = f. Nu ska vi välja C_{1}, C_{2} så att även ∥ σ(u^{∗} − u_{i}) ∥→ 0. För alla x ∈ [a, b] gäller

u(x) − u_{i}(x) = C_{1} − u_{i}(x_{1}) + (C_{2} − u′_{i}(x_{1}))(x − x_{1}) + ∫_{x_{1}}^{x}(g(ξ) − g_{i}(ξ))dξ

Integralen ovan går mot noll då i → ∞ eftersom g_{i} → g likformigt. Vi vet också att u_{i}(x_{k}) → η_{k}. Sätt nu in x = x_{n}. Då följer att lim u′_{i}(x_{1}) existerar. Välj C_{1} = lim_{i→∞} u_{i}(x_{1}) = η_{1} och C_{2} = lim_{i→∞} u′_{i}(x_{1}). Då gäller

u_{i}(x) → u^{∗}(x) för alla x ∈ [a, b]

och ur (8) fås

u^{∗}(x_{k}) = lim_{i→∞} u_{i}(x_{k}) = η_{k}

Nu återstår bara att visa att J(u^{∗}) = J^{∗}. Men det följer av att normerna i (3) är kontinuerliga, tillsammans med (7) och (8). Alltså är J(u^{∗}) = J^{∗} ≤ J(v), för alla v ∈ U. □

Nästa steg är att visa att denna funktion u^{∗} ∈ U i själva verket ligger i ett betydligt mindre rum.

3 Variationsräkning

Antag att u ∈ U minimerar J(u). Låt ϕ ∈ C^{∞}[a, b] vara godtyckligt, och betrakta I(α) = J(u + αϕ). Eftersom u är optimal och u + αϕ ∈ U för alla α så har I(α) ha ett minimum vid α = 0. Eftersom I(α) är ett andragradspolynom i α, så följer

0 = I′(0) = 2 ∫_{a}^{b} u′′(x)ϕ′′(x)dx + 2λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})ϕ(x_{k}) (9)

Då u′ är absolutkontinuerlig kan vi partialintegrera. Vi begränsar oss ett tag till sådana ϕ som har stöd i det öppna intervallet (a,b). Då ger (9) efter två partialintegrationer att

0 = ∫_{a}^{b} u(x)ϕ^{(4)}(x)dx + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})ϕ(x_{k})

Om u betraktas som en distribution på (a, b) betyder det att

u^{(4)} + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})δ_{x_{k}} = 0 i distributionsmening.

Integrera, och låt h beteckna funktionen

h(x) = \begin{cases}0 & om x < 0 \\ 1 & om x > 0 \end{cases}

Då är

(10) u^{(3)} + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})h(x − x_{k}) = konstant

Tydligen är u^{(3)} styckvis konstant, diskontinuerlig endast i punkterna x_{k}. Genom att integrera ännu ett par gånger ser vi att u′′ är kontinuerlig och styckvis linjär, så att u är kontinuerligt deriverbar två gånger, och på varje intervall (a, x_{1}), (x_{1}, x_{2}), …, (x_{n}, b) är u(x) ett tredjegradspolynom.

Lämna nu begränsningen att supp(ϕ) ⊂ (a, b). Låt ϕ ∈ C^{∞}[a, b] ha stöd i [a, x_{1}). Vi vet enligt ovan att u′′ är linjär på [a, x_{1}], så att u′′ = Cx + D för a ≤ x ≤ b. (9) ger

0 = ∫_{a}^{b} u′′ϕ′′ + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})ϕ(x_{k})

= ∫_{a}^{x_{1}}(Cx + D)ϕ′′(x)dx
= [(Cx + D)ϕ′(x) − Cϕ(x)]_{a}^{x_{1}}
= −(Ca + D)ϕ′(a) + Cϕ(a)

Eftersom likheten ska gälla för alla ϕ, och därmed för alla värden på ϕ(a) och ϕ′(a), så följer att C = D = 0. Alltså är u′′ = 0 för a ≤ x ≤ x_{1}. Samma resonemang ger att u′′ = 0 för x_{n} ≤ x ≤ b.

Vi har alltså följande

Lemma 2 Antag att u ∈ U, och att u minimerar J(u). Då gäller

• u ∈ C^{2}[a, b],

• u är på varje intervall [x_{k}, x_{k+1}] ett polynom av grad högst tre,

• u är linjär på intervallen [a, x_{1}] och [x_{n}, b].

• Speciellt är u′′(x_{1}) = u′′(x_{2}) = 0.

Med andra ord är u en kubisk spline på [x_{1}, x_{n}], med naturliga randvillkor.

Genom att sätta in u^{(3)}(x) = 0 för x < x_{1} och x > x_{n} i (10) ser vi också att

∑_{k=1}^{n} u(x_{k}) = ∑_{k=1}^{n} y_{k}

vilket är ganska naturligt. I själva verket ger (10) och lemma 2 tillräcklig information för att ställa upp ett linjärt ekvationssystem för tredjegradspolynomens koefficienter. Men för att besvara frågor om lösningens struktur, om dess entydighet, och om inverkan av olika val av λ, är denna “lösning” av problemet inte till stor hjälp.

4 Vektorrumsstruktur

Låt W = mängden av alla u ∈ C^{2}[a, b] sådana att

• u på varje intervall [x_{k}, x_{k+1}] är ett polynom av högst grad tre,

• u är linjär på intervallen [a, x_{1}] och [x_{n}, b].

W är ett linjärt underrum av U, och lemma 2 säger precis att varje optimalt u ∈ U i själva verket ligger i W. För att visa att den minimerande funktionen är entydigt bestämd räcker det därför visa entydighet under det extra villkoret att u ∈ W.

Ett element i u ∈ W är entydigt bestämt av de n talen

u(x_{1}), u′′(x_{2}), …, u′′(x_{n−1}), u(x_{n})

I avsnitt 7 används detta för att konstruera en bas för W med n element. Det följer då att dim(W) = n.

W ⊂ C^{2}[a, b], men vi kan också betrakta W som ett underrum av C^{2}[x_{1}, x_{n}], genom att identifiera varje element i W med dess restriktion till [x_{1}, x_{n}]. För om två element i W är lika på [x_{1}, x_{n}] så är de lika på hela [a, b], eftersom de är kontinuerligt deriverbara, och linjära på [a, x_{1}] och [x_{n}, b].

Nu bäddar vi in W i ett större vektorrum, genom att ta bort kravet att derivatorna ska vara kontinuerliga. Låt V = mängden av alla u ∈ C[x_{1}, x_{n}] sådana att u på varje intervall [x_{k}, x_{k+1}] är ett polynom av grad högst tre. V är ett ändligdimensionellt vektorrum; man ser lätt att dimensionen är högst 4(n − 1). W är ett linjärt underrum av V (om vi identifierar elementen i W med sina restriktioner till [x_{1}, x_{n}]). Däremot gäller förstås varken V ⊂ U eller U ⊂ V.

5 Skalärprodukt

Nästa mål är att beskriva funktionalen J med hjälp av en skalärprodukt på ett hilbertrum. Vi behöver därför en skalärprodukt på V. Denna skalärprodukt kommer att bero av λ, eftersom J gör det, men vi kan bygga upp den av λoberoende delar.

Vi börjar med att undersöka avbildningen σ : u → (u(x_{1}), …, u(x_{n}), som definierades i avsnitt 2. σ är förstås inte injektiv på V, eftersom ℝ^{n} har lägre dimension än V. Men dim W = dim ℝ^{n}. Detta, tillsammans med följande lemma, visar att restriktionen av σ till W är bijektiv.

Lemma 3 σ är injektiv på W.

Bevis Antag att v ∈ W och att v(x_{1}) = ⋯ = v(x_{n}) = 0. Vi behöver visa att detta medför v = 0.

Eftersom v ∈ W är v′′(x_{1}) = 0. Om även v′(x_{1}) = 0 måste v(x) = c(x − x_{1})^{3} för x_{1} ≤ x ≤ x_{2}, och v(x_{2}) = 0 ger att v = 0 på [x_{1}, x_{2}]. Då är även v′(x_{2}) = v′′(x_{2}) = 0, och induktion ger v = 0 på hela [x_{1}, x_{n}].

Om vi i stället har v′(x_{1}) ≠ 0 så måste v′(x_{2}) och v′′(x_{2}) vara nollskilda och ha samma tecken, eftersom v(x_{1}) = v(x_{2}) = 0 och v är ett tredjegradspolynom på [x_{1}, x_{2}], med inflexionspunkt i x_{1}. På [x_{2}, x_{3}] uppfyller därför v

• v(x_{2}) = v(x_{3}) = 0,

• v′(x_{2}) och v′′(x_{2}) är nollskilda, och har samma tecken,

• v är ett tredjegradspolynom.

Då måste även v′(x_{3}) och v′′(x_{3}) vara nollskilda och ha samma tecken. Induktion ger att v′′(x_{k}) ≠ 0 för alla k = 2, 3, …, n.

Men v′′(x_{n}) = 0, eftersom v ∈ W. Denna motsägelse visar att v′(x_{1}) ≠ 0 är omöjligt, och det följer att v′(x_{1}) = 0, och därmed v = 0, vilket avslutar beviset. □

Lemmat medför att det finns ett entydigt bestämt s ∈ W sådant att σ(s) = (y_{1}, …, y_{n}). Detta s är den kubiska spline med naturliga randvillkor som interpolerar punkterna (x_{k}, y_{k}). Vi kommer att behöva denna funktion igen i avsnitt 8.

Nu kan vi definiera ett par ickenegativa bilinjära former på V. Låt u, v ∈ V och sätt

u · v = σ(u) · σ(v) = ∑_{k=1}^{n} u(x_{k})v(x_{k})

(u, v)_{0} = ∫_{x_{1}}^{x_{n}} u′′(x)v′′(x)dx = ∑_{k=1}^{n−1} ∫_{x_{k}}^{x_{k+1}} u′′(x)v′′(x)dx

(u, v)_{λ} = (u, v)_{0} + λu · v
∥ u ∥_{λ}= \sqrt{(u, u)_{λ}}

Observera att det inte är samma norm och skalärprodukt som i avsnitt 2.

Lemma 3 ger att u · v är positivt definit på W. Dessutom är (u, v)_{λ} positivt definit på hela V. För anta att (u, u)_{λ} = 0. Då måste u · u = (u, u)_{0} = 0. Av (u, u)_{0} = 0 följer att u′′ = 0 och u är linjär på varje intervall [x_{k}, x_{k+1}]. Då medför u · u = 0 att u = 0 på [x_{1}, x_{n}].

Eftersom V har ändlig dimension så är V fullständigt med avseende på varje skalärprodukt. Alltså är V ett hilbertrum, med skalärprodukten (u, v)_{λ}. Av ändligdimensionaliteten följer också att W är ett slutet underrum i V.

Vi kommer flera gånger att använda följande egenskap hos (u, v)_{λ}: Om u′′ = 0 så är (u, v)_{λ} = λu · v.

Låt nu f vara den funktion som interpolerar punkterna (x_{k}, y_{k}) linjärt. Med h_{k} = x_{k+1} − x_{k} är

f(x) = y_{k} + \frac{1}{h_{k}} (x − x_{k})(y_{k+1} − y_{k}) för x_{k} ≤ x ≤ x_{k+1}

Eftersom f′′ = 0 får vi nu, för alla u ∈ W,

∥ u − f ∥_{λ}^{2} = ∫_{x_{1}}^{x_{n}}(u′′(x) − f′′(x))^{2}dx + λ ∑_{k=1}^{n}(u(x_{k}) − f(x_{k}))^{2}

= ∫_{x_{1}}^{x_{n}}(u′′(x))^{2}dx + λ ∑_{k=1}^{n}(u(x_{k}) − y_{k})^{2}

= J(u)

Denna beskrivning av J, tillsammans med lemma 1 och 2, ger

Lemma 4 Det finns ett entydigt bestämt u^{∗} ∈ U som minimerar J(u), nämligen den ortogonala projektionen av f i W med avseende på skalärprodukten (u, v)_{λ}.

6 Projektion

Vi kan dela upp W i två ortogonala underrum. Låt L = rummet av alla linjära funktioner på [x_{1}, x_{n}]. Låt, för varje u ∈ V, P (u) vara den ortogonala projektionen av u i L. P (u) är entydigt bestämd av att

(P (u), l)_{λ} = (u, l)_{λ} för alla l ∈ L

Men l′′ = 0 för alla l ∈ L, så att (v, l)_{λ} = λv · l för alla v och l, och villkoret ovan är därför ekvivalent med

P (u) · l = u · l för alla l ∈ L (11) vilket uppenbarligen inte beror av λ. Alltså beror P (u) ej av λ, vilket man kanske inte hade väntat sig.

Låt L^{⊥} vara det ortogonala komplementet av L i W. Eftersom L^{⊥} är nollrummet till P i W följer att även L^{⊥} är λ-oberoende. Med andra ord har L samma ortogonala komplement i W oavsett vilken av skalärprodukterna (u, v)_{λ} vi använder.

Projektionen av f i W kan förstås bestämmas som summan av projektionerna i L och i L^{⊥}. Beteckna dessa projektioner av f med f_{L} och f_{⊥}. Av (11) följer att f_{L} = P (f) är den räta linje som i minstakvadratmening bäst approximerar punkterna (x_{k}, y_{k}). f_{L} kan därför bestämmas med den “vanliga” minstakvadratmetoden. Det återstår att bestämma f_{⊥}. För att göra det behöver vi en bas.

7 En lämplig bas

Låt ϕ_{k} vara en funktion vars andraderivata har ett “tält” vid punkten x_{k+1},

ϕ′′_{k}(x) = \begin{cases}\frac{x−x_{k}}{h_{k}} & för x_{k} ≤ x ≤ x_{k+1} \\ \frac{x_{k+2}−x}{h_{k+1}} & för x_{k+1} ≤ x ≤ x_{k+2} \\ 0 & för övrigt \end{cases}

för k = 1, 2, … n − 2. ϕ_{k}:na är uppenbarligen linjärt oberoende, och givet en funktion i w ∈ W så finns en linjärkombination av {ϕ_{k}} som har samma andraderivata. Om vi kompletterar med de två linjära funktionerna l_{1}(x) = 1 och l_{2}(x) = x får vi därför en bas för W, vilket också bevisar att dim W = n, som utlovat.

Men ϕ_{k} är inte entydigt bestämd av att andraderivatan är ett tält; man kan lägga till godtycklig linjär funktion och fortfarande ha samma andraderivata. Till exempel kan man byta ut ϕ_{k} mot ϕ_{k} − P (ϕ_{k}), så vi kan förutsätta att ϕ_{k} ∈ L^{⊥}. Då är {ϕ_{k}}_{1}^{n−2} en bas för L^{⊥}.

Nu kan f_{⊥} bestämmas ur

(f_{⊥}, ϕ_{k})_{λ} = (f, ϕ_{k})_{λ} = λf · ϕ_{k} för k = 1, 2, … n − 2 (12)

Bilda två symmetriska matriser A och B genom

A_{ij} = (ϕ_{i}, ϕ_{j})_{0}
B_{ij} = ϕ_{i} · ϕ_{j}

så att A + λB är matrisrepresentationen av den bilinjära formen (u, v)_{λ} på L^{⊥}.

Ansätt f_{⊥} = α_{1}ϕ_{1} + ⋯ + α_{n−2}ϕ_{n−2}, och låt β vara vektorn β_{i} = f · ϕ_{i}. Då är (12) ekvivalent med det linjära ekvationssystemet

(A + λB)α = λβ (13)

A + λB och B är och positivt definita, eftersom (u, v)_{λ} och u · v är positivt definita på W, och därmed även på L^{⊥}. Även A är positivt definit, vilket man kan se genom at räkna ut matriselementen. A är tridiagonal, och de nollskilda elementen är

\begin{cases}A_{k,k} & = \frac{1}{3}(h_{k} + h_{k+1}) \\ A_{k,k+1} = A_{k+1,k} & = \frac{1}{6}h_{k+1} \end{cases}

A är samma matris som brukar användas för splineinterpolation. Den är diagonaldominant, och därmed positivt definit.

Alla inblandade matriser är inverterbara, eftersom de är positivt definita. Vi kan räkna ut

α = (A + λB)^{−1}λβ = (\frac{1}{λ}A + B)^{−1}β (14)

och den sökta minimerande funktionen u^{∗} blir

u^{∗}(x) = f_{L}(x) + ∑_{k=1}^{n−2} α_{k}ϕ_{k}(x)

8 Parameterberoende

Olika värden på λ ger olika approximerande kurvor. Stora λ straffar kurvor som avviker från funktionsvärdena y_{k}, medan små λ straffar kurvor med stor krökning. Vad händer i gränsfallen λ → 0 och λ → ∞?

För små λ ger (14) att

α = λA^{−1}β + O (λ^{2})

Det följer att u^{∗} → f_{L} likformigt då λ → 0.

Då λ → ∞ gäller, också enligt (14), att

α = B^{−1}β + \frac{1}{λ}B^{−1}AB^{−1}β + O (\frac{1}{λ^{2}} ) (15)

Vad är då B^{−1}β? Låt s vara den interpolerande splinefunktionen, vars existens och entydighet följer av lemma 3. Då är s det entydigt bestämda element i W som minimerar

∥ σ(u − f) ∥^{2}= ∑_{k=1}^{n}(u(x_{k}) − f(x_{k}))^{2}

av den enkla anledningen att σ(s − f) = 0. Det följer att

(s − f) · v = 0 för alla v ∈ W (16)

Genom att först ta v ∈ L och sedan v = ϕ_{k} ser vi att

P (s) = P (f) = f_{L}

s · ϕ_{k} = f · ϕ_{k} = β_{k} för k = 1, 2, …, n − 2 (17)

Vi kan utveckla s − P (s) i basen ϕ_{k}, eftersom s − P (s) ∈ L^{⊥}. Vi får

s(x) = f_{L}(x) + ∑_{k=1}^{n−2} γ_{k}ϕ_{k}(x)

Nu ger (17) och definitionen av B att Bγ = β. Men B är inverterbar, så att

B^{−1}β = γ

Tillsammans med (15) fås nu

\begin{aligned}u^{∗} & = f_{L} + ∑_{k=1}^{n−2}(γ_{k} + \frac{1}{λ}(B^{−1}Aγ)_{k})ϕ_{k} + O (\frac{1}{λ^{2}} ) \\ & = s + \frac{1}{λ} ∑_{k=1}^{n−2}(B^{−1}Aγ)_{k}ϕ_{k} + O (\frac{1}{λ^{2}} ) \end{aligned} (18)

så att u^{∗} → s likformigt då λ → ∞. Vi har alltså följande

Resultat Låt u_{λ} vara den funktion som minimerar J(u), för varje λ > 0. Låt l vara den räta linje som approximerar punkterna (x_{k}, y_{k}) i minstakvadratmening, och låt s vara den kubiska spline, med naturliga randvillkor, som interpolerar samma punkter. Då gäller

• u_{λ}(x) → l(x) likformigt då λ → 0

• u_{λ}(x) → s(x) likformigt då λ → ∞