Aproximace funkcí Aproximace je výpočet funkčních hodnot funkce z nějaké třídy funkcí, která je v určitém smyslu nejbližší funkci nebo datům, která chceme aproximovat. Tř ída funkcí, ze které volíme aproximace je nejčastěji: 2 m ● polynom c 0 c 1 xc 2 x c m x ● zobecněný polynom c 0 g 0 x c 1 g 1 xc m g m x , kde g 1 x ,, g m x je systém m1 jednoduchých, LN a dostatečně hladkých funkcí a 0 a 1 xa 2 x 2 a k x k ● racionální lomenou funkci ( m=kl ) b 0 b 1 xb 2 x 2 b l x l ● jiná Důvody aproximací jsou: ● přílišná složitost funkce ● nutnost počítání characteristik funkce (derivací, integrálů) ● neznalost funkce (resp. jejího analytického vyjádření) Typy aproximací: ● Interpolace, extrapolace – hledáme funkci ( m x ), která má v zadaných bodech x 0 ,, x n stejné hodnoty (případně i jejich derivace), jako funkce, kterou chceme aproximovat, dělá se buď globálně nebo lokálně na každém intervalu 〈 x k , x k1 〉 . ● Čebyševova aproximace hledáme funkci ( m x ) tak, abychom na zadaném intervalu 〈 a , b 〉 minimalizovali funkcionál max x ∈〈 a , b〉∣ m x − f x ∣ , m x pak nazýváme nejlepší stejnoměrnou aproximací. ● Metoda nejmenších čtverců – je založena na minimalizaci buď b
∫ w x [ f x −m x ] dx nebo a
2
n
∑ w x i [ f x i −m x i ]
2
(podle toho, jestli se
i=0
jedná o spojitý nebo diskrétní případ), v diskrétním případě je vždy nm a tak m neprochází všemi zadanými body, volí se často w x=1 , pokud je w x ≠const , mluví se o vážené metodě nejmenších čtverců a funkci w x se říká váha.
Interpolace a extrapolace Lagrangeův interpolační polynom Jsou zadány body x 0 ,, x n a k nim funkční hodnoty y 0 , , y n . Lagrangeův interpolační polynom L n je polynom stupně nejvýše n takový , že se jeho hodnota ve všech zadaných bodech x 0 ,, x n rovná zadaným funkčním hodnotám y 0 ,, y n .
Platí tedy L n x 0 =y 0 , L n x 1 =y 1 ,, L n x n =y n a bodům x 0 ,, x n říkáme uzly interpolace. 1 i= j K sestrojení Lagrangeova polynomu se definují funkce F i x tak, že F i x j = 0 i≠ j Tyto funkce mají tvar n x−x 0 ⋯ x−x i−1 ⋅ x−x i1 ⋯ x−x n x−x j F i x = = ∏ x i −x 0 ⋯ x i −x i−1 ⋅ x i −x i1 ⋯ x i −x n j=0, j≠i x i −x j a polynom L n x pak má tvar n
L n x =∑ y i⋅F i x . i=0
Pokud jsou uzly rozmístněny ekvidistantně (rovnoměrně s krokem h ), lze Lagrangeův polynom zapsat jako n t⋅t−1⋯t−n L n x=L n x 0 th=∑ y i . t−ii !n−i!−1n−i i=0 Lagrangeův vzorec se nehodí pro numerický výpočet interpolace. Koeficienty Lagrangeova polynomu lze spočítat řešením systému lineárních rovnic. Matice tohoto systému (Van der Mondova) je však (pro ekvidistantní uzly) často špatně podmíněná, vypočet koeficientů tedy není přesný. Pro implementaci Lagrangeovy interpolace se používá Nevillův algoritmus. Principem tohoto algoritmu je Lagrangeova interpolace na postupně se zvětšujícím počtu uzlů. L ik x=L i x , x i , x i−1 ,, x i−k , y i , y i−1 ,, y i−k je tedy interpolace Lagrangeovým polynomem stupně nejvýše i v k bodech x i ,, x i−k . L i0 x =y i , L i−1,0 x =y i−1 , x i −x L i−1,0 x − x i−1−x L i0 x L i1 x=L i x , x i , x i−1 , y i , y i−1 = = x i −x i−1 , atd. x i −x y i−1− x i−1−x y i = x i −x i−1 Pro tyto postupné aproximace platí totiž rekurentní vztah L ik x =
x i −x L i−1, k−1 x− x i−1−x L i , k−1 x x i −x i−k
x
Interpolace se tedy provádí podle následujícího schématu y k=0 k=1 k=i1 k=i
x0
y0
L 00
x1
y1
L 10
L 11
⋮
⋮
⋮
⋮
⋱
x i−1
y i−1
L i−1,0
L i−1,1
L i−1, i−1
xi
yi
L i0
L i1
L i , i−1
L ii
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋱
xn
yn
L n0
L n1
L n , i−1
Lni
k=n
Ln n
Polynomy se počítají zleva doprava a zezhora dolu (ve směru tmavnoucí barvy). Lagrangeův interpolační polynom je pak L n x =L n n x . Odhad chyby aproximace je dán buď rozdílem L n n x −L n−1, n−1 x nebo rozdílem L n n x −L n , n−1 x . Příklad v PASCALU DEMPOL.PAS . http://www.math.odu.edu/~bogacki/videnum/lagran.htm Při implementaci Nevillova algoritmu se počítají koeficienty C m , i =L im , m −L im−1, m−1 , D m , i =L im , m −L im , m−1 a x i −x C m , i1 −D m , i x im1 −x C m , i1−D m , i C m1, i = , D m1, i = x i −x im1 x i −x im1 Přitom v koeficientech C m , i a D m , i je odhad chyby. Odhad chyby polynomiální interpolace je dán n1
f , R x = f x −L n x = x−x 0 ⋯ x−x n n1! kde ∈〈 min x 0 ,, x n , max x 0 ,, x n 〉 .
Interpolace je hledání hodnoty funkce v bodě x ∈〈 x 0 , x n 〉 . Extrapolace je hledání hodnoty funkce v bodě xx 0 nebo xx n . Při extrapolaci při rostoucí vzdálenosti od krajních bodů intervalu chyba roste rychle. Pro ekvidistantní uzly a vysoké stupně polynomu, má polynom mezi uzly tendenci oscilovat. Proto se polynomiální interpolace s ekvidistantními uzly používá v praxi pro n≤7 .
Konvergence L n x f x s rostoucím n je zaručena pro analytické funkce, které mají derivaci v celém komplexním oboru s vyjímkou ∞ . Newtonův interpolační polynom Obdoba Lagrangeova interpolačního polynomu. N n x=a 0 a 1 x−x 0 a 2 x−x 0 x−x 1 a n x−x 0 ⋯ x−x n
Koeficienty a 0 ,, a n se počítají pomocí poměrných nebo obyčejných diferencí. Poměrná diference 1.řádu f x i , x i1 = ktého řádu f x i ,, x ik =
y i1 −y i
, x i1 −x i f x i1 ,, x ik − f x i ,, x ik−1 x ik −x i
Potom a 0 = f x 0 , a1 = f x 0 , x 1 ,, a n = f x 0 ,, x n Obyčejná diference 1. řádu 1 f i =y i1 −y i ktého řádu k f i =k−1 f i1 −k−1 f i Příklad v PASCALU ZKPOL.PAS , ZKPO1L.PAS , ZKPOL2.PAS .
Racionální lomená funkce
Pro některé funkce vhodnější než polynom. Je dáno m1 dobů, pak lze zadanou funkci aproximovat pomocí racionální lomené 2 i p 0 p1 x p 2 x pi x funkce R i , j x = , kde i j=m a lze q 0 =1 . 2 j q 0 q 1 xq 2 x q j x Algoritmus obdobný Nevillovu algoritmu – Burlisch a Stoer.
RATINT: procedure expose A. B. parse arg N, V Tiny = 1E25; Ns = 1; Hh = ABS(V A.1) do I = 1 to N H = ABS(V A.I) if H = 0 then do Y = B.I; Dy = 0 return Y Dy end else if H < Hh then do; Ns = I; Hh = H; end C.I = B.I; D.I = B.I + Tiny end Y = B.Ns; Ns = Ns 1 do M = 1 to N 1 do I = 1 to N M Ip1 = I + 1; W = C.Ip1 D.I IpM = I + M; H = A.IpM V T = (A.I V) * D.I / H Dd = T C.Ip1 if Dd = 0 then call ERROR "RATINT: Error ", "the interpolating function has", "a pole at the requested value V" Dd = W / Dd; D.I = C.Ip1 * Dd; C.I = T * Dd end if 2 * Ns < (N M) then do; Nsp1 = Ns + 1; Dy = C.Nsp1; end else do; Dy = D.Ns; Ns = Ns 1; end Y = Y + Dy end return Y Dy ERROR: say ARG(1); exit
Hermiteova interpolace
Kromě funkčních hodnot zadané i některé derivace, tedy x 0 ,, x n , y 0 ,, y n a a a y '0 , y ''0 ,, y 0 , , y 'n , y ''n ,, y n . 0
n
n=ma 0 a n pak podmínky dané zadáním splňuje funkce H n x , který lze zapsat jako H n x =L n x x−x 0 ⋯ x−x n H n−m−1 x . Pro polynom H n−m−1 x stupně n−m−1 platí y 'i −L 'n x i H n−m−1 x i = a tak se dopočítají pomocí x i −x 0 ⋯ x i −x i−1 x i −x i1 ⋯ x i −x n zadaných prvních derivací některé koeficienty. Pro další derivace se postupuje podobně. http://www.math.odu.edu/~bogacki/videnum/hermit.htm
Interpolační spline
Lokální interpolace taková, že kromě průchodu všemi uzly požadujeme v uzlech ještě spojitost první derivace (alespoň první). Máme tedy na funkci, aproximující zadané hodnoty na intervalu x i , x i1 4 podmínky y i , y i1 , y 'i =y 'i− , y 'i1 =y 'i1− . Proto funkce, použitá pro interpolaci musí mít nejméně 4 volitelné parametry. Velmi často se jako tato funkce používá kubický polynom – kubický spline. Mámeli tedy kubický polynom a 2 krát ho zderivujeme dostaneme lineární funkci. Pro tuto funkci máme zadány hodnoty v krajních bodech y ''i , y ''i1 . Lagrangeovou interpolací pak můžeme vyjádřit y ''= Ay ''i By ''i1 , kde A=
x i1 −x x i1 −x i
a B=
x−x i x i1−x i
.
Po dvojím integrování s podmínkami y x i =y i a y x i1 =y i1 dostaneme y= A y i B y i1C y ''i D y ''i1 , 1 3 2 kde C= A − A x i1−x i 6 1 3 2 a D= B −B x i1−x i . 6 (Integrují se koeficienty A x a B x ) Okrajové podmínky platí zkuste si dosadit do A x a B x body x i a x i1 . Hodnoty y ''i můžeme vypočítat z podmínky spojitosti 1. derivace a ze znalosti 1. derivací v krajních bodech intervalu ( tedy y '0 a y 'n ). Pokud derivace neznáme, zadává se y ''0 =0 a y ''n =0 při rozený spline. Provedením derivace y x a z podmínky spojitosti 1. derivace v krajních bodech (tj. pokud dosadíme do y ' x na intervalu 〈 x i−1 , x i 〉 bod x i , musí se to rovnat tomu, když se do y ' x na intervalu 〈 x i , x i1 〉 dosadí rovněž bod x i ). Z toho plyne následující soustava pro druhé derivace: x i −x i−1 6
y ''i−1
x i1−x i−1
y ''i
x i1−x i
y ''i1 =
y i1 −y i
3 6 x i1 −x i To je soustava s tridiagonální maticí.
−
y i −y i−1 x i −x i−1
Příklady v PASCALU SPLINE.PAS , SPLPRA.PAS , SPLPRB.PAS , SPLRUN.PAS .
Více dimenzí (2) lokální a globální
●
Lineární
Jsou zadány y1=y j , k , y 2 =y j1, k , y 3=y j1, k1 , y 4 =y j , k1 Definujeme t =
x 1 −x 1 j
a u=
x 2 −x 2 k
D1 úseček ohraničujících interval).
D2
(tedy v podstatě relativní vzdálenosti od
Lokální lineární interpolace je dána vztahem y x 1 , x 2 =1−t 1−u y1t 1−u y 2 tuy 31−t uy 4 .
●
Bikubická
Jde o lokální interpolaci Hermiteova typu. V každém bodě jsou zadány hodnoty fce. a jejích parciálních derivací. ∂y ∂y ∂2 y y x 1 , x 2 , , , ∂ x1 ∂ x2 ∂ x1 ∂ x2 V každém bodě jsou tedy dány 4 podmínky a ve dvourozměrném intervalu je jich tedy 16 Bikubická interpolace se hledá ve tvaru 4
4
y x 1 , x 2 =∑ ∑ c ij t i−1 u j−1 . i=1 j=1
●
Bikubický spline
Lokální interpolace se spojitými parciálními derivacemi v obou směrech na hranicích intervalu. Interpoluje se v jednom směru a pak ve druhém.
●
Globální interpolace
Vyšší řád přesnosti. Opět interpolace v jednom směru a pak ve druhém.
