[NumAn] PC-Oef3

2274 days ago by Niels.Neirynck

Oplossen van stelsel van lineaire vergelijkingen

Sage heeft een reeks van ingebouwde functies voor berekeningen uit de lineaire algebra. Een inleiding kan men vinden op http://www.sagemath.org/doc/tutorial/tour_linalg.html. Op de site http://wiki.sagemath.org/quickref kan men een overzichtsdocument vinden van de bewerkingen voor lineaire algebra: quickref-linalg.pdf

Je zal onmiddellijk heel wat vertrouwde namen herkennen als je het bovenstaande document bekijkt. Wil je meer weten over een bepaalde procedure, dan kan je de help-functie als volgt oproepen:

m = matrix([[1,2],[3,4]]); #bijvoorbeeld m.LU() 
       
(
[0 1]  [  1   0]  [  3   4]
[1 0], [1/3   1], [  0 2/3]
)
(
[0 1]  [  1   0]  [  3   4]
[1 0], [1/3   1], [  0 2/3]
)
m.LU? #voer uit om de helppagina te bekijken 
       

We zullen hier opnieuw het stelsel $Hx = b$ beschouwen, met $H$ de n-de orde Hilbert-matrix, en $b$ het rechterlid uit Oefening 3.2.

Er is (momenteel) geen ingebouwde constructiefunctie aanwezig in Sage voor de Hilbert-matrix. Gelukkig kan dit vrij eenvoudig in python:

def hilbert_matrix(n): return matrix([[1/i for i in [j..j+(n-1)]] for j in [1..n]]) 
       
We nemen als eerste geval n=10 (maar uiteraard kan je verderop met ander n-waarden experimenteren).
n = 10; H = hilbert_matrix(n); H 
       
[   1  1/2  1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10]
[ 1/2  1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11]
[ 1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12]
[ 1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13]
[ 1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14]
[ 1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15]
[ 1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16]
[ 1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17]
[ 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18]
[1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19]
[   1  1/2  1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10]
[ 1/2  1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11]
[ 1/3  1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12]
[ 1/4  1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13]
[ 1/5  1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14]
[ 1/6  1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15]
[ 1/7  1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16]
[ 1/8  1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17]
[ 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18]
[1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19]

We berekenen de vector $b$, waarbij $b_i$ de som is van de $i$de rij in $H$

b = vector([sum(H[i]) for i in range(n)]).column(); show(b) 
       

                                
                            

                                
x= vector([ 1 for _ in range(n)]).column(); show(x) 
       

                                
                            

                                
show(H*x - b) #voor controle 
       

                                
                            

                                

Dus $Hx = b$. We herberekenen de waarde van $x$ door  $x$ in het stelsel $Hx = b$ op te lossen met behulp van  ingebouwde functie solve_right:

H.solve_right(b) 
       
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]
[1]

Vervolgens doen we alle vorige berekeningen opnieuw in vlottende-puntvoorstelling.

Digits = 10; H1 = N(H , digits=Digits); H1 
       
[  1.000000000  0.5000000000  0.3333333333  0.2500000000  0.2000000000 
0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000]
[ 0.5000000000  0.3333333333  0.2500000000  0.2000000000  0.1666666667 
0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091]
[ 0.3333333333  0.2500000000  0.2000000000  0.1666666667  0.1428571429 
0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333]
[ 0.2500000000  0.2000000000  0.1666666667  0.1428571429  0.1250000000 
0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692]
[ 0.2000000000  0.1666666667  0.1428571429  0.1250000000  0.1111111111 
0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143]
[ 0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000
0.09090909091 0.08333333333 0.07692307692 0.07142857143 0.06666666667]
[ 0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091
0.08333333333 0.07692307692 0.07142857143 0.06666666667 0.06250000000]
[ 0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333
0.07692307692 0.07142857143 0.06666666667 0.06250000000 0.05882352941]
[ 0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692
0.07142857143 0.06666666667 0.06250000000 0.05882352941 0.05555555556]
[ 0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143
0.06666666667 0.06250000000 0.05882352941 0.05555555556 0.05263157895]
[  1.000000000  0.5000000000  0.3333333333  0.2500000000  0.2000000000  0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000]
[ 0.5000000000  0.3333333333  0.2500000000  0.2000000000  0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091]
[ 0.3333333333  0.2500000000  0.2000000000  0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333]
[ 0.2500000000  0.2000000000  0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692]
[ 0.2000000000  0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143]
[ 0.1666666667  0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143 0.06666666667]
[ 0.1428571429  0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143 0.06666666667 0.06250000000]
[ 0.1250000000  0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143 0.06666666667 0.06250000000 0.05882352941]
[ 0.1111111111  0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143 0.06666666667 0.06250000000 0.05882352941 0.05555555556]
[ 0.1000000000 0.09090909091 0.08333333333 0.07692307692 0.07142857143 0.06666666667 0.06250000000 0.05882352941 0.05555555556 0.05263157895]
parent(H1) 
       
Full MatrixSpace of 10 by 10 dense matrices over Real Field with 37 bits
of precision
Full MatrixSpace of 10 by 10 dense matrices over Real Field with 37 bits of precision
b1 = vector([sum(H1[i]) for i in range(n)]).column(); b1 
       
[ 2.928968254]
[ 2.019877345]
[ 1.603210678]
[ 1.346800422]
[ 1.168228993]
[ 1.034895660]
[0.9307289932]
[0.8466953798]
[0.7772509353]
[0.7187714032]
[ 2.928968254]
[ 2.019877345]
[ 1.603210678]
[ 1.346800422]
[ 1.168228993]
[ 1.034895660]
[0.9307289932]
[0.8466953798]
[0.7772509353]
[0.7187714032]
x1_calc = H1.solve_right(b1); show(x1_calc) 
       

                                
                            

                                

De onstabiliteit van de Hilbert-matrix wordt hierbij duidelijk. Je kan eventueel andere waarden aan n en aan Digits toekennen, en verder experimenten.

Terloops kunnen we ook even de conditiegetallen van enkele Hilbert-matrices berekenen. Er is geen standaard functie aanwezig in Sage om symbolisch het conditiegetal van een matrix te berekenen. Gebaseerd op   $\text{cond}(A) = ||A|| \cdot ||A^{-1}||$ , gebruiken we de volgende functie om het conditiegetal (op inefficiënte wijze)  symbolisch te berekenen:

def symcond(HH): return HH.norm(Infinity) * HH.inverse().norm(Infinity); symcondlist = [ symcond(hilbert_matrix(i)) for i in range(5,16)] symcondlist 
       
[943656.0,
 29070278.999999996,
 985194886.4999999,
 33872791094.999996,
 1099654541342.5,
 35357439251992.0,
 1233702357598850.2,
 4.115445402289639e+16,
 1.324409009034709e+18,
 4.53775784394382e+19,
 1.539191562955312e+21]
[943656.0,
 29070278.999999996,
 985194886.4999999,
 33872791094.999996,
 1099654541342.5,
 35357439251992.0,
 1233702357598850.2,
 4.115445402289639e+16,
 1.324409009034709e+18,
 4.53775784394382e+19,
 1.539191562955312e+21]

Sage kan beroep doen op het pakket NumPy om zuivere numerieke berekeningen uit te voeren . Deze berekeningen gebeuren met machinegetallen (64 bit fp-getallen). We gebruiken de functie cond() uit het pakket NumPy om het conditiegetal van Hilbert-matrices numeriek te bepalen.

import numpy from numpy.linalg.linalg import cond numcondlist = [cond(hilbert_matrix(i),p=numpy.inf) for i in range(5,16)]; numcondlist 
       
[943655.99999993353,
 29070279.002940644,
 985194889.71984804,
 33872790819.494709,
 1099650991701.0522,
 35353724553756.43,
 1230369938308719.0,
 37983201226912136.0,
 4.2759533532683162e+17,
 5.9615764639409674e+18,
 8.0290366187793702e+17]
[943655.99999993353,
 29070279.002940644,
 985194889.71984804,
 33872790819.494709,
 1099650991701.0522,
 35353724553756.43,
 1230369938308719.0,
 37983201226912136.0,
 4.2759533532683162e+17,
 5.9615764639409674e+18,
 8.0290366187793702e+17]
show( matrix(vector(symcondlist)).stack(vector(numcondlist)).transpose()) 
       

                                
                            

                                

Het conditiegetal van de Hilbert-matrix stijgt dramatisch wanneer $n$ toeneemt. De numerieke berekeningen van de conditiegetallen zijn zelf onderhevig aan fouten naarmate $n$ toeneemt.

Bij hogere waarden van $n$ is de Hilbert-matrix zodanig onstabiel dat het numeriek berekende conditiegetal enorme fouten vertoont. Vaak is het zo dat men het conditiegetal niet exact nodig heeft, een schatting is voldoende.

Je kan je natuurlijk de vraag stellen: waarom nog berekeningen uitvoeren in een programmeertaal zoals Java, Pascal, Fortran, C, enz. (waar men met vlottende-puntvoorstelling werkt) als men de mogelijkheid heeft (zoals in Sage, Maple, Mathematica, enz.) om exact te werken? Hiertoe zijn verschillende redenen. Vooreerst is het zo dat de gegevens voor een praktisch probleem (zoals een lineair stelsel oplossen) meestal uit een experiment of meting komen, en dus slechts als vlottende-puntgetallen kunnen beschouwd worden. Verder is het vooral een kwestie van ruimte en tijd. In werkelijke toepassingen kunnen stelsels enorm groot zijn (n=100, 1000, 10000, ...), en een exacte voorstelling van de gegevens zou te veel geheugenruimte opeisen. Bovendien vereist de praktijk meestal dat de berekening snel gebeurt (b.v., als het resultaat dient om een machine bij te sturen kan dit niet eeuwig duren, of de bijsturing gebeurt te laat), en dit is niet het geval voor exacte berekeningen. Beschouw bij voorbeeld n=100 (wat dan nog relatief klein is), en voer opnieuw bovenstaande berekening uit.

Men kan Sage berekeningen die te lang duren onderbreken door bovenaan:  "Action.." >  "Interrupt"     of    "Action.." >  "restart worksheet"  uit te voeren.

#opbouw van matrix H en vector b: n = 100; H = hilbert_matrix(n); b = vector([sum(H[i]) for i in range(n)]).column(); 
       

Symbolisch oplossen naar $x$ via solve_right

time x_sym=H.solve_right(b); 
       
Time: CPU 0.22 s, Wall: 0.22 s
Time: CPU 0.22 s, Wall: 0.22 s

Numeriek oplossen naar $x$ via solve in het pakket NumPy. De getallen worden door de functie eerst omgezet naar machinegetallen.

from numpy import linalg time x_num=linalg.solve(H,b); 
       
Time: CPU 0.08 s, Wall: 0.08 s
Time: CPU 0.08 s, Wall: 0.08 s