6.1 Interpolasjon — Lagrange og Newtons dividerte differanser

$$(x_0,f_0),\ (x_1,f_1),\ \dots,\ (x_n,f_n),$$

finn et polynom $p$ av grad høyst $n$ som treffer alle punktene:

$$p(x_k)=f_k \qquad \text{for } k=0,1,\dots,n.$$

Merk tellingen. Med tre punkter er $n=2$, og du leter etter et polynom av grad høyst 2. Fire punkter gir grad høyst 3. Antall punkter er alltid én mer enn graden — det er den vanligste kilden til rot i feilformelen senere.

Nodeverdiene er de tilhørende $f_0,\dots,f_n$. De kan være hva som helst — de kan komme fra en måling, en tabell eller en funksjon $f$ vi kjenner, men bare vil regne billig på.

Nodene trenger ikke ligge med jevn avstand. Er avstanden konstant, kaller vi nodene ekvidistante, og da forenkles flere formler.

$$\ell_k(x)=\prod_{j\ne k}\frac{x-x_j}{x_k-x_j} = \frac{(x-x_0)\cdots\widehat{(x-x_k)}\cdots(x-x_n)}{(x_k-x_0)\cdots\widehat{(x_k-x_k)}\cdots(x_k-x_n)},$$

der hatten betyr at faktoren med $j=k$ er utelatt.

Den avgjørende egenskapen er at $\ell_k$ er 1 i sin egen node og 0 i alle de andre:

$$\ell_k(x_j)=\begin{cases}1,& j=k\\ 0,& j\ne k.\end{cases}$$

Grunnen er lett å se: setter du inn $x=x_j$ med $j\ne k$, står faktoren $(x_j-x_j)$ i telleren og gjør hele produktet null. Setter du inn $x=x_k$, blir hver brøk lik 1.

Hver $\ell_k$ har grad nøyaktig $n$ (det er $n$ faktorer i telleren). Formelen står på det utdelte formelarket — tren oppslaget ved å øve på hvilke faktorer som skal utelates.

$$p(x)=\sum_{k=0}^{n} f_k\,\ell_k(x).$$

Hvorfor det virker: sett inn $x=x_j$. Alle leddene med $k\ne j$ forsvinner (fordi $\ell_k(x_j)=0$), og leddet med $k=j$ gir $f_j\cdot 1=f_j$. Altså er $p(x_j)=f_j$ for alle $j$, som er nøyaktig det vi ba om.

Styrken: ingen likninger å løse, svaret er ferdig med én gang.
Svakheten: kommer det et nytt punkt til, må alle kardinalfunksjonene skrives om fra bunn.

Grad høyst $n$, ikke nødvendigvis grad $n$.

Ligger alle punktene på en rett linje, er interpolasjonspolynomet gjennom tre av dem førstegradspolynomet — koeffisienten foran $x^2$ blir null. Oppgaveteksten sier derfor som regel «polynomet av minste grad», og det er det samme polynomet: entydighetssetningen under garanterer at det bare finnes ett.

Praktisk: blir den ledende koeffisienten null i regningen din, er det ikke en feil. Det betyr at punktene ligger på en enklere kurve enn du trodde.

$$f[x_i,x_{i+1}]=\frac{f_{i+1}-f_i}{x_{i+1}-x_i}.$$

Det er ikke annet enn et gjennomsnittlig stigningstall — den samme brøken du bruker for å finne stigningen til en rett linje.

Navnet «dividert differanse» kommer av at det er en differanse (i telleren) dividert med en differanse (i nevneren).

$$f[x_i,\dots,x_{i+k}]=\frac{f[x_{i+1},\dots,x_{i+k}]-f[x_i,\dots,x_{i+k-1}]}{x_{i+k}-x_i}.$$

Merk nevneren: den er avstanden mellom den ytterste noden til høyre og den ytterste til venstre i den gruppen du regner på — ikke avstanden mellom to nabonoder. Det er den klassiske feilen i tabellen.

En dividert differanse er symmetrisk i argumentene sine: rekkefølgen på nodene spiller ingen rolle for verdien. Det er en fin kontroll — regner du tabellen ovenfra og ned eller nedenfra og opp, skal du få det samme.

$$\begin{array}{c|c|c|c} x_0 & f_0 & & \\ & & f[x_0,x_1] & \\ x_1 & f_1 & & f[x_0,x_1,x_2] \\ & & f[x_1,x_2] & \\ x_2 & f_2 & & \end{array}$$

Koeffisientene i Newton-formen er den øverste skråkanten — altså $f_0$, $f[x_0,x_1]$, $f[x_0,x_1,x_2]$ og så videre.

Selve utfyllingen må kunnes. Formen på tabellen finnes på det utdelte formelarket, men rekkefølgen du fyller den ut i, og hvilke tall som havner i nevneren, må sitte i fingrene.

$$p(x)=f_0+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\dots$$
$$\dots+f[x_0,\dots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1}).$$

Styrken: hvert ledd forsvinner i alle nodene til venstre for seg. Setter du inn $x=x_0$, står bare $f_0$ igjen; setter du inn $x=x_1$, står de to første leddene igjen, og så videre. Derfor legger et nytt punkt bare til ett nytt ledd — alt du har regnet før, står.

Formen står på det utdelte formelarket — tren oppslaget. Det som må kunnes, er å fylle ut tabellen og lese av riktig skråkant.

$$\text{koeffisienten foran } x^n \ =\ f[x_0,x_1,\dots,x_n].$$

Det ser du direkte av Newton-formen: bare det siste leddet inneholder $x^n$, og faktoren foran produktet er $f[x_0,\dots,x_n]$.

To nyttige følger. (1) Er den siste dividerte differansen null, har polynomet lavere grad enn du trodde. (2) Har du regnet ut polynomet på Lagrange-form, kan du sammenlikne den ledende koeffisienten med den siste dividerte differansen som en rask kontroll — de må være like.

$$\omega(x)=(x-x_0)(x-x_1)\cdots(x-x_n)=\prod_{k=0}^{n}(x-x_k).$$

Det har grad $n+1$ og er null i hver eneste node — akkurat som feilen må være, siden $p$ treffer nøyaktig der.

Nodepolynomet er den delen av feilen du selv styrer, gjennom valget av noder. Funksjonsdelen er gitt; plasseringen av nodene er din.

$$\varepsilon_n(x)=f(x)-p(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\,\prod_{k=0}^{n}(x-x_k),$$

der $\xi$ er et ukjent punkt i det minste intervallet som inneholder både $x$ og alle nodene.

Tre ting å merke seg.

1. Du kjenner ikke $\xi$ — og du trenger det ikke. Du erstatter $f^{(n+1)}(\xi)$ med en øvre grense $\max|f^{(n+1)}|$ over intervallet, og får et feilanslag.
2. Ordenen på den deriverte er $n+1$, én mer enn graden. Med tre noder ($n=2$) er det den tredjederiverte, og fakultetet er $3!=6$.
3. Feilen er null i nodene, som den skal være.

Formelen står på det utdelte formelarket — tren oppslaget. Det du må gjøre selv, er å telle riktig: hvor mange noder, hvilken $n$, hvilken derivert, hvilket fakultet.

$$|\varepsilon_n(x)|\ \le\ \frac{M_{n+1}}{(n+1)!}\,\max_{x\in[a,b]}\left|\omega(x)\right|,\qquad M_{n+1}=\max_{[a,b]}\left|f^{(n+1)}\right|.$$

Skal du bare ha feilen i ett punkt $x$, setter du inn $|\omega(x)|$ direkte og får et skarpere anslag.

Anslaget er en garanti, ikke en prognose. Den virkelige feilen er som regel klart mindre. Ber oppgaven om «et anslag på feilen», er det denne du skal levere — med begge maksimeringene skrevet ut.

For tre ekvidistante noder er

$$\max\left|\omega(x)\right|=\frac{2}{9\sqrt{3}}\,h^{3}\approx 0{,}1283\,h^{3}$$

over hele nodeintervallet. Det viktige er potensen: halverer du $h$, faller feilgrensen med en faktor 8 for tre noder, og mer for flere.

Ekvidistante noder er det naturlige valget når dataene kommer fra en tabell eller en logger. De er derimot ikke det beste valget hvis du står fritt — se Chebyshev-punktene under.

At interpolasjon med mange ekvidistante noder kan bli dårligere, ikke bedre, når du legger til punkter.

Bruker du 15 ekvidistante noder på en funksjon som $1/(1+25x^2)$, svinger interpolasjonspolynomet voldsomt nær kantene av intervallet, med feil som vokser med antall noder.

Årsaken ligger i nodepolynomet: med ekvidistante noder blir $|\omega(x)|$ dramatisk mye større nær kantene enn på midten. Feilformelen forutsier dette direkte.

Botemidlet er å flytte nodene tettere sammen mot kantene — det er nettopp det Chebyshev-punktene gjør.

$$x_k=\cos\left(\frac{(2k+1)\pi}{2n+2}\right),\qquad k=0,1,\dots,n.$$

De ligger tettere ved kantene enn på midten, og det er akkurat motsatt av hva intuisjonen sier — men det er der ekvidistante noder svikter.

Skal du bruke dem på et annet intervall $[a,b]$, flytter du dem med den samme transformasjonen som i kap. 6.2:

$$x=\frac{(b-a)\xi+(a+b)}{2}.$$

Punktene står på det utdelte formelarket — tren oppslaget på indekseringen: det er $2k+1$ i telleren og $2n+2$ i nevneren, og $n$ er graden, ikke antall punkter.

Begrunnelsen er entydigheten. La $p$ være interpolasjonspolynomet, og sett $q(x)=-p(-x)$. Da har $q$ samme grad som $p$, og i hver node er

$$q(x_k)=-p(-x_k)=-f(-x_k)=-\left(-f(x_k)\right)=f(x_k).$$

Altså interpolerer $q$ de samme punktene. Etter entydighetssetningen er $q=p$, det vil si $-p(-x)=p(x)$ — og det er nettopp definisjonen på at $p$ er odde.

Konsekvensen er praktisk: alle koeffisienter foran like potenser er null. Med fire symmetriske noder vet du på forhånd at $p(x)=ax^3+bx$, og du har spart deg halve regningen.

Det samme argumentet med like data ($f(-x_k)=f(x_k)$) gir at $p$ er like: bare like potenser overlever.

$$p_{\text{ny}}(x)=p(x)+f[x_0,\dots,x_{n+1}]\,(x-x_0)(x-x_1)\cdots(x-x_n).$$

Det gamle polynomet står uendret; du regner bare én ny dividert differanse og legger på ett ledd.

Dette er Newton-formens hovedgrunn til å eksistere. På Lagrange-form måtte alle $n+1$ kardinalfunksjonene bygges om fra bunn, fordi hver av dem inneholder alle nodene.

$$\begin{pmatrix}1&x_0&\cdots&x_0^{n}\\ 1&x_1&\cdots&x_1^{n}\\ \vdots& & &\vdots\\ 1&x_n&\cdots&x_n^{n}\end{pmatrix}\begin{pmatrix}c_0\\ c_1\\ \vdots\\ c_n\end{pmatrix}=\begin{pmatrix}f_0\\ f_1\\ \vdots\\ f_n\end{pmatrix}.$$

Matrisen kalles Vandermonde-matrisen, og den er inverterbar nettopp fordi nodene er forskjellige — det er en annen måte å se entydigheten på.

Metoden er korrekt, men treg, og den blir numerisk ustabil for mange punkter. Den er verdt å kjenne fordi den forklarer hvorfor problemet har én løsning, og fordi et lite $2\times2$- eller $3\times3$-system av og til er den raskeste veien for hånd.

Feilformelen viser hvorfor det andre er farlig: nodepolynomet $\omega(x)=\prod(x-x_k)$ vokser raskt så snart $x$ kommer utenfor nodene, siden alle faktorene da vokser samtidig. Innenfor har faktorene ulike fortegn og delvis utlikner hverandre.

Praktisk regel: hold deg innenfor. Må du utenfor, si det, og bruk feilformelen til å vise hvor fort anslaget forfaller.