1. No category

Teknisk Beräkningsvetenskap I
Tema 4: Förbränningsstrategier för raketer
modellerade som begynnelsevärdesproblem
Del 1: Interpolation och kvadratur
Eddie Wadbro
2 december 2015
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(1 : 27)
Innehåll
I
Interpolation
I
I
I
Polynominterpolation
Splines
Kvadratur
I
I
Newton–Cotes kvadratur (mittpunkts, trapets och Simpsons formel)
Noggrannhet och feluppskattningar
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(2 : 27)
Interpolation
Vanlig uppgift: vi behöver skapa en
trevlig kurva som passerar ett antal
givna punkter (till exempel i datorgrafik)
◦
◦
◦
◦
◦
Interpolationsproblemet:
givet n + 1 par (xk , yk ), k = 0, 1, . . . , n, där xk 6= xl om k 6= l, hitta en
funktion f så att f (xk ) = yk
Funktionen f är en interpolant till punktmängden {(xk , yk )}nk=0
Ett klassiskt val: polynominterpolation
f (x) = a0 + a1 x + a2 x 2 + · · · + an x n
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(3 : 27)
Polynominterpolation
Polynomets koefficienter kan enkelt bestämmas om man skriver polynomet på
följande form:
p(x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 ) + . . .
+ bn (x − x0 )(x − x1 ) · · · (x − xn−1 )
Villkoren p(xk ) = yk , k = 0, . . . , n, ger Newtons interpolationformel:
y0 = b 0
y1 = b0 + b1 (x1 − x0 )
y2 = b0 + b1 (x2 − x0 ) + b2 (x2 − x0 )(x2 − x1 )
..
.
. . . + bn (xn − x0 ) · · · (xn − xn−1 )
yn = b 0 + . . .
Ett undertriangulärt linjärt ekvationssystem för koefficienterna b0 , . . . , bn .
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(4 : 27)
Polynominterpolation
Sats
Låt (xk , yk )nk=0 vara en godtycklig koordinatmängd, där alla xk är skilda.
Det finns ett entydligt polynom p av grad ≤ n så att
p(xk ) = yk
k = 0, . . . , n
Observera: Antalet koefficienter i polynomet = antalet punkter som ska
interpoleras.
I Matlab:
p = polyfit(x, y, n);
I
Vekotorerna x, y (längd n + 1) innehåller koordinaterna
I
n: polynomgrad
I
p: vektor med interpolationpolynomets koefficienter
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(5 : 27)
Interpolationsfel
Antag att vi till approximera funktionen f med en interpolant
n
on
Med hjälp av koordinaterna xk , f (xk )
skapar vi ett polynom p n
k=0
som interpolerar f i dessa punkter p n (xk ) = f (xk ), k = 0, 1, . . . , n
Minskar felet när n → ∞?
Ja: Om punkterna {xk } är “valda väl”.
Nej: För givna {xk } så kan felet växa utan gräns!
Sats (Weierstrass approximationssats)
Varje funktion f ∈ C (a, b) kan approximeras till varje önskad noggrannhet
med polynom
Bernsteinpolynom
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(6 : 27)
Polynominterpolation
Hur beter sig interpolanten?
Test: interpolera funktionen
f (x) =
1
1 + 25x 2
Den blå funktionen f interpoleras i n + 1 ∗ markerade punkter jämt
fördelade över x-axeln av det gröna pollynomet p av grad n
2
2
1.5
1.5
1
1
∗
0.5
0∗
−0.5
−1
∗
∗
−0.5
∗
0.5
∗
0
∗
0.5
∗
1
6 punkter n = 5
∗
∗ ∗
∗ ∗
∗ ∗
0∗ ∗
−0.5
−1
−0.5
0
0.5
1
11 punkter n = 10
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(7 : 27)
Polynominterpolation
Runges fenomen: Interpolation med polynom med jämt fördelade
punkter tenderar att generera en interpolant som visar oscilationer nära
interpolationsintervallets kanter. Dessa oscillationer blir kraftigare när
polynomgraden ökar
Slutsats:
I
Interpolation med polynom av hög grad är ofta en fruktansvärd ide!
I
Genererar ofta stora svängningar mellan interpolationspunkterna
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(8 : 27)
Polynominterpolation: Två möjliga lösningar
Lösning 1:
I
Flytta interpolationspunkterna så att de är mer koncentrerade nära
kanterna
I
Ett bra val: Chebychevpunkter
Lösning 2: (det vanligaste tillvägagångssättet!)
I
Skapa en funktion (spline) genom att sammanfoga polynom av låg
grad
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(9 : 27)
Lösning 1: interpolera i Chebyshevpunkterna
1
Om vi vill använda ett polynom av grad n
för att interpolera på intervallet [−1, 1], så
använd punkterna
xk = cos
kπ
,
n
0.5
0 ∗∗ ∗ ∗
k = 0, . . . , n
−1
1
1
0.5
0.5
∗
∗
∗
−0.5
∗
0
∗ ∗ ∗∗
0.5
1
∗
1
∗ ∗
∗
∗
∗
∗
0∗
−1 −0.5
0.5
∗∗
0
0.5
1
∗ ∗
0 ∗∗
−1 −0.5
6 punkter n = 5
∗
∗
∗ ∗ ∗∗
0
0.5
∗∗ ∗
0 ∗∗∗
∗
∗
∗∗
−1 −0.5
1
11 punkter n = 10
0
∗∗∗∗∗
0.5
1
21 punkter n = 20
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(10 : 27)
Lösning 2: splines
I
I
I
I de flesta fall kan vi inte välja våra
interpolationspunkter! Ex.: De kan vara
givna av mätningar gjorde med ett
instrument som generarar data i N Hz.
Den vanligaste interpolationsmetoden:
splines: styckvisa polynom av låg grad.
Lämpligare än polynominterpolation i de
flesta fall
Enklaste valet, en linjär spline är en
kontinuerlig, styckvis linjär interpolant
4
3
2
1
0
0
1
2
3
4
Definition: En spline är en funktion som består av styckvisa polynom av
grad m sådan att den är kontinuerligt deriverbar m − 1 gånger
Vanligaste valet, förutom linjära splines: kubiska splines (CAD-standard)
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(11 : 27)
Kubiska splines
I
I
I
Antag att vi har n + 1 talpar (xk , yk ),
k = 0, . . . , n, sorterade så att xi < xj
om i < j
Den globala funktionen s definieras
styckvis på intervallen [xk−1 , xk ],
k = 1, . . . , n
För k = 1, . . . , n, bestäm de n kubiska
funktionerna
(k)
(k)
(k)
(k)
sk (x) = a0 + a1 x + a2 x 2 + a3 x 3
4
3
s4(x)
2
s (x)
2
s1(x)
s3(x)
1
0
0
1
2
3
4
på intervallen [xk−1 , xk ]
I
Det finns 4n koefficienter
I
Vi behöver således 4n ekvationer (villkor) för att bestämma dessa
koefficienter
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(12 : 27)
Kubiska splines
4
I
Interpolering i ändpunkterna för varje
interval
(
sk (xk−1 ) = yk−1
k = 1, . . . , n
sk (xk ) = yk
s4(x)
2
s (x)
2
s1(x)
s3(x)
1
0
0
2n villkor
I
3
1
2
3
4
Kontinuerliga första och andra derivator där två kubiska polynom
sätts samman
(
0
sk0 (xk ) = sk+1
(xk )
i = 1, . . . , n − 1
00
00
sk (xk ) = sk+1 (xk )
2(n − 1) villkor
I
Totalt 2n + 2(n − 1) = 4n − 2 villkor. Ytterligare två villkor behövs!
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(13 : 27)
Kubiska splines
Det finns många val för de två sista villkoren:
I “Not-a-knot” spline. Default i Matlab’s spline.
Lägger till villkoret om att tredje derivatan i x1
och xn−1 är kontinuerlig (s1000 (x1 ) = s2000 (x1 ) och
000
sn−1
(xn−1 ) = sn000 (xn−1 )).
I “Naturlig spline”. Kräv att kurvans krökning är
noll i ändpunkterna
(s100 (x0 )
=
sn00 (xn )
= 0).
4
3
2
1
natural
prescribed slopes
gL = gR = 0
0
−1
I Påtvinga givna lutningar gL , gR vid
ändpunkterna (s10 (x0 ) = gL och sn0 (xn ) = gR ).
−2
0
not-a-knot
1
2
3
Observera att splines är “globala”: små förändringar (exempelvis villkoren vid
ändpunkterna) kan påverka funktionens utseende överallt!
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(14 : 27)
Numerisk integration (Kvadratur)
Kvadrator, också kallat numerisk integration, består i att numerisk
beräkna integralen
Z b
I (f ) =
f (x) dx,
a
där a, b ∈ R.
När behövs numerisk integration?
Z
Z
dx
2
Ex.:
e −x dx
ln x
Z p
1 + cos 2 xdx
Saknar primitiva funktioner (iallafall i form av elementära funtioner)
(Liouville, Liouville–Risch)
(En elementär funktion är en funktion som kan uttryckas med ändligt många
algebraiska operationer (+, -, *, /), konstanter, exponentialfunktionen,
sammansättningar av elementära funktioner, samt inverser till elementära funktioner.)
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(15 : 27)
4
Kvadratur
Ex.: Vi vill beräkna integralen
b
Z
I =
2
e −x dx
a
I
Eftersom integranden är “enkel” kan vi i Matlab definera den som
som anonym funktion och därefter approximera integralen
func = @(x) exp(-x.*x);
I = quad(func, a, b);
Variablen I kommer att innehålla en uppskattning av
Rb
a
2
e −x dx
I
Observera: Den funktion som använd i quad måste acceptera en
vektor som inargument och då returnera en vektor med motsvarande
funktionsvärden!
I
Om integranden är komplicerad, så är det oftast smidigare att istället
definiera den i en m-fil
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(16 : 27)
Grundiden bakom många kvadraturmetoder
Integraler av polynom kan beräknas exakt!
Approximera f (x) ≈ P(x), där polynomet P interpolerar f i punkterna.
Då
Z b
Z b
f (x)dx ≈
P(x)dx.
a
a
Den beräkningsbara integralen till höger approximerar den potentiellt
icke-beräkningsbara integralen till vänster
I
Hög noggrannhet kräver polynom av hög grad
I
Alternativt, dela upp intervalet [a, b] i många delintervall där
funktionen approximeras med polynom av låg grad!
(styckvis approximation)
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(17 : 27)
Newton–Cotes kvadratur
Approximera integranden med ett polynom P av grad N,
Z
b
Z
f (x)dx =
a
b
P(x)dx +
R
|{z}
a
rest term
(fel)
med trunkeringsfel
b
Z
(x − x0 )(x − x1 ) . . . (x − xN )
R=
a
f N+1 (ξ(x))
dx
(N + 1)!
Öppna Newton–Cotes metoder:
Inkludera intervalländpunkterna bland interpolationspunkterna
N = 1:
N = 2:
Trapetsformeln
Simpsons formel
Slutna Newton–Cotes metoder:
Intervalländpunkterna tas ej med bland interpolationspunkterna
N = 0:
Mittpunktsformeln
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(18 : 27)
Mittpunktsformeln
I
Dela intervallet [a, b] i n delintervall
(I illustrationen till höger är n = 6.)
I
I varje delinterval, approximera f med
dess värde i intervallmittpunkten
I
Summera rektanglarnas areor
x0
x1
x2
x3
x4
x5
x6
Detta ger oss att
n−1
X
b
Z
f (x)dx ≈
a
(xk+1 − xk )f
k=0
x + x
k
k+1
2
Om vi använder en ekvidistant uppdelning, xk+1 − xk = ∆x för alla k, så
Z
b
(∆x)
f (x)dx ≈ IM
a
(a, b) = ∆x
n−1
X
där
f (xk+1/2 ),
xk+1/2 =
k=0
xk + xk+1
2
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(19 : 27)
Trapetsformeln
I
Dela intervallet [a, b] i n delintervall
(I illustrationen till höger är n = 6.)
I
I varje delintervall approximera f med en
kontinuerlig, styckvis linjär funktion
I
Summera parallelltrapetsernas areor
x0
x1
x2
x3
x4
x5
x6
Detta ger oss att
b
Z
f (x)dx ≈
a
n−1
X
(xk+1 − xk )
k=0
f (xk ) + f (xk+1 )
2
Om vi använder en ekvidistant uppdelning, xk+1 − xk = ∆x för alla k, så
Z
a
b
(∆x)
f (x)dx ≈ IT
(a, b) =
∆x f (x0 )+2f (x1 )+2f (x2 )+. . .+2f (xn−1 )+f (xn )
2
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(20 : 27)
Simpsons formel
I
Dela intervallet [a, b] i N delintervall
(I illustrationen till höger är N = 3.)
I
I varje delintervall approximera f med en
kontinuerlig, styckvis kvadratisk funktion
I
Beräkna och summera integralerna av
interpolanten över varje delinterval
x0
x1
x2
x3
x4
x5
x6
Detta ger oss att
Z
b
f (x)dx ≈
a
N−1
X
(x2k+2 − x2k )
k=0
f (x2k ) + 4f (x2k+1 ) + f (x2k+2 )
6
Om vi använder en ekvidistant uppdelning, xk+1 − xk = ∆x för alla k, så
Z
a
b
(∆x)
f (x)dx ≈ IS
(a, b) =
∆x f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) . . .
3
+ 2f (xn−2 ) + 4f (xn−1 ) + f (xn )
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(21 : 27)
Noggrannhet:
Kvadraturfelet är ett diskretiseringsfel
Sats
För en två gånger kontinuerligt deriverbar f gäller
b
Z
a
b
Z
(∆x)
(a, b) +
∆x 2
(b − a)f 00 (ξ), för något ξ ∈ [a, b],
24
(∆x)
(a, b) −
∆x 2
(b − a)f 00 (ξ), för något ξ ∈ [a, b].
12
f (x) dx = IM
f (x) dx = IT
a
Fel: O(∆x 2 ): dessa metoder har ordning p = 2
För en fyra gånger kontinuerligt deriverbar f gäller
Z
a
b
(∆x)
f (x) dx = IS
∆x 4
(b − a)f (4) (ξ), för något ξ ∈ [a, b].
180
(a, b) −
Fel: O(∆x 4 ): metoden har ordning p = 4
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(22 : 27)
Feluppskattning
Informationen från satsen ovan kan användas för att uppskatta felet utan
att veta det exakta värdet på integralen!
b
∆x 2
(b − a)f 00 (ξ1 )
12
a
Z b
4∆x 2
(2∆x)
f (x) dx = IT
(a, b) −
(b − a)f 00 (ξ2 )
12
a
Z
(∆x)
f (x) dx = IT
(a, b) −
ξ1 ∈ [a, b]
(1)
ξ2 ∈ [a, b]
(2)
2
(∆x)
00
Antag att f 00 (ξ1 ) ≈ f 00 (ξ2 ) och låt ET
= − ∆x
12 (b − a)f (ξ1 ). Genom
att subtrahera uttryck (2) från uttryck (1) erhålls
(∆x)
ET∆x
=
IT
(2∆x)
(a, b) − IT
3
(a, b)
(“tredjedelsregeln” )
Alltså, genom att utföra två integralberäkningar med olika steglängd,
2∆x and ∆x, så kan vi uppskatta felet vid användning av steglängden ∆x
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(23 : 27)
Feluppskattning
En liknande analys av Simpons formel ger uppskattningen
(∆x)
ES∆x
=
IS
(2∆x)
(a, b) − IS
15
(a, b)
(“femtondelsregeln” )
Generellt, för en inegrationsformel med felterm O(∆x p ):
E ∆x =
I (∆x) (a, b) − I (2∆x) (a, b)
2p − 1
Dessa uppskattningar kan användas i en adaptiv process för att lokalt
förfina där det behövs eller i så kallade extrapolationsprocedurer
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(24 : 27)
Relation mellan trapetsformeln och Simpsons formel
(∆x)
IS
4 (∆x)
1 (2∆x)
(a, b) = IT (a, b) − IT
(a, b)
3
3
Richardsonextrapolation:
Diskretisering med steglängd h approximerar y (0) med y (h)
Antag att vi har en asymptotisk felutveckling
y (h) = y (0) + c1 hp1 + c2 hp2 + c3 hp3 + . . .
där p1 < p2 < p3 < . . .
Beräkna y (qh) och y (h) (vanligvis väljs q = 2). Då
y (h) =
y (qh) =
y (0) +
c1 hp1 + O(hp2 )
y (0) + c1 (qh)p1 + O(hp2 )
q p1 y (h) − y (qh) = (q p1 − 1)y (0)
+ O(hp2 )
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(25 : 27)
Adaptiv Simpson: en rekursiv algoritm
Vi vill ha
Z
a
b
f (x) dx − I (a, b) ≤ 0. Låt arbetsintervalet vara [a, b] och sätt ∆x = (b − a)/2
1. Beräkna IS∆x och IS2∆x på nuvarande arbetsinterval
2. Uppskatta felet med femtondelsregeln
längd
3. om det uppskattade felet > arbetsintervallets
b−a
I
I
I
(∆x)
Förkasta IS
Dela arbetsintervallet i två lika delinterval
För dessa fortsätt från punkt 1. med ∆x ← ∆x/2
Annars
I
I
(∆x)
Acceptera IS
Om det finns obehandlade delinterval, uppskatta integralen över
nästa, annar klar!
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(26 : 27)
Multipelintegraler
Z
f (x, y , z)dxdzdz →
Ω
X
wix wjy wkz f (xi , yj , zk )
i,j,k
Många funktionsevalueringar behövs (“the curse of dimensions”)
För väldigt många dimensioner kan behöva tillgå exempelvis Monte Carlo
metoder, statistiska uppskattningar
Eddie Wadbro, Tema 4: Interpolation och kvadratur, 2 december 2015
(27 : 27)