Programování pro fyziky 2021/22 – Numerické recepty č. 4

Náš průvodní text je především v poznámkách k přednášce. Pomáhat si budeme i odkazy do knih Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR), které jsou (v Second Edition) na webu volně k dispozici. Vidět je ovšem lze jen z Adobe Readeru po instalaci pluginu FileOpen. Interní PDF prohlížeče Firefoxu a Chromu nám nepomohou.

Hledání kořenů

(NR kap. 9) Úlohou je řešit rovnici \(f(x)=0\) neboli hledat nulové body (kořeny) reálné funkce reálné proměnné. Při konečné přesnosti reálné aritmetiky jde spíše o řešení úlohy \(|f(x)|<\delta\). Pozor však na situace podle NR obr. 9.1.1: blízkými kořeny, nespojitostmi či singularitami funkce \(f(x)\) mohou být níže uvedené metody zmateny. Obvyklým požadavkem metod je lokalizace kořenu v intervalu \(\langle a,b\rangle\) nebo, v některých případech, odhad polohy kořenu \(x_0\). Kořeny se lokalizují analýzou funkce nebo jejím vykreslením (např. v gnuplotu) nebo postupy navrženými v NR kap. 9.1 (procedury zbrac pro expanzi intervalu a zbrak pro rozdělení intervalu). Některé metody lze zobecnit pro funkce komplexní proměnné nebo pro soustavy rovnic. Octave: fzero pro 1 rovnici, fsolve pro soustavu.

Metoda půlení intervalu (bisekce)

(NR kap 9.1) Není-li o kořenu k dispozici jiná informace než meze intervalu \(a_1\) a \(b_1\), pro které ovšem platí \(f(a_1)f(b_1)<0\), omezujeme kořen posloupností intervalů \(\langle a_1,b_1\rangle \supset \langle a_2,b_2\rangle \supset \dots \supset \langle a_k,b_k\rangle \supset \dots\) tak, že v každém kroku aktuální interval rozpůlíme a vybereme tu polovinu, pro kterou zůstává zachováno \(f(a_k)f(b_k)<0\). Půlení zastavíme při dosažení \(|f(x)|<\delta\) nebo \(|b_k-a_k|<\varepsilon\). Chybu v \(k\)-tém kroku lze ztotožnit s aktuální délkou intervalu, \(\varepsilon_k=b_k-a_k\). Zjevně platí \(\varepsilon_{k+1}=\varepsilon_k/2\), chyba metody se vyvíjí lineárně. Každé půlení dodá mantise 1 bit přesnosti, není tedy třeba více než 24, resp. 53 půlení při reálné přesnosti 4, resp. 8 bytů. Není-li funkce spojitá, metoda může lokalizovat místo kořenu nespojitost či singularitu. Metoda nemůže lokalizovaný kořen ztratit, je vždy konvergující.

format long; f=@(x) x^3-1; a=.5; b=2; eps=1e-7;
do x=(a+b)*.5, if f(x)*f(a)<0, b=x; else a=x; end until b-a<eps

Metoda regula falsi

(NR kap 9.2) Proložíme sečnou body omezující kořen, \((a,f(a)),\ (b,f(b)),\ f(a)f(b)<0\), vypočteme její kořen a ztotožníme jej s jednou z předchozích mezí \(a,b\) tak, aby spolu s druhou mezí i nadále platilo \(f(a)f(b)<0\) (NR obr. 9.2.2). Jde proto opět o vždy konvergující metodu. Rovnici sečny můžeme zapsat ve tvaru

\[y(x)=f(a)+(f(b)-f(a))(x-a)/(b-a),\]

odkud pro její kořen \(x\), v němž \(y(x)=0\), plyne

\[x=(af(b)-bf(a))/(f(b)-f(a)).\]
format long; f=@(x) x.^3-1; a=.5, b=2, eps=1e-7; x=a; fa=f(a); fb=f(b);
do xo=x; x=(a*fb-b*fa)/(fb-fa), fx=f(x); if fx*fa<0, b=x; fb=fx; else a=x; fa=fx; end until abs(x-xo)<eps

% včetně kreslení...
format long; f=@(x) x.^3-1; a=.5, b=2, eps=1e-7; x=a; fa=f(a); fb=f(b);
hold off; z=a:(b-a)/100:b; plot(z,f(z),'g',z,zeros(size(z)),'k'); xlim([a b]);
do
  hold on;
  plot([a b],[fa fb],'r',[a b],[fa fb],'r.');
  xo=x; x=(a*fb-b*fa)/(fb-fa), fx=f(x); if fx*fa<0, b=x; fb=fx; else a=x; fa=fx; end
until abs(x-xo)<eps

Metoda sečen

(NR kap 9.2) Proložíme sečnou body omezující kořen, \((x_1,f_1)),\ (x_2,f_2)\), kde \(f_k=f(x_k)\), vypočteme její kořen a označíme jej \(x_3\). Pokračujeme dál spojováním naposledy získaných dvojic bodů sečnami až k dosažení konvergence nebo divergence. Není zajištěno, aby každý interval \(\langle x_k,x_{k+1}\rangle\) omezoval kořen \(f(x)\), metoda proto není vždy konvergující (NR obr. 9.2.1). Může však konvergovat rychleji než předchozí dvě metody. Rovnice sečny je zde

\[y(x)=f_{k-1}+(f_k-f_{k-1})(x-x_{k-1})/(x_k-x_{k-1}),\]

odkud pro její kořen \(x_{k+1}\), v němž \(y(x_{k+1})=0\), plyne

\[x_{k+1}=(x_{k-1}f_k-x_kf_{k-1})/(f_k-f_{k-1}).\]
format long; f=@(x) x.^3-1; a=.5, b=2, eps=1e-7; x1=a; f1=f(a); x2=b; f2=f(b);
do x=(x1*f2-x2*f1)/(f2-f1), x1=x2; f1=f2; x2=x; f2=f(x); until abs(x2-x1)<eps

% včetně kreslení...
format long; f=@(x) x.^3-1; a=.5, b=2, eps=1e-7; x1=a; f1=f(a); x2=b; f2=f(b);
hold off; z=a:(b-a)/100:b; plot(z,f(z),'g',z,zeros(size(z)),'k'); xlim([a b]);
do
  hold on; plot([x1 x2],[f1 f2],'r',[x1 x2],[f1 f2],'r.')
  x=(x1*f2-x2*f1)/(f2-f1), x1=x2; f1=f2; x2=x; f2=f(x);
until abs(x2-x1)<eps

Newtonova metoda tečen

(NR kap 9.4) Funkci \(f(x)\) linearizujeme neboli aproximujeme v okolí odhadu polohy kořenu \(x_0\) nultým a prvním členem Taylorova rozvoje neboli vedeme tečnu k \(f(x)\) bodem \((x_0, f_0)\) a nalezneme kořen této aproximace. Ten označíme \(x_1\), bodem \((x_1, f_1)\) vedeme další tečnu, nalezneme její kořen a pokračujeme takto až k dosažení konvergence nebo divergence. Rovnice pro kořen tečny vedené bodem \(x_k\) je tedy

\[f(x)\approx f(x_k)+f'(x_k)(x-x_k)=0,\]

odkud plyne vzorec Newtonovy metody,

\[x_{k+1}=x_k-f(x_k)/f'(x_k).\]

Pro chybu \(\varepsilon_k=x_k-x\), kde \(x\) je hledaný kořen, zřejmě platí stejný vztah, \(\varepsilon_{k+1}=\varepsilon_k-f(x_k)/f'(x_k)\). Dosadíme-li sem Taylorovy rozvoje pro \(f(x_k)\) a \(f'(x_k)\) kolem kořenu \(x\), kde ovšem \(f(x)=0\), tedy \(f(x_k)\approx \varepsilon_kf'(x)+\frac12\varepsilon_k^2f''(x)\) a \(f'(x_k)\approx f'(x)+\varepsilon_k f''(x)\), dostaneme odhad

\[\varepsilon_{k+1}\approx \varepsilon_k-\frac{\varepsilon_kf'(x)+\frac12\varepsilon_k^2 f''(x)}{f'(x)+\varepsilon_k f''(x)} \approx \varepsilon_k^2 \frac{f''(x)}{2f'(x)}\,.\]

Chyba metody se vyvíjí kvadraticky, pro \(\varepsilon_k=10^{-2}\) může pokračovat přes \(10^{-4}\) a \(10^{-8}\) k \(10^{-16}\), tedy dosáhnout plné 15místné přesnosti 8bytového datového typu jen několika kroky. Je-li lineární aproximace (tečna) nevhodná, metoda selže. (Kořen tečny bude mimo oblast.) Praktickým postupem je užívání rychle konvergující Newtonovy metody, od které se při jejím pokusu o opuštění omezujícího intervalu ustoupí ke vždy konvergující bisekci. Aproximujeme-li ve vzorci Newtonovy metody \(f'(x_k)\) numericky, \(f'(x_k)=[f(x_k)-f(x_{k-1})]/(x_k-x_{k-1})\), získáme vzorec metody sečen.

format long; f=@(x) x.^3-1; df=@(x) 3*x.^2; a=.5; b=2; eps=1e-7; x=(a+b)*.5;
x, do xold=x; x=xold-f(x)/df(x) until abs(x-xold)<eps

% včetně kreslení...
format long; f=@(x) x.^3-1; df=@(x) 3*x.^2; a=.5; b=2; eps=1e-7; x=(a+b)*.5;
hold off; z=a:(b-a)/100:b; plot(z,f(z),'g',z,zeros(size(z)),'k'); xlim([a b]);
x
do
  hold on; plot(z,f(x)+df(x)*(z-x),'r',[x],[0],'k.',[x],[f(x)],'r.')
  xold=x; x=xold-f(x)/df(x)
until abs(x-xold)<eps

Newtonovou metodou (i metodou sečen) lze hledat také kořeny funkcí komplexní proměnné. Graficky vděčnou úlohou je vykreslit ty body komplexní roviny, z nichž Newtonova metoda pro funkci \(f(z)=z^3-1\) konverguje ke zvolenému kořenu, např. \(z=1\). (NR obr. 9.4.4 sice v PDF souboru chybí, ale vykreslí jej následující skript.)

ndot=250000; nstep=50; eps=1e-5;
f=@(z) z.^3-1;
df=@(z) 3*z.^2; root=1;
x1=-2; x2=2; y1=-2; y2=2;
z0=x1+(x2-x1)*rand(ndot,1)+i*(y1+(y2-y1)*rand(ndot,1));
z=z0;
for n=1:nstep z=z-f(z)./df(z); end; result=z0(abs(z-root)<eps);
disp(['Plotting ',num2str(prod(size(result))),' dots']);
plot(result,'.','markersize',2)

Newtonova metoda pro soustavu rovnic

(NR kap. 9.6) Uvažujme soustavu \(N\) nelineárních rovnic

\[\mathbf{F}(\mathbf{x})=\mathbf{0} \quad\mbox{neboli}\quad F_i(x_1,x_2,\dots,x_N)=0,\ i=1,\dots,N.\]

Linearizujeme Taylorovým polynomem,

\[\mathbf{F}(\mathbf{x}+\delta\mathbf{x})\approx \mathbf{F}(\mathbf{x})+\frac{\partial\mathbf{F}(\mathbf{x})}{\partial\mathbf{x}}\cdot\delta\mathbf{x} \quad\mbox{neboli}\quad F_i(\mathbf{x}+\delta\mathbf{x})\approx F_i(\mathbf{x})+\sum_{j=1}^N\frac{\partial F_i(\mathbf{x})}{\partial x_j}\delta x_j,\]

kde budeme značit Jacobiovu matici \(\mathbf{J}(\mathbf{x})=\partial\mathbf{F}/\partial\mathbf{x}\), \(J_{ij}=\partial F_i/\partial x_j\). Podobně jako v Newtonově metodě pro jednu rovnici hledáme \(\delta\mathbf{x}\) tak, aby \(\mathbf{F}(\mathbf{x}+\delta\mathbf{x})=\mathbf{0},\) tedy řešíme v každém kroku soustavu lineárních algebraických rovnic

\[\mathbf{J}\cdot\delta\mathbf{x}=-\mathbf{F}.\]

Vzorec Newtonovy metody pro soustavu rovnic má tak tvar

\[\mathbf{x}_{k+1}=\mathbf{x}_k-\mathbf{J}^{-1}(\mathbf{x}_k)\cdot\mathbf{F}(\mathbf{x}_k)\,.\]

Počáteční hodnotu \(\mathbf{x}_0\) musíme odhadnout.

Řešení soustavy nelineárních rovnic \(\begin{array}{rcl}2x_1x_2+3x_1x_3-\ \,93 &=& 0\\ 4x_1x_2+5x_2x_3-235 &=& 0\\ 6x_1x_3+7x_2x_3-371 &=& 0 \end{array}\)

kmax=7; x=[1; 2; 3];
printf("%2d %15.10f %15.10f %15.10f\n",0,x)
for k=1:kmax
  F=[2*x(1)*x(2)+3*x(1)*x(3)-93; 4*x(1)*x(2)+5*x(2)*x(3)-235; 6*x(1)*x(3)+7*x(2)*x(3)-371];
  J=[2*x(2)+3*x(3) 2*x(1) 3*x(1); 4*x(2) 4*x(1)+5*x(3) 5*x(2); 6*x(3) 7*x(3) 6*x(1)+7*x(2)];
  x=x-J\F;
  printf("%2d %15.10f %15.10f %15.10f\n",k,x)
end

Octave: Pro hledání nulových bodů jedné funkce je určen řešič fzero(f,x0), kde f je skalární funkce skalárního argumentu a x0 odhad polohy nulového bodu (skalár nebo dvouprvkové pole mezí).

format long; f=@(x) sin(x); x0=3; fzero(f,x0), x0=[3 4]; fzero(f,x0)  % vrátí pí

Pro soustavy slouží řešič fsolve(F,x0,options), kde F je vektorová funkce a x0 vektor s počátečním odhadem.

F=@(x) [2*x(1)*x(2)+3*x(1)*x(3)-93; 4*x(1)*x(2)+5*x(2)*x(3)-235; 6*x(1)*x(3)+7*x(2)*x(3)-371];
x0=[1; 2; 3]; fsolve(F,x0)                                            % vrátí totéž jako skript výše

V MATLABu je fsolve součástí Optimization toolboxu.