MFF UK / Ústav teoretické fyziky / Tomáš Ledvinka
Přednášky
. . . . . . . . . . . . . . . . . . . . . . . . . .
Programování pro fyziky (1.r)
  Úvod
  Programy a algoritmy
  Píšeme program
  Píšeme program II
  Procedury a funkce
  Malujeme funkce
  Chyby. Typy I.
  Typy II. Pole a Záznamy
  Pole II.Řetězce.Soubory.
  Gnuplot.Interpolace...
    Funkce v tabule
    Interpolace
    Polynomy
  Matice. Velké O...
  Fronta,Zásobník. Postscript
  Bin. soubory, ...
  Ukazatele,Objekty, ...
Počítačová algebra
Klasická elektrodynamika (2.r)
Vybrané partie OTR

Cvičení
. . . . . . . . . . . . . . . . . . . . . . . . . .
Programování pro fyziky (1.r)
Teoretická mechanika (2.r)
Klasická elektrodynamika (2.r)


Věda
. . . . . . . . . . . . . . . . . . . . . . . . . .
Diskové zdroje v OTR
Hyperbolické systémy v OTR


Kontakt
. . . . . . . . . . . . . . . . . . . . . . . . . .
Email
Konzultační hodiny


Ostatní
. . . . . . . . . . . . . . . . . . . . . . . . .
Mallorca
Ze společnosti

Interpolace

V výše uvedeném případě je rozložení hodnot nezávislé proměnné x rovnoměrné, tzv. ekvidistantní. Velmi často se s hodnotami takto rovnoměrně poskládanými setkáme u veličin měřených v závislosti na čase (např. vzorky zvuku na CD...)  Data v proměnných HonotyX a HodnotyY můžeme znázornit grafem

 

Jak máme ale určit funkční hodnotu mezi jednotlivými body grafu? Co třeba proložit sousedními body úsečky?

 

Vidíme, že výsledek není příliš dokonalý, ale přesto

Cvičení: (důležité) napište funkci, která mi pro reálné x v rozsahu -3 .. 3 vrátí takto vypočtenou hodnotu a bude mít hlavičku např.

function FzTabulky(x : real) : real;

a namalujte její graf. Pokud se ke cvičení dostanete později, nezapomeňte, že v uspořádaném poli se dá hledat velmi rychle.

Pokud se nám tato hranatá funkce nezdá dost dobrá, můžeme zkusit lepší. Víme, že N+1 body můžeme proložit právě jeden polynom N-tého stupně, takže máme hned kandidát. A kolika body budeme polynom prokládat? No přece všemi, když už je máme. Tady je výsledek:

Že nevypadá stále dost dobře? No tak zkusíme někde sehnat víc bodů a to musí pomoci

Nepomohlo...

Přes toto počáteční varování jsou polynomiální interpolace velmi užitečné, a tak si o nich povíme více.
V příkladu byla ale použit lehce zákeřná funkce arcustangens a navíc aspoň uprostřed intervalu vypadá výsledek hezky.

Především problém je lineární a tak stačí zjistit jak vypadá funkce proložená implusem a tu pak oškálujeme a sečteme přes všechny vzorky. Takto vypadá impuls v x=5:


Funkce proložená vyznačenými puntíky bude určitě úměrná polynomu
Y5(x) = (x - 1)(x - 2)(x - 3)(x - 4)(x - 6)(x - 7)(x - 8)(x - 9)(x - 10)
Pokud ji vydělíme Y5(5) máme funkci na obrázku.

Takovouto úvahou dostáváme Lagrangeův interpolační vzorec


Pokud si povšimneme jmenovatelů, tak vidíme, že pro obecné rozložení bodů xi budeme muset provést N*(N-1) násobení. Protože ale není důvod předpokládat, že N bude nabývat nějakých závratných výšek, nemluvíme zde ani tak o náročnosti výpočtu jako o faktu, že kód bude obsahovat dva cykly. Samozřejmě existují i jiná schémata pro výpočet interpolace.

Cvičení: Kromě obecné procedury pro interpolaci polynomem obecného stupně, si zkuste  napsat některou z řady užitečných interpolačních funkcí jako např:

function Interp2(x, x1,x2,y1,y2 : real) : real;
function Interp3(x, x1,x2,x3,y1,y2,y3 : real) : real;

atp. které určitě někde využijeme.

Příklad:

toto je funkce pro výpočet interpolované hodnoty přímo z Lagrangova vzorečku:

function LInterp(t:real; const X,Y : array of real):real;
var i,j,n : integer;
s,f : real;
begin
  n := High(X);
  {assert( n = High(Y) ); tam dame az budeme vedet co to je}

  s := 0;
  for i := 0 to n do begin
  f:=1;
  for j := 0 to n do if i<>j then f := f*(t-X[j])/(X[i]-X[j]);
    s := s+f*Y[i];
  end;
  LInterp := s;
end;

přičemž jde o přímý přepis vzorečku, bez jakýchkoli triků.

.