Přechod na menu, Přechod na obsah, Přechod na patičku
Početní praktikum Ústav teoretické fyziky a astrofyziky, Přírodovědecká fakulta Masarykovy univerzityMgr. Ing. arch. Petr Kurfürst, Ph.D.

C.3 Numerické metody výpočtů funkcí jedné proměnné

C.3.1 Hledání kořene funkce jedné proměnné – Newtonova metoda

Kořeny obecně nelineární funkce (rovnice) $f(x)=0$ často nelze vyjádřit formou explicitního analytického vzorce. K nalezení řešení takové rovnice musíme potom použít některou z numerických (iteračních) metod, kdy pomocí určitého počtu počátečních aproximací hledaného kořene $x_0$ generujeme posloupnost $x_1,\,x_2,\,x_3,\,\ldots\,,$ která ke kořenu $x_0$ konverguje. V některých případech je třeba zadat interval $a,\,b$, který podle předběžného předpokladu obsahuje hledaný kořen, čím lépe se k němu na počátku přiblížíme, tím rychleji daná metoda konverguje. V následujících příkladech předpokládejme reálnou spojitou funkci $f(x)$ s odpovídajícím počtem spojitých derivací na vymezeném intervalu, s hledaným kořenem $f(x_0)=0$.

Počáteční odhad intervalu (intervalů) kde se kořen (kořeny) mohou nalézat provedeme například grafickou metodou: pomocí vhodného výpočetního programu nebo vypisováním funkčních hodnot do tabulky vykreslíme funkci $f(x)$ a vyhledáme její přibližné průsečíky s osou $x$. Například u funkce, dané předpisem

$$ f(x)=x^3–3x^2+2x-3 $$
C.17

snadno zjistíme že existuje jeden reálný kořen, který musí s určitostí ležet uvnitř intervalu $x_{0}\in (2,3)$.

Existuje celá řada možných numerických postupů (například metoda sečen, atd.), asi nejznámější je tzv. Newtonova metoda neboli metoda tečen. Vyjdeme z počáteční aproximace $x_0$ a postupně počítáme $x_1,\,x_2,\,x_3,\,\ldots\,.$ Známe-li určitou aproximaci $x_k$ a chceme určit lepší aproximaci $x_{k+1}$, proložíme bodem $[x_k,f(x_k)]$ tečnu ke křivce $y=f(x)$, průsečík této tečny s osou $x$ považujeme potom za hodnotu $x_{k+1}$. Dostáváme tak rovnici popsané tečny ve tvaru

$$ f^\prime(x_k )=\{3x_k^2–6x_k+2\}=\frac{f(x_k)}{x_k-x_{k+1}}, $$
C.18

z níž odvodíme vztah pro výpočet každého následujícího kroku (iterace),

$$ x_{k+1}=x_k-\frac{f(x_k)}{f^\prime(x_k )}. $$
C.19

Ukončení výpočtu s tím, že hodnota poslední iterace $x_{k+1}$ je prohlášena za hledaný kořen $x_0$ s požadovanou přesností, nastane (podle velikosti malého čísla $\epsilon$ které stanovuje požadovanou přesnost) například v těchto případech:

$ |x_{k+1}-x_k|\leq\epsilon$ nebo $|f(x_{k+1})|\leq\epsilon$.
C.20
Příklad možného způsobu naprogramování rovnice C.17 – program fortran95:
program newton                          !deklarace názvu programu 	
implicit none 	
 	
integer :: i                            !deklarace celočíselných proměnných v záhlaví 
                                        !programu: i = pořadové číslo prostorového kroku 	
double precision :: x, dx, f, df 	!deklarace reálných veličin, kde výrazy x $=x_k$, 
                                        !dx $=x_k-x_{k+1}, f=f(x_k)$, df $=f^\prime(x_k)$
                                        !s dvojitou přesností
 	
x=3.d0 	                                !odhad vstupní hodnoty $x_k$ 	
do                                      !výpočetní cyklus 	
  i = i+1 	
  f = x**3 – 3.d0*x**2 + 2.d0*x – 3.d0 	!vlastní rovnice C.17 	
  df = 3.d0*x**2 – 6.d0*x + 2.d0        !derivace funkce C.17 	
  dx = f/df                             !rovnice C.19 	
  x = x-dx                              !nová hodnota $x_k$ 	
  if (dabs(dx).lt.1.d-12) exit          !stop kritérium: $|x_k-x_{k+1}|<10^{-12}$
end do 	

write (100,*) x                         !zápis kořene funkce f do souboru fort.100

stop                                    !zastavení celého procesu 	

end program newton                      !konec programu
Tabulka výsledků programu výpočtu kořene funkce $f(x_k)$ pro jednotlivé iterace $k$
$k$$x_k$$f(x_k)$$x_k-x_{k+1}$
03.00000000000000003.00000000000000000.27272727272727271
12.72727272727272750.425995492111198365.3581553581554045E-002
22.67369117369117331.4723079585858834E-0021.9886039436713345E-003
32.67170256974750191.9848200396133109E-0052.6880854327577796E-006
42.67169988166206903.6242120415863610E-0114.9083680000275664E-012
52.67169988165716041.7763568394002505E-0152.4057679206275180E-016
nahoru

C.3.2 Numerické derivování

Nejjednodušší numerická aproximace 1. derivace má podobu tzv. dopředné diference,

$$ f^\prime(x)\approx\frac{f(x+h)−f(x)}{h}, $$
C.21

kde $h=\Delta x>0$, s chybou aproximace $\delta f^\prime(x)$ vyjádřenou pomocí Taylorova rozvoje $x$, tj. $\delta f^\prime(x)=-({h}/{2})f^{\prime\prime}(\xi)$, kde $\xi\in(x,x+h)$. Podobně jednoduchá je i tzv. zpětná diference,

$$ f^\prime(x)\approx\frac{f(x)−f(x-h)}{h}, $$
C.22

s chybou stejného řádu. Aproximací s vyšší přesností je tzv. centrální diference,

$$ f^\prime(x)\approx\frac{f(x+h)−f(x-h)}{2h}, $$
C.23

s chybou aproximace $\delta f^{\prime}(x)=-({h^2}/{6})f^{\prime\prime\prime}(\xi)$, kde $\xi\in(x-h,x+h)$, zabírající ovšem dva prostorové kroky (buňky) výpočetní sítě. Analogickým způsobem můžeme odvodit i 2. derivaci ve tvaru

$$ f^{\prime\prime}(x)\approx\frac{f(x+h)−2f(x)+f(x-h)}{h^2}, $$
C.24

s chybou aproximace $\delta f^{\prime\prime}(x)=-({h^2}/{12})f^{(4)}(\xi)$, kde $\xi\in(x-h,x+h)$.

Existují i přesnější a propracovanější diferenční schémata (viz například van Leer, 1977, 1982; Vitásek, 1987; LeVeque, 2002), například jednostranná aproximace 1. derivace, která je 2. řádu přesnosti,

$$ f^\prime(x)\approx\frac{−3f(x)+4f(x+h)−f(x+2h)}{2h}, $$
C.25

s chybou aproximace $\delta f^{\prime}(x)=({h^2}/{3})f^{\prime\prime\prime}(\xi)$, kde $\xi\in(x,x+2h)$, nebo aproximace 1. derivace, která je 4. řádu přesnosti, ve tvaru

$$ f^\prime(x)\approx\frac{−f(x+2h)+8f(x+h)−8f(x-h)+f(x-2h)}{12h}, $$
C.26

tedy s chybou, která je řádu $h^4$, atd. Jejich nevýhodou ovšem je, že zabírají několik prostorových intervalů (buněk) výpočetní sítě a také při složitých a objemných výpočtech díky nim mohou narůstat nároky na dobu výpočtu. Je proto vždy nutné zvážit výpočetní schéma, adekvátní dané úloze a její požadované přesnosti, odpovídající ale reálným možnostem používaného výpočetního zařízení.

nahoru

C.3.3 Numerické integrování

je vždy založené na nahrazení složitě ohraničeného geometrického útvaru (plochy pod křivkou dané funkce v případě jedné proměnné) jednodušším útvarem, nebo součtem takových útvarů. Používá se také název numerická kvadratura, ve smyslu konstrukce plošných (tedy dvourozměrných, kvadraturních) útvarů. Ukážeme zde příklady pouze nejběžnějších (většinou ovšem zcela dostačujících) způsobů numerické integrace funkce jedné proměnné pomocí tzv. Newton-Cotesových vzorců, existuje samozřejmě celá řada jiných metod numerické integrace, například Gaussovy kvadraturní vzorce, Rombergova kvadratura, atd. Nebudeme zde uvádět ani přesnosti a způsob stanovení chyb, atd., vše je standardně dostupné v literatuře.

nahoru

Newton-Cotesovy (kvadraturní) vzorce

  • Obdélníková metoda

    Tato metoda se formálně nepočítá mezi tzv. Newton-Cotesovy vzorce, představuje sice nejjednodušší ale zároveň nejméně přesnou numerickou integrační metodu, kdy se určitý integrál dané funkce (tj. velikost plochy pod křivkou grafu funkčních hodnot funkce $f(x)$ v rámci intervalu $\langle a,b\rangle$) aproximuje obdélníkem. Tuto aproximaci můžeme zpřesnit, rozdělíme-li například interval $\langle a,b\rangle$ na zvolený počet $n$ dílčích stejných intervalů, vypočítáme obdélníkovou aproximaci pro každý interval zvlášť a výsledky sečteme, tedy

    $$ \int_a^bf(x)\,\text{d} x\approx\frac{b-a}{n}\sum_{k=0}^{n-1}f\left(a+k\,\frac{b-a}{n}\right).$$
    C.27

  • Lichoběžníková metoda

    představuje přesnější numerickou integrační metodu, kdy se určitý integrál dané funkce aproximuje lichoběžníky (body funkce se spojí úsečkami). Rozdělíme-li interval $\langle a,b\rangle$ na zvolený počet $n$ dílčích stejných intervalů, vypočítáme lichoběžníkovou aproximaci opět pro každý interval zvlášť a výsledky sečteme, tedy

    $$ \int_a^bf(x)\,\text{d} x\approx\frac{b-a}{2n}\sum_{k=0}^{n-1}\left[f\!\left(x_{k+1}\right)+f\!\left(x_k\right)\right], $$
    C.28

    kde $x_k=a+k(b-a)/n$.

  • Simpsonovo pravidlo

    založené na kvadratické (parabolické) interpolaci dílčích intervalů integrované funkce.

    V případě integrace polynomů dává tato metoda velmi přesné výsledky. Složená aproximace Simpsonovým pravidlem, kdy interval $\langle a,b\rangle$ je rozdělen na sudý počet $n$ dílčích intervalů, má tvar (kdy $x_k=a+k(b-a)/n$)

    $$ \int_a^bf(x)\,\text{d} x\approx\frac{b-a}{3n}\sum_{k=1}^{n/2}\left[f\!\left(x_{2k-2}\right)+4f\!\left(x_{2k-1}\right)+f\!\left(x_{2k}\right)\right]. $$
    C.29

  • Simpsonovo $3/8$ pravidlo

    (nebo také druhé Simpsonovo pravidlo): založené na kubické interpolaci dílčích intervalů integrované funkce.

    Složená aproximace Simpsonovým 3/8 pravidlem, kdy interval $\langle a,b\rangle$ je rozdělen na $n$ dílčích intervalů, kde $n$ je dělitelné třemi, má tvar (kdy opět $x_k=a+k(b-a)/n$)

    $$ \int_a^bf(x)\,\text{d} x\approx\frac{3(b-a)}{8n}\sum_{k=1}^{n/3}\left[f\!\left(x_{3k-3}\right)+3f\!\left(x_{3k-2}\right)+3f\!\left(x_{3k-1}\right)+f\!\left(x_{3k}\right)\right]. $$
    C.30

    Obdobným způsobem lze sestrojit Newton-Cotesovy formule libovolně vyšších řádů, vyšší řády než 4 (Booleovo pravidlo) se nicméně používají málo, jejich nevýhodou je velmi rychle (až exponenciálně) narůstající chyba integrace.

nahoru

C.3.4 Jednoduché numerické metody řešení obyčejných diferenciálních rovnic

Existuje opět celá řada způsobů numerického řešení obyčejných diferenciálních rovnic, například řešení rovnic Eulerovou metodou nebo tzv. Runge-Kuttovou metodou (metodami), atd., viz například (Humlíček, 2009). Ukážeme zde pouze často používanou jednoduchou tzv. metodu střelby a také příklad řešení obyčejné diferenciální rovnice 2. řádu pomocí tridiagonální matice (viz odstavec C.2.1):

nahoru

Metoda střelby

je jednoduchá metoda, založená na vystřelení dané funkce v závislosti na okrajových podmínkách z pevného bodu daným směrem. Hledáme potom takové koeficienty funkce, které zajistí dopadnutí funkce do požadovaného bodu. Metodu si ukážeme na příkladu řešení tzv. Lane-Emdenovy rovnice, což je obyčejná diferenciální rovnice 2. řádu, obvykle zapsaná v implicitním tvaru

$ \dfrac{1}{x^2}\dfrac{\text{d}}{\text{d} x}\left({x^2 \dfrac{\text{d} y}{\text{d} x}}\right)+y^n=0$, a tedy $y^{\prime\prime}+\dfrac{2}{x}y^{\prime}+y^n=0$,
C.31

kde $x$ je nezávisle proměnná, $y$ je závisle proměnná a $n$ je konstanta. Tato rovnice je řešitelná analyticky pouze pro $n=0,\,1,\,5$, pro všechna ostatní $n$ musí být řešena numericky. Pro $n=0$ získáme řešení přímou integrací se zahrnutím uvedených okrajových podmínek, pro $n=1$ řešíme sférickou Besselovu diferenciální rovnici (viz rovnice B.117), Pro $n=5$ dostáváme řešení prostřednictvím tzv. Emdenovy transformace, kde $y=Ar^\omega s$, kde $r,s$ jsou nové proměnné a $\omega=2/(n-1)$.

Budeme počítat rovnici C.31 pro $n=1,\!5$ a s Newtonovými okrajovými podmínkami $y(0)=1,\,y^\prime(0)=0$. Algoritmus je tedy vystřelen z bodu [0,1] vodorovně, nová hodnota $y$ je v každém prostorovém kroku vypočítána jako $y=y+(\Delta y/\Delta x)\Delta x$. Druhá derivace $y^{\prime\prime}$ je rozepsána jako $(y^{\prime})^\prime$, tedy $(\Delta y/\Delta x)/\Delta x$ a každá nová hodnota $y^{\prime\prime}$ je v každém prostorovém kroku počítána jako $(\Delta y/\Delta x)/\Delta x=(\Delta y/\Delta x)/\Delta x-[(2/x)(\Delta y/\Delta x)+y^n]$, což můžeme po vynásobení celé rovnice $\Delta x$ přepsat jako $(\Delta y/\Delta x)=(\Delta y/\Delta x)-[(2/x)(\Delta y/\Delta x)+y^n]\Delta x$. Numerický algoritmus rovnice C.31 lze tedy zapsat například takto (fortran 95): 

program Emden                         !deklarace názvu programu 	
implicit none 	
double precision :: x, y, dydx, n, dx !deklarace reálných veličin x, y, dydx, n, dx 
                                      !s dvojitou přesností, kde dydx = $y^\prime = \Delta y/\Delta x$ 
                                      !a dx = $\Delta x$ 	
x=0.d0                                !deklarace parametrů a okrajových podmínek 	
y=1.d0 	
dydx=0.d0 	
dx=1.d-3 	
n=1.5d0 	
 	
do                                    !výpočetní cyklus 	
  x=x+dx 	
  y=y+dydx*dx 	 	
  dydx=dydx-(2.d0*dydx/x+y**n)*dx 	
  if(x$>$15.d0) exit                    !stop kritérium, dané předběžným odhadem 	
  write(1,*) x,y,dydx                 !zápis do souboru fort.1 	
end do                                !konec cyklu 	
end program Emden
nahoru

Příklad řešení obyčejné diferenciální rovnice 2. řádu pomocí tridiagonální matice

Řešení okrajové úlohy obyčejné diferenciální rovnice 2. řádu $p(x)y^{\prime\prime}(x)+q(x)y^\prime(x)+r(x)y(x)=f(x)$, na intervalu $a,b$, s Dirichletovými okrajovými podmínkami $y(a)=A$, $y(b)=B$, lze snadno řešit pomocí tridiagonální matice (viz odstavec 3.2.1). Jako příklad řešení uvedeme jednoduchou rovnici s konstantními koeficienty $y^{\prime\prime}+3y^\prime+2y=(20x+29)\,\text{e}^{3x}$ (snadno řešitelnou i analyticky), na intervalu $\langle 0,1\rangle$, s Dirichletovými okrajovými podmínkami $y(0)=0$, $y(1)=1$.

Po přepsání dané rovnice do diferenčního schématu s $n+2$ body prostorové sítě, kdy jednotlivé derivace levé strany rozepíšeme podle rovnic C.23 a C.24, dostáváme rovnici ve tvaru

$$ \left(p-\frac{h}{2}q\right)y_{i-1} + \left(h^2r-2p\right)y_i+\left(p+\frac{h}{2}q\right)y_{i+1}=h^2f_i, $$
C.32

kde $i=0,\,1,\,\ldots\,,n,\,n+1$, s koeficienty $p=1,\,q=3,\,r=2$ a s konstantním prostorovým krokem $h=(x_{n+1}-x_0)/(n+1)=[x(1)-x(0)]/(n+1)$, při $n=99$ tak bude $h=0,01$. Označíme-li jednotlivé závorky na levé straně rovnice C.32 postupně $P,\,Q,\,R$, dostáváme rovnici s tridiagonální maticí na levé straně ve tvaru

$$ \begin{pmatrix} Q & R\\[4pt] P & Q & R\\[4pt] & & \ddots & \ddots & \ddots \\[4pt] & & & P & Q & R\\[4pt] & & & & P & Q \end{pmatrix} \begin{pmatrix} y_1\\[4pt] y_2\\[4pt] \vdots\\[4pt] y_{n-1}\\[4pt] y_n \end{pmatrix}= \begin{pmatrix} h^2f_1-Py_0\\[4pt] h^2f_2\\[4pt] \vdots\\[4pt] h^2f_{n-1}\\[4pt] h^2f_n-Ry_{n+1} \end{pmatrix}. $$
C.33

Numerický algoritmus zapíšeme následujícím způsobem (fortran 95): 

program tridiag                       !deklarace názvu programu 	
implicit none 
integer :: i, ni, N, LDB, NRHS, INFO  !deklarace celočíselných proměnných: 
                                      !i = pořadové číslo nezávisle proměnné,
                                      !ni = celkový počet diskrétních hodnot,
                                      !ostatní jsou parametry podprogramu DGTSV
                                      !balíčku LAPACK – viz odstavec C.1 	
parameter(ni=99,N=ni,LDB=N,NRHS=1)    !hodnoty celočíselných proměnných 	
double precision :: x(ni) ,y(ni), h(ni), p(ni), q(ni), r(ni), f(ni), B(LDB,NRHS), 
DL(N-1), DU(N-1), D(N)                !deklarace reálných proměnných jako pole
                                      !s dvojitou přesností 	
double precision, parameter :: x0=0.d0, xn=1.d0, y0=0.d0, yn=1.d0 	
                                      !deklarace reálných konstant s dvojitou přesností
do i=1,N                              !výpočetní cyklus prostorového kroku 	
  x(i)=x0+(xn-x0)*i/dfloat(ni+1)  
                                      !příkaz dfloat mění celočíselnou proměnnou na reálnou 	
end do 	

do i=1,N                              !hlavní výpočetní cyklus 	
  h(i)=x(i)-x(i-1)                    !konstantní krok $h$ je zde zapsán obecně 	
  p(i)=1.d0                           !konstantní koeficienty jsou zapsány obecně 	
  q(i)=3.d0 	
  r(i)=2.d0 	
  f(i)=(20.d0*x(i)+29.d0)*dexp(3.d0*x(i)) 	

  if(i.eq.1) then                     !spodní okrajová podmínka 	
    B(i,1)=h(i)**2.d0*f(i)-(p(i)-q(i)*h(i)/2.d0)*y0 	
  elseif((i.gt.1).and.(i.lt.N)) then  !hlavní pole 	
    B(i,1)=h(i)**2.d0*f(i) 	
  else 	                              !horní okrajová podmínka 	
    B(i,1)=h(i)**2.d0*f(i)-(p(i)+q(i)*h(i)/2.d0)*yn 	
  endif 	
end do 	

do i=1,N-1                            !zadání spodní diagonály 	
  DL(i)=p(i)-q(i)*h(i)/2.d0 	
end do 	

do i=1,N                              !zadání hlavní diagonály 	
  D(i)=r(i)*h(i)**2.d0-2.d0*p(i) 	
end do 	

do i=2,N                              !zadání horní diagonály 	
  DU(i)=p(i)+q(i)*h(i)/2.d0 	
end do 	
 	
!volání podprogramu DGTSV balíčku LAPACK (viz oddíl C.1): 	
call DGTSV(N, NRHS, DL, D, DU, B, LDB, INFO) 	
if(INFO.ne.0) write(*,*) ``INFO='',INFO,``!!!'' 	
 	
write(1,*) x0, y0 	
do i=1,N 	
  y(i)=B(i,1) 	
  write(1,*) x(i), y(i)               !zápis do souboru fort.1 	
end do 	
write(1,*) xn, yn 	
 	
end program tridiag