Initiation au Calcul Scientifique Intensif

  • Rappels sur les matrices

  • Structure des matrices

Illustration du remplissage par la multiplication matricielle

In [1]:
import matplotlib.pylab as plt
import scipy.sparse as sps
import scipy.linalg as spl
import numpy as np

n=100
A = sps.rand(n,n, density=0.05)
show('creux : ', (len(sps.find(A)[0])/n^2*100).n(),'%')
A=A.todense()

plt.spy(A,precision=0.01, markersize=1)
plt.show();
In [2]:
B=np.dot(A,A)
show('creux : ', (len(sps.find(B)[0])/n^2*100).n(),'%')

plt.spy(B,precision=0.01, markersize=1)
plt.show();
In [3]:
C=np.dot(B,A)
show('creux : ', (len(sps.find(C)[0])/n^2*100).n(),'%')

plt.spy(C,precision=0.01, markersize=1)
plt.show();
In [4]:
D=np.dot(C,A)
show('creux : ', (len(sps.find(D)[0])/n^2*100).n(),'%')

plt.spy(D,precision=0.01, markersize=1)
plt.show();

Illustration de l'effet de la renumérotation

In [5]:
# Matrice aléatoire au format CSR
ran=sps.random(100,100,.05, format='csr')
A = sps.csr_matrix((np.ones(ran.data.shape).astype(np.uint8),ran.indices, ran.indptr))
plt.spy(A,precision=0.01, markersize=1)
plt.show();

# Largeur de bande
i,j=A.nonzero()
show(np.max(i-j))

# Tableau de renumérotation après algorithme Reverse Cuthill MacKee
perm = sps.csgraph.reverse_cuthill_mckee(A, False)
# Permutation
B=A[np.ix_(perm,perm)].A

plt.spy(B,precision=0.01, markersize=1)
plt.show();

# Largeur de bande
i,j=B.nonzero()
show(np.max(i-j))

Illustration de l'effet du conditionnement

In [6]:
A= matrix(4,4,[10,7,8,7,7,5,6,5,8,6,10,9,7,5,9,10])
b1=vector([32,23,33,31])
show('A=',A,'b=',b1)
In [7]:
sol1=A\b1
show(sol1)
In [8]:
b2=vector([32.1,22.9,33.1,30.9])
show('b=',b2)
In [9]:
sol2=A\b2
show(sol2.n())
In [10]:
err1=(b1-b2).norm()/b1.norm()
show(err1.n())
err2=(sol1-sol2).norm()/sol1.norm()
show(err2.n())

show(err2.n()/err1.n())

L'erreur commise sur la solution du système perturbée est 2460 fois plus importante que l'erreur commise sur le second membre !

Calculons le conditionnement de la matrice : A est symétrique et ses valeurs propres sont positives, son conditionnement en norme 2 est le rapport des la valeur propre maximale sur la valeur propre minimale.

In [11]:
show(A.eigenvalues())
condA=max(A.eigenvalues())/min(A.eigenvalues())
show(condA)

Le résultat était donc à prévoir car la matrice est mal conditionnée...

1 - Inversion de systèmes linéaires

1-1- Méthodes directes par factorisation

1-1-2- Décomposition LU

Illustration du remplissage des matrices par la méthode LU

In [12]:
import numpy as np
A=np.eye(n,n);
for i in range(n):
    A[i,n-1]=1
    A[n-1,i]=1
    
plt.spy(A)
plt.show()

P, L, U = spl.lu(A)

plt.spy(L)
plt.show()

plt.spy(U)
plt.show()
In [13]:
import numpy as np
A=np.eye(n,n);
for i in range(n):
    A[i,0]=1
    A[0,i]=1
    
plt.spy(A)
plt.show()

P, L, U = spl.lu(A)

plt.spy(L)
plt.show()

plt.spy(U)
plt.show()

Illustration de l'importance du pivotage

In [14]:
M=matrix([[0,1],[1,1]])
y=matrix([[1],[2]])
show(M,'x=',y)

x=M\y
show('x=',x)
In [15]:
var('eps')

A=matrix([[eps,1],[1,1]])
a=matrix([[1],[2]])
show(A,'x1=',a)

B=matrix([[1,1],[eps,1]])
b=matrix([[2],[1]])
show(B,'x2=',b)
In [16]:
x1=(A\a); show('x1=',x1)
x2=(B\b); show('x2=',x2)

En terme de représentation machine, le premier résultat est (0,1), tandis que le deuxième est bien (1,1) !

1-1-3- Décomposition QR

Illustration du procédé de Gram-Schmidt modifié

Procédure de Gram-Schmidt pour une matrice 3x3.

Voir exercice 2 pour les algorithmes.

In [17]:
def GramSchmidt(A):
    v1=vector(A.transpose()[0])
    v2=vector(A.transpose()[1])
    v3=vector(A.transpose()[2])
    
    u1=v1; 
    n1=u1.norm(); 
    q1=(u1/n1)

    u2=v2 - u1.dot_product(v2)/n1^2*u1;
    n2=u2.norm()
    q2=(u2/n2)

    u3=v3 - u1.dot_product(v3)/n1^2*u1 - u2.dot_product(v3)/n2^2*u2;
    n3=u3.norm()
    q3=(u3/n3)
    
    r1=u1.dot_product(v3)/n1^2
    r2=u2.dot_product(v3)/n2^2
    
    # Base orthogonale obtenue : (q1,q2,q2)
    return q1, q2, q3, r1, r2
In [18]:
var('epsilon')
assume(epsilon>0)
A=matrix([[1,1,1],[epsilon,0,0],[0,epsilon,0]]);
show('A=',A)
In [19]:
[q1,q2,q3,r1,r2]=GramSchmidt(A)
show(q1.simplify_full(),q2.simplify_full(),q3.simplify_full())
show('Verification de l orthogonalite')
show('(q1,q2)=',(q1.dot_product(q2)).simplify_full())
show('(q1,q3)=',(q1.dot_product(q3)).simplify_full())
show('(q2,q3)=',(q2.dot_product(q3)).simplify_full())

Regardons de plus près les valeurs suivantes

r1 = proj_{u1}(v3) = < u1,v3 >/< u1,u1 >

r2 = proj_{u2}(v3) = < u2,v3 >/< u2,u2 >

In [20]:
show('proj_{u1}(v3)')
show(r1)
show(limit(r1,epsilon=0))
show(r1(epsilon=1e-10))
In [21]:
show('proj_{u2}(v3)=')
show(r2)
show(r2.simplify_full())
show(limit(r2,epsilon=0))
show(r2(epsilon=1e-10))

Le terme de projection de v3 sur u2 est censé valoir 1/2 mais est numériquement nul !

Le vecteur proj_{u2}(v3) ne contribue donc pas au calcul du vecteur de base comme il le devrait.

In [22]:
epsilon=1e-10
A=matrix([[1,1,1],[epsilon,0,0],[0,epsilon,0]]);

[q1,q2,q3,R1,R2]=GramSchmidt(A)
show(q1,q2,q3)
show('Verification de l orthogonalite')
show('(q1,q2)=',q1.dot_product(q2))
show('(q1,q3)=',q1.dot_product(q3))
show('(q2,q3)=',q2.dot_product(q3))

La base otenue n'est pas orthogonale !

In [23]:
def GramSchmidt_Modif(A):
    v1=vector(A.transpose()[0])
    v2=vector(A.transpose()[1])
    v3=vector(A.transpose()[2])
    
    u1=v1; 
    n1=u1.norm(); 
    q1=(u1/n1)

    u2=v2 - u1.dot_product(v2)/n1^2*u1;
    n2=u2.norm()
    q2=(u2/n2)

    u31=v3 - u1.dot_product(v3)/n1^2*u1; 
    u3= u31 - u2.dot_product(u31)/n2^2*u2;
    n3=u3.norm()
    q3=(u3/n3)
    
    R1=u1.dot_product(v3)/n1^2
    R2=u2.dot_product(u31)/n2^2
    
    return q1, q2, q3, R1, R2
In [24]:
[q1,q2,q3,R1,R2]=GramSchmidt_Modif(A)

show(q1,q2,q3)
show('Verification de l orthogonalite')
show('(q1,q2)=',q1.dot_product(q2))
show('(q1,q3)=',q1.dot_product(q3))
show('(q2,q3)=',q2.dot_product(q3))

Avec l'algorithme modifié, la base obtenue est bien orthogonale.

In [25]:
show(r2.simplify_full())
show(R2)

En effet, dans ce cas, le terme de projection proj_{u2}(v3) est correctement approché.

1-2- Méthodes itératives classiques

Illustration de l'influence de omega sur la convergence de l'algorithme de la relaxation

In [26]:
A=matrix(5,[1., 0.1, 0.01, 0.001, 0.0001, 1., 1., 1., 1., 1., 1., 1.5, 2.25, 3.375, 5.0625, 1., 2., 4., 8., 16., 1., 3., 9., 27., 81.])
b=vector([1., 1.5, 2.25, 3.375, 5.0625])
show('A=',A)
show('b=',b)
In [27]:
A\b
Out[27]:
(0.943657635467981, 0.640855911330046, -0.871733032293372, 0.994868637110014, -0.207649151614668)

Algorithme de la relaxation

In [28]:
def relax(A,b,x0,niter,omega):
    n=(A.dimensions())[0];
    
    D=matrix(RR,n)
    E=matrix(RR,n)
    F=matrix(RR,n)
    
    for i in range(n):
        D[i,i]=A[i,i]
        for j in range(i):
            E[i,j]=-A[i,j];
            F[j,i]=-A[j,i];
    sol=A\b;
    
    u=x0;
    res=[x0];
    resn=[norm(u-sol)];
    
    for k in range(niter):
        u=(1/omega*D-E)\(((1-omega)/omega*D+F)*u+b);
        res=res+[u];
        resn=resn+[norm(u-sol)];
        
    return resn;

Convergence logarithmique de la méthode de la relaxation pour différentes valeurs de omega

In [29]:
x0=vector(RR,[1,5,1,5,1])
n=50

resn05=relax(A,b,x0,n,0.5);
resn08=relax(A,b,x0,n,0.8);
resn12=relax(A,b,x0,n,1.2);
resn14=relax(A,b,x0,n,1.4);
resn16=relax(A,b,x0,n,1.6);
resn19=relax(A,b,x0,n,1.9);

p=point((1,log(resn05[0])),color='green')
p=point((1,log(resn08[0])),color='red')
p=point((1,log(resn12[0])),color='blue')
p=point((1,log(resn14[0])),color='cyan')
p=point((1,log(resn16[0])),color='magenta')
p=point((1,log(resn19[0])),color='black')
for i in range(n):
    p=p+point((i,log(resn05[i])),color='green')
    p=p+point((i,log(resn08[i])),color='red')
    p=p+point((i,log(resn12[i])),color='blue')
    p=p+point((i,log(resn14[i])),color='cyan')
    p=p+point((i,log(resn16[i])),color='magenta')
    p=p+point((i,log(resn19[i])),color='black')
p.show()

La meilleur valeur dans ce cas est donc omega=1.6

1-3- Méthodes du gradient

Illustration du gradient à pas constant

In [30]:
def f(x,y):
    return x^2+2*y^2+x*y+x-y

def df(x,y):
    return [2*x+y+1, x+4*y-1]

Le minimium est (-5/7,3/7), où le gradient est nul

In [31]:
import numpy as np
import matplotlib.pylab as plt

X=np.arange(-3,2,0.01)
Y=np.arange(-2,2.5,0.01)
[X,Y]=np.meshgrid(X,Y)
Z=f(X,Y)

plt.axis('equal')
plt.contour(X,Y,Z,[i*0.1 for i in srange(-100,40,5)])
plt.scatter(-5/7,3/7)
plt.show()
In [32]:
def grad(f,df,rho,X0,Y0,tol,N):
    err=1; it=1;
    X=X0; Y=Y0; Xs=[X]; Ys=[Y];
    while ((err>tol) and (it<N)):
        Xt = X - rho*df(X,Y)[0]
        Yt = Y - rho*df(X,Y)[1]
        err=abs(f(Xt,Yt)-f(X,Y))
        X=Xt; Y=Yt;
        Xs=Xs+[X]; Ys=Ys+[Y]; 
        it=it+1;
    if it==N:
        show('Algorithm stopped before reaching convergence')
    return Xs,Ys

[Xs,Ys]=grad(f,df,0.1,0,0,1E-10,100)
show([Xs[-1]+5/7,Ys[-1]-3/7])

Le résultat numérique est proche de la solution exacte

In [33]:
plt.axis('equal')
plt.contour(X,Y,Z,[i*0.1 for i in srange(-100,40,5)])
plt.scatter(Xs,Ys)
plt.show()

1-4- Méthodes de projection sur des espaces de Krylov

Illustration de la dégénérescence numérique de la base naturelle des espaces de Krylov

In [34]:
D=2*matrix([[1,0,0,0,0],[0,10,0,0,0],[0,0,100,0,0],[0,0,0,1000,0],[0,0,0,0,10000]])
P=matrix([[1,2,0,0,1],[0,1,0,-1,1],[2,0,1,0,1],[0,2,0,1,1],[1,-2,1,1,-1]])
A=P*D*P^(-1)
show('A=',A)
In [35]:
r0=vector([1,1,1,1,1])

show(r0)
show(A*r0)
show(A^2*r0)
show(A^3*r0)
show(A^4*r0)

Le nombre de chiffres significatifs de la représentation des réels en machine est dangereusement proche...

2 - Calcul de valeurs propres et de vecteurs propres

2-1- Méthodes QR

In [36]:
def MethodeQR(A,n):
    for i in range(n):
        [Q,R]=A.QR();
        A=R*Q;
    return A;
In [37]:
A=matrix(QQbar,[[-17,36,27],[-33,34,33],[-84,102,94]])
show('A=',A)
In [38]:
show('Sp(A)=',A.eigenvalues())
In [39]:
An=MethodeQR(A,5);
show('An=',An)

3- Résolution de systèmes non linéaires

Illustration du bassin de convergence de l'algorithme de Newton

On s’intéresse à la recherche des solutions complexes de l’équation z^3 = 1. Pour cela nous allons utiliser la méthode de Newton appliquée à la fonction

F : (x,y)->(Re((x+i y)^3-1),Im((x+i y)^3-1)) qui s'annule exactement lorsque (x+i y)^3=1.

Calculons explicitement la fonction F et sa différentielle DF

In [40]:
var('x y')
F=vector([-1+x^3-3*x*y^2,3*x^2*y-y^3])
DF=matrix([[3*x^2-3*y^2,-6*x*y],[6*x*y,3*x^2-3*y^2]])
show('F=',F,',  DF=',DF)

L'algorithme de Newton en 2D s'écrit de la manière suivante

In [41]:
def Newton(F,DF,x0,y0,eps):
    z=vector([x0,y0]); 
    err=[norm(DF(x=z[0],y=z[1])\F(x=z[0],y=z[1]))];
    n=0;
    while ((norm(F(x=z[0],y=z[1]))>eps) and (n<100)):
        z = z - (DF(x=z[0],y=z[1])\F(x=z[0],y=z[1]));
        err=err+[norm(DF(x=z[0],y=z[1])\F(x=z[0],y=z[1]))];
        n=n+1;
    return n,z;

On travaille sur une grille couvrant le carré [−1.5, 1.5] × [−1.5, 1.5] de pas h = 3/n .

On va résoudre F(x, y) = 0 par la méthode de Newton initialisée successivement en chacun des points (x_i , y_j ) de cette grille.

Dans un tableau, on mémorisera le numéro k entre 1 et 3 de la racine exp(2ikπ/3) vers lequel l’algorithme aura convergé.

In [42]:
n=30;
T=matrix(n,n);
l=3;
h=l/n;
racines=[(cos(0),sin(0)),(cos(2*pi/3).n(),(sin(2*pi/3)).n()),(cos(4*pi/3).n(),(sin(4*pi/3)).n())];
tol=0.001;

for i in range(n):
    xg=-l/2+h*(i-0.5);
    for j in range(n):
        yg=-l/2+h*(j-0.5);
        [iter,z]=Newton(F,DF,xg,yg,tol);
        for k in range(3):
            if ((abs(racines[k][0]-z[0])<tol) and (abs(racines[k][1]-z[1])<tol)):
                T[i,j]=k;
In [43]:
p=plot((0,0),color='red')
for i in range(n):
    xg=-l/2+h*(i-0.5);
    for j in range(n):
        yg=-l/2+h*(j-0.5);
        if (T[i,j]==0):
            p=p+point((xg,yg),color='red')
        if (T[i,j]==1):
            p=p+point((xg,yg),color='green')
        if (T[i,j]==2):
            p=p+point((xg,yg),color='yellow')

show(p)

En raffinant, on obtiendrait une figure fractale comme celle-ci !

In [44]:
from IPython.display import Image
Image("Newton3.png")
Out[44]: