Numerické výpočty

Antonín Slavík
Matematicko-fyzikální fakulta 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 přesnost čísla (bude podrobněji vysvětleno dále).

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]

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]

Graphics[{Line[{{0,0},{100,0}}],
MapThread[{Text[#1,{#2,3},{-1,0},{0,1}],Line[{{#2,-1},{#2,1}}]}&,Transpose@{{0,0},{"$MinMachineNumber",2},{"$MaxMachineNumber",100},{"1/2",10},{"1/2+$MachineEpsilon/2",11},{"1/2+$MachineEpsilon",12},{1,20},{"1+$MachineEpsilon",22},{"1+2$MachineEpsilon",24},{2,40},{"2+2$MachineEpsilon",44},{"2+4$MachineEpsilon",48},{4,80}}]},ImageSize->Full,BaseStyle->12,PlotLabel->Style["Nezáporná čísla reprezentovatelná ve strojové přesnosti",FontSize->24,FontWeight->Bold]]

{$MinMachineNumber,$MaxMachineNumber}

Log[2,%]

$MinMachineNumber = (Subscript[(1.0⋯0), 2])·(2^(-1022)), $MaxMachineNumber = (Subscript[(1.1⋯1), 2])·(2^1023) ≐ (2^1024)

?$MachineEpsilon

$MachineEpsilon

Log[2,$MachineEpsilon]

$MachineEpsilon = (Subscript[(1.0⋯01), 2])-(Subscript[(1.0⋯0), 2])=(2^(-52))

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.

Další příklady

Závislost na pořadí operací:

10.^20+1-10.^20

10.^20-10.^20+1

Vysvětlení:

IntegerDigits[10^20,2]

Length[%]

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í nezáporná čísla, se kterými lze počítat v arbitrary precision:

{$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:

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

Příklad: Výpočet iterací funkce f(x)=4x(1-x) 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]

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,π}]

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

Stejný výsledek, ale rychleji:

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

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,π},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]

Varovné příklady

WorkingPrecision->n znamená, že výpočet probíhá s čísly v přesnosti n. V průběhu výpočtu však obvykle dochází k zaokrouhlovacím chybám a přesnost finálního výsledku může být mnohem nižší než n.

Integrate[Sin[x]/x,{x,0,∞}]

N[(π/2),20]

NIntegrate[Sin[x]/x,{x,0,∞},WorkingPrecision->20]

NIntegrate[Sin[x]/x,{x,0,∞},WorkingPrecision->30]

NIntegrate[Sin[x]/x,{x,0,∞},WorkingPrecision->35]

Pokud počítáme s vyšší než strojovou přesností, nesmí se na vstupu objevit čísla ve strojové přesnosti.

NSum[(1/(n^(3/2))),{n,1,∞},WorkingPrecision->30]

NSum[(1/(n^1.5)),{n,1,∞},WorkingPrecision->30]

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,4 π},Mesh->All,MeshStyle->PointSize[0.015],PlotStyle->{Thick,Red}]

Najdeme je pomocí EvaluationMonitor, Reap a Sow:

Reap[Plot[Sin[x],{x,0,4 π},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

Myšlenka: Necháme funkci Plot sestrojit graf zadané funkce. Procházíme výslednou lomenou čáru a hledáme úsečky, které protínají osu x (opačná znaménka funkčních hodnot v krajních bodech). Tím získáme přibližné polohy kořenů, které použijeme jako počáteční odhady ve FindRoot.

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]}]

Varování: Algoritmus nenajde kořeny, ve kterých se graf funkce pouze dotýká osy x.

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

Clear[solns,FindAllRoots]

Cvičení

Nad vodorovnou úsečkou o délce l=5280 cm je sestrojen oblouk kružnice o délce t=5281 cm. Zjistěte, v jaké výšce h nad úsečkou se nachází nejvyšší bod kružnice. Najděte výsledek s přesností minimálně 10 platných cifer.
Návod. Vyjádřete vztahy mezi veličinami t, l, h, d, r, θ. Hodnotu θ získejte numerickým řešením vhodné rovnice, poté dopočítejte r a následně h.

Plot[(Sqrt[625-(x^2)]),{x,-15,15},PlotRange->{{-20,20},{-1,30}},Axes->False,Epilog->{Line[{{0,0},{-15,20},{15,20},{0,0},{0,25}}],Circle[{0,0},4.5,{.94,(π/2)}],
Text["h",{0.6,22.5},{-1,0}],Text["t",{0,27}],Text["θ",{1,3}],Text["r",{9,10}],Text["l/2",{7.5,18.6},{0,1}],Text["d",{.9,11}]},BaseStyle->12, ImageSize->Medium]

Uvažujme funkci f(a)=(2+sin(10 a)) (∫_0^2)(x^a) sin((a/(2-x)))dx. Sestrojte její graf na intervalu [0,5]. Najděte její globální maximum na tomto intervalu s přesností minimálně 12 platných cifer. Mějte na paměti, že WorkingPrecision->n znamená, že výpočet probíhá s čísly v přesnosti n, avšak přesnost finálního výsledku může být mnohem nižší než n (viz dokumentaci).