Stabiele en onstabiele berekeningen: demonstratie in Sage

In Sage gebeuren bewerkingen met gehele of rationale getallen automatisch symbolisch, dus "exact", tenzij anders gespecifieerd. Stel dat we dezelfde functies $f(x)$ en $g(x)$ uit voorbeeld 1 kiezen, m.a.w.

var('x') f(x) = 1/(1+2*x) - (1+x)/(1+3*x); g(x) = -2*x^2/((1+2*x)*(1+3*x));

Als we $f(x)$ en $g(x)$ berekenen voor een zekere $x$ krijgen we uiteraard dezelfde waarden:

show([f(10^(-2)) , g(10^(-2))]) f(10^(-2)) - g(10^(-2))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{5253}, -\frac{1}{5253}\right]$

We vragen het type op van de functiewaarde en we verifiëren dat de bewerkingen inderdaad symbolisch worden uitgevoerd:

val = f(10^(-4)) parent(val)

Symbolic Ring

Symbolic Ring

Als we $f(x)$ en $g(x)$ berekenen voor $x = 10^{-k}$ met $k = 1, 2, \dots, 20 \quad$ worden de uitdrukkingen echter heel omslachtig.

show([f(10^(-i)) for i in range(20)])

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{6}, -\frac{1}{78}, -\frac{1}{5253}, -\frac{1}{502503}, -\frac{1}{50025003}, -\frac{1}{5000250003}, -\frac{1}{500002500003}, -\frac{1}{50000025000003}, -\frac{1}{5000000250000003}, -\frac{1}{500000002500000003}, -\frac{1}{50000000025000000003}, -\frac{1}{5000000000250000000003}, -\frac{1}{500000000002500000000003}, -\frac{1}{50000000000025000000000003}, -\frac{1}{5000000000000250000000000003}, -\frac{1}{500000000000002500000000000003}, -\frac{1}{50000000000000025000000000000003}, -\frac{1}{5000000000000000250000000000000003}, -\frac{1}{500000000000000002500000000000000003}, -\frac{1}{50000000000000000025000000000000000003}\right]$

Zoals je ziet, in dit eenvoudig voorbeeld valt het nog mee; maar in het algemeen worden exacte rationale uitdrukkingen al snel heel ingewikkeld. Om een beter idee te hebben over het gedrag van een rij getallen kan men overgaan naar de vlottende-puntvoorstelling. In Sage kan men dit doen met de functie N() (andere benamingen: n(), numerical_approx()). De vorige rij wordt dan:

list = [N(f(10^(-i))) for i in range(20)]; list

[-0.166666666666667,
 -0.0128205128205128,
 -0.000190367409099562,
 -1.99003787042067e-6,
 -1.99900037987004e-8,
 -1.99990000379987e-10,
 -1.99999000003800e-12,
 -1.99999900000038e-14,
 -1.99999990000000e-16,
 -1.99999999000000e-18,
 -1.99999999900000e-20,
 -1.99999999990000e-22,
 -1.99999999999000e-24,
 -1.99999999999900e-26,
 -1.99999999999990e-28,
 -1.99999999999999e-30,
 -2.00000000000000e-32,
 -2.00000000000000e-34,
 -2.00000000000000e-36,
 -2.00000000000000e-38]

[-0.166666666666667,
 -0.0128205128205128,
 -0.000190367409099562,
 -1.99003787042067e-6,
 -1.99900037987004e-8,
 -1.99990000379987e-10,
 -1.99999000003800e-12,
 -1.99999900000038e-14,
 -1.99999990000000e-16,
 -1.99999999000000e-18,
 -1.99999999900000e-20,
 -1.99999999990000e-22,
 -1.99999999999000e-24,
 -1.99999999999900e-26,
 -1.99999999999990e-28,
 -1.99999999999999e-30,
 -2.00000000000000e-32,
 -2.00000000000000e-34,
 -2.00000000000000e-36,
 -2.00000000000000e-38]

Laten we een element uit die lijst nemen (het eerste) en het type inspecteren:

parent(list[0]);

Real Field with 53 bits of precision

Real Field with 53 bits of precision

Het getal behoort tot de structuur van vlottende-puntgetallen met 53 bits in de mantisse. Dit is de standaard voor een 64-bit computer. In java gebruiken we daarvoor het type double. Oudere 32-bit computers gebruiken standaard vlottende-puntgetallen met 23 bits in de mantisse (float in java). Bij de omzetting van symbolische uitdrukkingen naar vlottende-puntgetallen kunnen we het aantal precisiebits vrij kiezen in sage (prec staat voor precision):

num = N(f(10^(-2)),prec=112) parent(num)

Real Field with 112 bits of precision

Real Field with 112 bits of precision

We kunnen ook de precisie van een omzetting vastleggen door op te geven hoeveel cijfers na de komma (na normalizatie) we wensen in een decimale voorstelling. Stel dat we omzetting willen met 10 cijfers na de komma:

num = N(f(10^(-2)),digits=10) show(num) parent(num)

$\newcommand{\Bold}[1]{\mathbf{#1}}-0.0001903674091$

Real Field with 37 bits of precision

Real Field with 37 bits of precision

Om een precisie te bekomen van 10 decimale cijfers plaatst Sage dit getal in een vlottende-puntsvoorstelling met 37 bits.

Laten we nu eens simuleren hoe de berekening van $f(x)$ in vlottende-puntvoorstelling zou gebeuren in een klassieke programmeeromgeving. We beginnen met $f(x)$ symbolisch te berekenen en vervolgens zetten we de waarde om naar een vlottende-puntvoorstelling. We vergelijken die waarde met de waarde waarbij we eerst $x$ omzetten naar een vlottende-puntvoorstelling en vervolgens $f$ op toepassen. We bekijken het resultaat hiervan voor verschillende mantisselengtes. Eerst stellen we een hulpfunctie $ff(x)$ op die we gebruiken voor het simuleren van een klassieke programmeeromgeving:

De functie constructie die we hebben gebruikt bij het aanmaken van de functie $f$ is eerder geschikt voor symbolische bewerkingen. Sage zal dan ook proberen opgegeven waarden symbolisch te behandelen. Voor het simuleren van een klassieke programmeeromgeving gebruiken we de python def-constructie voor functies:

def ff(z): return 1/(1+2*z) - (1+z)/(1+3*z) val = ff(N(10^(-3))) parent(val)

Real Field with 53 bits of precision

Real Field with 53 bits of precision

Digits = 10; #aantal cijfers na de komma. nrElements = 20; symlist = [N( f(10^(-i)) ,digits=Digits) for i in range(nrElements)]; numlist = [ff( N(10^(-i),digits=Digits) ) for i in range(nrElements)]; for i in range(nrElements): print('{:>30} {:>30}'.format(symlist[i], numlist[i]))

                 -0.1666666667                  -0.1666666667
                -0.01282051282                 -0.01282051282
              -0.0001903674091               -0.0001903674056
               -1.990037870e-6                -1.990039891e-6
               -1.999000380e-8                -1.999433152e-8
              -1.999900004e-10               -2.110027708e-10
              -1.999990000e-12                   0.0000000000
              -1.999999000e-14                   0.0000000000
              -1.999999900e-16                1.455191523e-11
              -1.999999990e-18                   0.0000000000
              -1.999999999e-20                   0.0000000000
              -2.000000000e-22                   0.0000000000
              -2.000000000e-24                   0.0000000000
              -2.000000000e-26                   0.0000000000
              -2.000000000e-28                   0.0000000000
              -2.000000000e-30                   0.0000000000
              -2.000000000e-32                   0.0000000000
              -2.000000000e-34                   0.0000000000
              -2.000000000e-36                   0.0000000000
              -2.000000000e-38                   0.0000000000

                 -0.1666666667                  -0.1666666667
                -0.01282051282                 -0.01282051282
              -0.0001903674091               -0.0001903674056
               -1.990037870e-6                -1.990039891e-6
               -1.999000380e-8                -1.999433152e-8
              -1.999900004e-10               -2.110027708e-10
              -1.999990000e-12                   0.0000000000
              -1.999999000e-14                   0.0000000000
              -1.999999900e-16                1.455191523e-11
              -1.999999990e-18                   0.0000000000
              -1.999999999e-20                   0.0000000000
              -2.000000000e-22                   0.0000000000
              -2.000000000e-24                   0.0000000000
              -2.000000000e-26                   0.0000000000
              -2.000000000e-28                   0.0000000000
              -2.000000000e-30                   0.0000000000
              -2.000000000e-32                   0.0000000000
              -2.000000000e-34                   0.0000000000
              -2.000000000e-36                   0.0000000000
              -2.000000000e-38                   0.0000000000

Experimenteer nu met andere waarden van Digits en vergelijk.

#De weergavecode in de vorige blok faalt bij hogere waarden van 'Digits' #Voor hoge waarden kan je de waarden weergeven als een matrix: #We steken de twee lijsten als vector in een matrix en gebruik de show() functie: m = matrix(vector(symlist)); m = m.stack(vector(numlist)); show(m.transpose()); #gebruik dit indien de vorige weergavecode faalt bij hoge waarden van 'Digits'

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} -0.1666666667 & -0.1666666667 \\ -0.01282051282 & -0.01282051282 \\ -0.0001903674091 & -0.0001903674056 \\ -1.990037870 \times 10^{-6} & -1.990039891 \times 10^{-6} \\ -1.999000380 \times 10^{-8} & -1.999433152 \times 10^{-8} \\ -1.999900004 \times 10^{-10} & -2.110027708 \times 10^{-10} \\ -1.999990000 \times 10^{-12} & 0.0000000000 \\ -1.999999000 \times 10^{-14} & 0.0000000000 \\ -1.999999900 \times 10^{-16} & 1.455191523 \times 10^{-11} \\ -1.999999990 \times 10^{-18} & 0.0000000000 \\ -1.999999999 \times 10^{-20} & 0.0000000000 \\ -2.000000000 \times 10^{-22} & 0.0000000000 \\ -2.000000000 \times 10^{-24} & 0.0000000000 \\ -2.000000000 \times 10^{-26} & 0.0000000000 \\ -2.000000000 \times 10^{-28} & 0.0000000000 \\ -2.000000000 \times 10^{-30} & 0.0000000000 \\ -2.000000000 \times 10^{-32} & 0.0000000000 \\ -2.000000000 \times 10^{-34} & 0.0000000000 \\ -2.000000000 \times 10^{-36} & 0.0000000000 \\ -2.000000000 \times 10^{-38} & 0.0000000000 \end{array}\right)$

[NumAn] PC-Oef1

2274 days ago by Niels.Neirynck

Stabiele en onstabiele berekeningen: demonstratie in Sage