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...
  Matice. Velké O...
  Fronta,Zásobník. Postscript
  Bin. soubory, ...
  Ukazatele,Objekty, ...
    Ukazatele
    Přetypování
    Num. kvadratura
    Objekty
    Přehled probraných témat
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

Numerická kvadratura ( výpočet určitého integrálu funkce jedné proměnné)

Nejdříve si ukážeme jak spočítat tzv. lichoběžníkovým pravidlem přibližnou hodnotu určitého integrálu. Řekněme, že počítáme

Pro rozdělení integračního intervalu na 2, 4, 8 a 16 podintervalů (kde funkci nahradíme lineární interpolací) vypadá lichoběžníkové pravidlo takto:

Protože určitý integrál je lineární zobrazení z prostoru funkcí na intervalu <a,b> do R, nepřekvapí, že lichoběžníkové pravidlo se redukuje na lineární kombinaci funkčních hodnot, konkrétně

Následující program spočítá podle tohoto schématu jak výše uvedený intergrál tak dvojnásobek plochy pod půlkružnicí sqrt(1-x^2).

program TrapInt;

type tRFunkce = function (x:real) : real;

function LichobeznikovePravidlo(f: tRFunkce; a,b : real; N:integer):real;
var dx : real;
    s_kraj,s_vnitr : real;
    i : integer;
begin
  dx := (b-a)/N;
  s_kraj  := f(a) + f(b);

  s_vnitr := 0;
  for i := 1 to N-1 do s_vnitr := s_vnitr + f(a+dx*i);

  LichobeznikovePravidlo:=(s_vnitr+s_kraj/2.0)*dx;
end;

function fce(x:real):real;
begin
  fce := sin(x);
end;

function gce(x:real):real;
begin
  gce := sqrt(1-x*x);
end;

const nmax=200000;
var n: integer;

begin
  n:=2;
  while n<nmax do begin
    Writeln(n, LichobeznikovePravidlo(fce, 0,Pi/2, n):18:14);
    n:=n*2;
  end;
  Writeln;Writeln;

  n:=2;
  while n<nmax do begin
    Writeln(n, 2*LichobeznikovePravidlo(gce, -1,1, n):18:14);
    n:=n*2;
  end;
  Readln;
end.

Zde je výstup programu:
2 0.94805944896852
4 0.98711580097278
8 0.99678517188617
16 0.99919668048507
32 0.99979919432002
64 0.99994980009210
128 0.99998745011753
256 0.99999686253529
512 0.99999921563419
1024 0.99999980390857
2048 0.99999995097714
4096 0.99999998774429
8192 0.99999999693607
16384 0.99999999923402
32768 0.99999999980850
65536 0.99999999995213
131072 0.99999999998803


2 2.00000000000000
4 2.73205080756888
8 2.99570906810244
16 3.08981914435717
32 3.12325303782774
64 3.13510242287713
128 3.13929691277968
256 3.14078079239662
512 3.14130558295723
1024 3.14149115271966
2048 3.14155676653901
4096 3.14157996541144
8192 3.14158816760775
16384 3.14159106754971
32768 3.14159209283888
65536 3.14159245533423
131072 3.14159258349584


Tento výstup vhodně zpracován gnuplotem dá graf závislosti chyb na velikosti intervalu dx.

Na přednášce:

  • Komentář o hezkých a ošklivých funkcích.
  • Koměntář o sčítání (-1)k-1/k , aneb jak co nejsložitěji spočíst ln(2).
    K tomu potřeba následující obrázek, kde se má ilustrovat, že sečny konvergují k tečně v x=0 a má-li
    funkce Taylorův rozvoj v x=0 musí zelené body jít k limitě kvadraticky atd.

Lichoběžníkové pravidlo splňuje tři základní valstnosti integrálů, lineariatu při násobení konstantou, translační nezávislost a škálování při lineární transformaci měnící šířku intervalu. To znamená, že můžeme zkoumat, jak umí integrovat funkce x^k v intervalu 0..1 a výsledky pak zobecnit pro všechny polynomy a intervaly. Při rozdělení na N stejných intervalů jsou hodnoty spočtené lichoběžníkovým pravidlem uvedeny v tabulce.

1   1
x  
x^2  
x^3  
x^4  
x^5  
x^6  
x^7  
x^8  
 

Proto tzv. inženýrskou indukcí dostáváme, že pro každý polynom f bude platit, že

I( f , dx ) = SpravnaHodnotaIntegralu + a*q + b*q^2 + ...

kde q = 1/N^2 případně q = dx^2 (na konstantě nesejde)

Ilustrace na přednášce: likvidace členu 1/N^2.

Nyní vezměme vektor sedmnácti funkčních hodnot f(x) vzniklých dělením integračního intervalu na šestnáct podintervalů:
[f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16].
Lichoběžníkové pravidlo odpovídá skalárnímu součinu tohoto vektoru (ekvidistantních) funkčních hodnot s vektorem
dx.[1/2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/2] (a budeme mu odteď říkat forma).
Lichoběžníkové pravidlo pro dvojnásobný krok odpovídá formě
(2 dx).[1/2,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1/2] a pro čtyřnásobný krok je to
(4 dx).[1/2,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1/2] atd.

Provedeme-li v grafu závislosti I na q lineární extraterpolaci pro dx=0 podobně jako v případě sčítání řady, ale tentokráte prokládáme přímku body I(dx) a I(2 dx), dostneme v případě polynomu odhad

(4 I(dx) - I(2 dx) ) / 3 = SpravnaHodnotaIntegralu + O(q^2) .

Konkrétně pro sinus (i když to není polynom...) prokládáme body
[1/8^2, 0.99678517188617] a [1/16^2, 0.99919668048507]
a získáme odhad 1.0000005166847065. V principu můžeme pokračovat a prokládat paraboly třemi body [1,I(dx)] ,[4, I(2 dx)] a [16, I(4 dx)], kubické křivky čyřmi body atd... Pro hezké funkce to vede na efektivní způsob výpočtu integrálu (Romberg). Také proto, že při výpočtu lze opakovaně použít již spočtené hodnoty.

Druhá varianta spočívá v tom, že místo abychom prováděli extrapolaci výsledných hodnot integrálu, extrapolujeme přímo ony formy. Protože Lagrangeův interpolační vzorec byl lineární ve funkčních hodnotách, můžeme extrapolovat i vektorovou funkci, neboť umíme sčítat a škálovat vektory. Pak prokládáme dva "body" [dx^2, dx.[1/2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/2] ] a [(2 dx)^2, (2 dx).[1/2,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1/2]] .

Ilustrace: Nalezli jsme extrapolovanou hodnotu lineární závislosti, jež prochází body [1,f1] a [4,f2] v bodě x=0 a to (4 I(dx) - I(2 dx) ) / 3. Použijme výsledek na f1= [...] a f2 = [...] viz výše:
(4 f1-f2)/3 = (1/3)*(4 dx.[1/2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/2] - 2 dx [1/2,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1/2]) =(dx/3).[1, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 1] ,
tedy (pro n dělitelné dvěmi)


což je známé Simpsonovo pravidlo. Podobně proložením paraboly "body"
[dx^2, dx.[1/2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/2] ] ,
[(2 dx)^2, (2 dx).[1/2,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1/2]] a
[(4 dx)^2, (4 dx).[1/2,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1/2]]
a extrapolací do dx=0 získáme formu
(2 dx / 45).[7, 32, 12, 32, 14, 32, 12, 32, 14, 32, 12, 32, 14, 32, 12, 32, 7].

Poznámka: Kromě toho jaké koeficienty mají být ve výše uvedených formách se můžeme také ptát, které z funkčních hodnot vzít, aby byl výsledek co nejpřesnější. Ukazuje se, že to není ekvidistantní vzorek funkčních hodnot. Zde je pro zajímavost uveden program používající osmibodový Gausův kvadraturní vzorec, kdybyste chtěli experimentovat...

program GaussInt;

type tRFunkce = function (x:real) : real;

function Gauss8(f: tRFunkce; a,b : real):real;
const t : array [1..4] of real = (0.183434642495649805,0.525532409916328986,0.796666477413626740,0.960289856497536232);
      w : array [1..4] of real = (0.181341891689180991,0.156853322938943644,0.111190517226687235,0.0506142681451881296);
var s,K,L : real;
    i : integer;
begin
  s := 0;
  K := (b+a)/2.0;
  L := (b-a)/2.0;
  for i:=low(t) to high(t) do s := s + ( f(K+L*t[i]) + f(K-L*t[i]) )*w[i];
  Gauss8:=s*(b-a);
end;

function fce(x:real):real;
begin
  fce := sin(x);
end;

function gce(x:real):real;
begin
  gce := sqrt(1-x*x);
end;

begin
  Writeln(Gauss8(fce, 0,Pi/2):1:16);
  Writeln(2*Gauss8(gce,-1,1 ):1:16);
  Readln;
end.

Jeho výstup, tedy čísla
1.00000000000000
3.14431044825165

jasně ukazují, že pro hezkou funkci, jakou je sinus, lze z osmi fukčních hodnot spočíst výsledek integrace přesně na 16 míst, zatímco pro funkce ošklivé máme potíže na čtvrtém desetinném místě.

Protože existuje mnoho různých způsobů, jímž může být funkce ošklivá, jedinou obecnou radou může být, že před integrací rodělíme v problematických bodech interval na více dílů a v každém podintervalu pak substitucí převedeme integrand na hezkou fuknci. (Samozřejmě, i tuto operaci můžeme více či méně převést na kvadraturní vzorce, pokud k tomu budmeme mít důvod.)

Náhodná čísla v počítači

Tady měla být zmínka o tom co to jsou (pseudo) náhodná čísla v počítači a jak je použít. Z časových důvodů jen uvedu že v Pascalu je k dispozici funkce
function random: real,
která vrátí (pseudo) náhodné číslo v [0,1) přičemž platí, že pravděpodobnost, že padne do daného podintervalu je rovna jeho délce.

.