Numerické výpočty

Antonín Slavík
MFF UK

Machine-precision / arbitrary-precision numbers

Machine precision

= desetinná čísla tak, jak jsou implementována v procesoru počítače

Vysoká rychlost

Omezený počet cifer

Chybí informace o přesnosti výsledku

Závislost na hardwaru

Arbitrary precision

= softwarová implementace operací s desetinnými čísly

Libovolný počet cifer

Informace o přesnosti výsledku

Nezávislost na hardwaru

Nižší rychlost

Příklady

MachineNumberQ[x] testuje, zda x je číslo v machine precision.

N[výraz] vrací výsledek v machine precision.

x=N[Sqrt[2]]

MachineNumberQ[x]

Mathematica nevypisuje všechny cifry. Lze je získat pomocí FullForm:

FullForm[x]

Clear[x]

Desetinné číslo je obvykle chápáno jako číslo v machine precision.

MachineNumberQ[3.141592653]

Pouze v případě, že kapacita machine precision nestačí na uložení daného desetinného čísla, se použije arbitrary precision.

MachineNumberQ[3.1415926535897932385]

N[výraz, počet platných cifer] vrací výsledek v arbitrary precision (i v případě, že by se výsledek vešel do machine precision).

N[Sqrt[2],20]

MachineNumberQ[%]

Čísla v arbitrary precision lze také zadávat ve tvaru x`p, kde p je požadovaná přesnost.

20.0`10

MachineNumberQ[%]

Technické detaily

Kolik platných cifer může mít číslo v machine precision?

?$MachinePrecision

$MachinePrecision

Reprezentace desetinných čísel na typickém PC (double precision, 64 bitů):

Znaménko: 1 bit
Exponent: 11 bitů
Zlomková část: 53 bitů (explicitně uloženo 52)

Reprezentované číslo: (((-1))^(znaménko))(1.(Subscript[b, 51])...(Subscript[b, 0]))(2^(e-1023))

Přípustné hodnoty exponentu jsou e∈{1,...,2046}; pro e=0 je číslo chápáno jako nula, pro e=2047 jako nekonečno.

(Viz např. Wikipedia, Double precision floating-point format a Floating point)

53 platných cifer ve dvojkové soustavě:

Log[10.0,2^53]

Minimální a maximální čísla, se kterými lze počítat v machine precision (v absolutní hodnotě):

{$MinMachineNumber,$MaxMachineNumber}

Log[2,%]

Hodnoty bitů uložených v paměti lze získat pomocí RealDigits.

?RealDigits

{d,e}=RealDigits[7.,2]

Length[d]

RealDigits[7.*4,2]

(10^23) ve dvojkové soustavě:

IntegerDigits[10^23,2]

Length[%]

V machine precision nelze (10^23) uložit zcela přesně (příliš mnoho cifer) :

{d,e}=RealDigits[10.^23,2]

FromDigits[{d,e},2]

Ekvivalentně pomocí SetPrecision[x,Infinity]:

SetPrecision[10.^23,Infinity]

0.1 v machine precision (ve dvojkové soustavě má nekonečný rozvoj):

RealDigits[.1,2]

SetPrecision[0.1,Infinity]

N[%,30]

?$MachineEpsilon

$MachineEpsilon

Log[2,$MachineEpsilon]

Tj. rozdíl mezi dvěma sousedními machine-precision čísly z intervalu [1,2) je (2^(-52)), rozdíl mezi dvěma sousedními machine-precision čísly z intervalu [2,4) je (2^(-51)) atd.

Závislost na pořadí operací:

10.^20+1-10.^20

10.^20-10.^20+1

Vysvětlení:

IntegerDigits[10^20,2]

RealDigits[10.^20,2]

Zcela chybný výsledek:

N[Sin[(10^25)]]

N[Sin[(10^25)],20]

Vysvětlení:

Sin[N[(10^25)]]

SetPrecision[N[10^25],Infinity]

Minimální a maximální čísla, se kterými lze počítat v arbitrary precision (v absolutní hodnotě):

{$MinNumber,$MaxNumber}

2$MaxNumber

Kontrola chyb

Přesnost u arbitrary-precision numbers

Číslo x reprezentované v arbitrary precision je vždy doplněno údajem 
o své přesnosti p, tj. vnitřní reprezentace takového čísla je x`p.

N[Pi,20]

FullForm[%]

Každé arbitrary precision číslo x je chápáno jako reálné číslo, které neznáme zcela přesně, ale je zatíženo jistou chybou Δx a jeho přesná hodnota leží někde v intervalu [x-Δx/2,x+Δx/2]. Přesnost čísla x je definována vztahem p=-Log[10,Δx/|x|].

Obráceně: Známe-li přesnost p, pak Δx=(10^(-p))|x|.

a=1.2`2;

10^(-2)*1.2

Přesnost tedy zhruba odpovídá počtu platných cifer daného čísla.

Při výpočtech se přesnost automaticky přepočítává, lze ji zjistit pomocí FullForm nebo Precision.

a^2

FullForm[a^2]

Precision[a^2]

FullForm[Sqrt[a]]

Precision[Sqrt[a]]

Ilustrace:

x=1;δ=.1;f[x_]:=x^2;
p1=Plot[f[x],{x,0,1.5},Epilog->{Line[{{x,0},{x,f[x]},{0,f[x]}}],Dashed,Line[{{x+δ,0},{x+δ,f[x+δ]},{0,f[x+δ]}}],Line[{{x-δ,0},{x-δ,f[x-δ]},{0,f[x-δ]}}]},AspectRatio->Automatic,PlotStyle->Thick];
x=1;δ=.1;f[x_]:=Sqrt[x];
p2=Plot[f[x],{x,0,1.5},Epilog->{Line[{{x,0},{x,f[x]},{0,f[x]}}],Dashed,Line[{{x+δ,0},{x+δ,f[x+δ]},{0,f[x+δ]}}],Line[{{x-δ,0},{x-δ,f[x-δ]},{0,f[x-δ]}}]},AspectRatio->Automatic,PlotStyle->Thick];
GraphicsRow[{p1,p2},ImageSize->Large,Spacings->20]

Příklad: Výpočet iterací funkce f(x)=4x(1-x) (logistické zobrazení) v bodě 0.1.

f[x_]:=4*x*(1-x)

machine=NestList[f,0.1,70];

arbitrary50=NestList[f,0.1`50,70];

ListLinePlot[{machine,arbitrary50},PlotStyle->{{Thickness[0.01],Orange},{Black}}]

Výsledek v machine precision je nesprávný (kumulace chyb při velkém množství aritmetických operací).

Jak víme, že je výpočet v arbitrary precision správně?

Map[Precision,arbitrary50]

Clear[f,machine,arbitrary50]

Intervalová aritmetika

Místo sledování přesnosti je možné počítat přímo s intervaly, ve kterých jsou daná reálná čísla obsažena.

Interval[{9,11}]^2

Exp[Interval[{1,2}]]

Interval[{1,2}]+Interval[{0,2}]

Interval[x] vyrobí interval s koncovými body v machine precision, který zaručeně obsahuje zadané číslo.

0.1 v machine precision:

N[SetPrecision[0.1,Infinity],30]

0.1 reprezentováno jako interval:

a=Interval[0.1]

FullForm[a]

a^2

FullForm[a^2]

Sin[Interval[N[(10^25)]]]

Interval[N[(10^25)]]//FullForm

Funkce pro numerické výpočty, volba WorkingPrecision

({{N[expr], Převod přesného čísla na desetinné}, {NIntegrate[f[x],{x,(Subscript[x, min]),(Subscript[x, max])}], Numerická integrace (∫_(Subscript[x, min])^(Subscript[x, max])) f}, {NSum[f[i],{i,(Subscript[i, min]),(Subscript[i, max])}], Součet (∑_(i=(Subscript[i, min]))^(Subscript[i, max]))f[i]numericky}, {NProduct[f[i],{i,(Subscript[i, min]),(Subscript[i, max])}], Součin (∏_(i=(Subscript[i, min]))^(Subscript[i, max]))f[i] numericky}, {FindRoot[lhs==rhs,{x,(Subscript[x, 0])}], Přibližné řešení rovnice blízko (Subscript[x, 0])}, {NSolve[lhs==rhs,x], Všechna řešení (algebraické) rovnice}, {FindMinimum[f[x],{x,(Subscript[x, 0])}], Nalezení lokálního minima blízko (Subscript[x, 0])}, {NMinimize[f[x],x], Nalezení globálního minima}, {NDSolve[eqns,y,{x,(Subscript[x, min]),(Subscript[x, max])}],  Numerické řešení diferenciální rovnice }})

Tento integrál nelze počítat symbolicky:

Integrate[Sin[x^3]/(x^3+x^2+1),{x,0,Pi}]

N[Integrate[Sin[x^3]/(x^3+x^2+1),{x,0,Pi}]]

Stejný výsledek, ale rychleji:

NIntegrate[Sin[x^3]/(x^3+x^2+1),{x,0,Pi}]

Varovný příklad:

NIntegrate[x^20*Exp[x], {x, 0, 1}]

N[Integrate[x^20*Exp[x], {x, 0, 1}]]

N[Integrate[x^20*Exp[x], {x, 0, 1}],25]

Vysvětlení:

Integrate[x^20*Exp[x], {x, 0, 1}]

Expand[%]

N[{2432902008176640000,895014631192902121 E},25]

Integrace přes neomezený interval:

NIntegrate[Exp[-x^2],{x,-Infinity,Infinity}]

Dvojný integrál přes trojúhelníkovou oblast:

NIntegrate[Sin[x y],{x,0,1},{y,0,x}]

NIntegrate (i ostatní funkce) standardně počítá v machine precision. Jinou přesnost můžeme nastavit pomocí WorkingPrecision.

NIntegrate[Sin[x^3]/(x^3+x^2+1),{x,0,Pi},WorkingPrecision->30]

Numerické řešení rovnice (e^x)=5x:

Plot[{Exp[x],5x},{x,-1,4},PlotStyle->Thick]

FindRoot[Exp[x]==5 x,{x,0}]

FindRoot[Exp[x]==5 x,{x,2.5}]

FindRoot[Exp[x]==5 x,{x,2.5},WorkingPrecision->40]

Příklad: Taylorův polynom kosinu stupně n

poly[n_]:=Normal[Series[Cos[x], {x,0,n}]];

poly[8]

Manipulate[Plot[Evaluate[{Cos[x],poly[i]}],{x,-15,15},PlotRange->{-2,2},PlotStyle->{{Thick,Red},{Thick,Blue}}],{i,2,40,2}]

Pokus o vykreslení polynomu stupně 120 (neúspěšný):

p=poly[120];

Plot[p,{x,0, 50}, PlotRange->{-2, 2}]

Kde je problém?

Hodnota polynomu pro x=40:

p/.x->40

N[%]

Ale Plot při vyhodnocování používá machine precision!

p/.x->40.0

Náprava: Pomocí volby WorkingPrecision můžeme požádat o výpočet v arbitrary precision:

Plot[p,{x,0, 50}, PlotRange->{-2, 2},WorkingPrecision->40]

Clear[p,poly]

Příklad: Šikmý vrh (NDSolve + FindMaximum)

Jak daleko doletí předmět vystřelený pod úhlem α rychlostí v? Předpokládáme, že odpor vzduchu je přímo úměrný velikosti okamžité rychlosti (koeficient úměrnosti je k).

g=9.81;m=1;

Jde o soustavu diferenciálních rovnic (Newtonův pohybový zákon), které vyřešíme numericky pomocí NDSolve. Použitím WhenEvent (podrobnosti příště) zastavíme výpočet v okamžiku, kdy y=0.

trajectory[α_,v_,k_]:={x[t],y[t]}/.First[NDSolve[
{x''[t]==-k x'[t],
y''[t]==-m g-k y'[t],
x[0]==0,y[0]==0,x'[0]==v Cos[α],y'[0]==v Sin[α],WhenEvent[y[t]==0,"StopIntegration"]},{x[t],y[t]},{t,0,100}]]

t45=trajectory[45°,20,0.05]

ParametricPlot[t45,{t,0,2.82}]

distance[α_,v_,k_]:=Module[{x=trajectory[α,v,k][[1]]},x/.t->x[[0,1,1,2]]]

distance[45°,20,0.05]

Pro jaký úhel α doletí předmět nejdále? Stačí použít FindMaximum.

f[α_?NumericQ]:=distance[α,20,0.05]

Plot[f[α],{α,0,90°},AxesOrigin->{0,0},Ticks->{Table[{i Degree,StringForm["``°",i]},{i,0,90,10}],Automatic}]

FindMaximum akceptuje i funkce definované „programem” (tj. nikoliv vzorcem v uzavřeném tvaru):

FindMaximum[f[α],{α,45°}]

(α/.%[[2]])/Degree

FindMaximum implicitně používá metodu tečen a potřebuje znát derivace zkoumané funkce.
Pro funkce definované programem nejsou derivace k dispozici (jako argumenty funkce f připouštíme jen čísla, Mathematica počítá derivaci numericky), proto je lepší použít metodu sečen – stačí zadat dvě počáteční hodnoty (viz dokumentaci).

FindMaximum[f[α],{α,40°,45°}]

(α/.%[[2]])/Degree

Jaká je závislost optimálního úhlu na koeficientu tření k?

f[α_?NumericQ,k_?NumericQ]:=distance[α,20,k]

optimalAngle[k_?NumericQ]:=(α/.FindMaximum[f[α,k],{α,40°,45°}][[2]])/Degree

optimalAngle[0.05]

Plot[optimalAngle[k],{k,0,0.2},AxesOrigin->{0,0}]

Clear[f,optimalAngle,distance,g,m,trajectory,t45]

Sledování výpočtů pomocí StepMonitor a EvaluationMonitor

Volby StepMonitor a EvaluationMonitor umožňují sledovat průběh numerického výpočtu.

Výpočet se často skládá z mnoha kroků. Pomocí volby StepMonitor lze zadat výraz nebo příkaz, který má být v každém kroku vyhodnocen.

Příklad: Numerické řešení rovnice s výpisem všech postupných aproximací.

FindRoot[Cos[x]==x,{x,1},StepMonitor:>Print[x]]

Poznámka: Při zadávání StepMonitor je potřeba místo -> použít :> (lokální přepisovací pravidlo s odloženým vyhodnocením).

Chceme-li s daty dále pracovat, je vhodné je ukládat do seznamu. K tomu dobře poslouží dvojice funkcí Sow a Reap.

?Sow

?Reap

Reap[FindRoot[Cos[x]==x,{x,1},StepMonitor:>Sow[x]]]

EvaluationMonitor specifikuje, co se má provést při každém vyhodnocování vstupních dat.

StepMonitor vs. EvaluationMonitor: Jeden krok může zahrnovat více vyhodnocení.

Reap[FindRoot[Cos[x]==x,{x,5},StepMonitor:>Sow[x]]]

Reap[FindRoot[Cos[x]==x,{x,5},EvaluationMonitor:>Sow[x]]]

Příklad: Numerická integrace a nalezení všech bodů, ve kterých je integrovaná funkce vyhodnocována.

Integrate[1/Sqrt[x],{x,0,1}]

NIntegrate[1/Sqrt[x],{x,0,1}]

xHodnoty=Reap[NIntegrate[1/Sqrt[x],{x,0,1},EvaluationMonitor:>Sow[x]]][[2,1]]

ListPlot[xHodnoty]

Příklad: Nalezení všech kontrolních bodů, které použije Plot při vykreslování grafu funkce.

Mesh->All zvýrazní kontrolní body, ale jejich souřadnice neznáme:

Plot[Sin[x],{x,0,4Pi},Mesh->All,MeshStyle->PointSize[0.015],PlotStyle->{Thick,Red}]

Najdeme je pomocí EvaluationMonitor, Reap a Sow:

Reap[Plot[Sin[x],{x,0,4Pi},EvaluationMonitor:>Sow[x]]]

xHodnoty=Sort[%[[2,1]]];

Length[xHodnoty]

yHodnoty=Sin[xHodnoty];

Zkouška: Lomená čára procházející nalezenými body.

ListLinePlot[Transpose[{xHodnoty,yHodnoty}]]

Clear[xHodnoty,yHodnoty]

Aplikace: Nalezení všech nulových bodů funkce v zadaném intervalu

FindAllRoots[f_,{x_,a_,b_}]:=Module[{xValues,seedPositions,roots},xValues=Sort[Reap[Plot[f,{x,a-(b-a)/100,b+(b-a)/100},EvaluationMonitor:>Sow[x]]][[2,1]]];
seedPositions=Select[Range[Length[xValues]-1],(Sign[f/.x->xValues[[#]]]Sign[f/.x->xValues[[#+1]]]<0)&];
roots=Map[x/.FindRoot[f,{x,xValues[[#]]}] & ,seedPositions];
Select[Union[roots],a≤#≤b &]]

solns=FindAllRoots[(x/5)  Cos[x]+Sin[(x^3)],{x,-3,3}]

Plot[(x/5)  Cos[x]+Sin[(x^3)],{x, -3, 3}, Epilog->{PointSize[Large],Map[Point[{#,0}] &,solns]}]

solns=FindAllRoots[Sin[(x^2)]-(E^((Cos[x])/2))+1,{x,-10,0}]

Plot[Sin[(x^2)]-(E^((Cos[x])/2))+1,{x,-10,0},Epilog->{PointSize[Large],Map[Point[{#,0}] &,solns]}]

FindAllRoots[x^2,{x,-10,10}]

Clear[solns,FindAllRoots]

Ukládání výsledků

Mathematica si pamatuje výsledky některých výpočtů:

N[Pi,1000000];//Timing

N[Pi,1000000];//Timing

Komentář z dokumentace k Timing:

Timing[expr] may give different results on different occasions within a session. A typical reason is the use of internal system caches, which can be cleared with ClearSystemCache.