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

Polynomy jako pole koeficientů

Protože teď umíme spočíst fuknční hodnotu polynomu proloženého body Xi,Yi, můžeme chápat tuto dvojici vektorů také jako zápis polynomu. Obvyklé ale polynomem rozumíme spíše pole jeho koeficientů.

Následující program ukazuje jak definujeme pro takový polynom základní operace a jak  je můžeme použít ke konstrukci koeficientů interpolačního polynomu.  Jde o  rozdílné obou přístupy a nakonec vykreslíme jejich rozdíl. (Komentář: Spočíst nejdříve koeficienty interpolovaného polynomu a pak používat Hornerovo schéma k výpočtu hodnot interpolované funkce není nejlepší.)

Program Lagr2;
Uses GnuPlot01;

type tPolynom=array of real;
{
  Pozor formalni parametry funkci jsou typu array of real
  a nikoli tPolynom, abych mohl psat
  W := SoucinPolynomu(A,[-1,0,1]); pro W:=(x^2-1)*A
}

Function HodnotaPolynomu(const P:array of real;x:real): real;
var s : real;
    i : integer;
{Spocte funkcni hodnotu polynomu pomoci Hornerova schematu}
begin
  s := P[High(P)];
  for i := High(P)-1 downto 0 do s:=s*x+P[i];
  HodnotaPolynomu := s;
end;


Function KopiePolynomu(const P:array of real): tPolynom;
var Q : tPolynom;
    i : integer;
begin
  SetLength(Q,High(P)+1);
  for i := 0 to High(P) do Q[i]:=P[i];
  KopiePolynomu := Q;
end;

Function SoucetPolynomu(const A,B:array of real): tPolynom;
var Q : tPolynom;
    i,n : integer;
    s : real;
begin
  n:=High(A);
  if nj then f := SoucinPolynomu(f,[-X[j]/(X[i]-X[j]),1/(X[i]-X[j])]);
     s := SoucetPolynomu(s,f);
   end;
   InterpolacniPolynom := s;
end;



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) );

   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;



const N = 12;
type tIndexTabulky = 0 .. N;
     tPolicko = array[tIndexTabulky] of real;
const HodnotyX : tPolicko = (-3,-2.5,-2.0,-1.5,-1.0,-0.5, 0,0.5,1.0,1.5,2.0,2.5,3.0);
var   HodnotyY : tPolicko ;

const M = 1000;
var iX,iY,iDelta : array [0..M] of real;
    var i : integer;
    var xa,xb:real;
    IP : tPolynom;
begin

  HodnotyY[6]:=1;

  MalujPole(HodnotyX,HodnotyY,'Hodnoty','title "Data" with points pointtype 2 pointsize 2');
  //MalujPole(HodnotyX,HodnotyY,'>>Hodnoty',' smooth cspline title "cspline"');

  // nyni spoctu nejdriv interpolacni polynom
  IP := InterpolacniPolynom(HodnotyX,HodnotyY);

  // potrebuji vzit nasledujici rozsah promenne x
  xa := HodnotyX[Low(HodnotyX)];
  xb := HodnotyX[High(HodnotyX)];

  // a do poli nacpu hodnoty x_i a  y_i = IP(x_i)
  for i :=0 to M do begin
     iX[i] := xa + i/M*(xb-xa);
     iY[i] := LInterp(iX[i],HodnotyX,HodnotyY);
     iDelta[i] := HodnotaPolynomu(IP,iX[i])-iY[i];
  end;

  MalujPole(iX,iY,'>>Hodnoty','title "Lagrange" with lines');
  MalujPole(iX,iDelta,'Rozdily','title "Chyba IP" with lines');

  ZobrazGraf('Hodnoty','[][-1:2]');
  ZobrazGraf('Rozdily');

  { Nasledujici prikaz nas muze presvedcit, ze problem je ve vypoctu polynomu:
   Writeln(HodnotaPolynomu(IP,-3),' ',LInterp(-3,HodnotyX,HodnotyY));
   readln;
  }
end.

Cvičení: Vyzkoušejte jaké obrázky program vykreslí. Měňte počet bodů a prokládané funkce a pozorujte co se stane.

Cvičení: V jedné z minulých přednášek jsme použili rekurentní vztah pro Legendrovy polynomy. Zkuste jej pomocí funkcí uvedených v programu použít k výpočtů koeficientů Pn(x).

Cvičení: Jak je to s přesností pro polynomy vyššího řádu. Porovnejte podobně jako ve výše uvedeném programu pro interpolace.

.