Bonjour,

SVP pourriez-vous m'expliquer par quelle méthode la fonction numpy.linalg.solve résout un système d'équations linéaires (triangularisation puis méthode de remontée, diagonalisation, ou quoi exactement) ? 

Merci

Un coup d'oeil à la doc Numpy:

The solutions are computed using LAPACK routine _gesv.

Il semblerait que LAPACK soit écrit en Fortran. Je n'ai pas trouvé de site officiel mais on a une documentation dgesv sur un site éducatif américain:

DGESV computes the solution to a real system of linear equa-
      tions
         A * X = B, where A is an N-by-N matrix and X and B are
      N-by-NRHS matrices.

      The LU decomposition with partial pivoting and row inter-
      changes is used to factor A as
         A = P * L * U,
      where P is a permutation matrix, L is unit lower triangular,
      and U is upper triangular.  The factored form of A is then
      used to solve the system of equations A * X = B.

Bonjour,

Je vois, merci pour votre réponse

En fait, j'ai besoin de résoudre par la méthode des différences finies l'équation :

a d²u /dx² + b du/dx +c u = 0 pour x appartenant à [0,L], avec u(0) = U0 et u(L) = UL. 

Ceci à revient à résoudre l'équation matricielle A*U = B 

 

Voici ce que j'obtiens pour l'équation :  -d²u /dx² + u = 0 pour x appartenant à [0,1], avec u(0) = 2,  u(1) = 1 et N (nombre des noeuds) = 10.

La solution analytique est correcte, contrairement à la solution numérique obtenue par  numpy.linalg.solve ... 

Je ne comprends pas non plus pourquoi j'obtiens plusieurs courbes pour l'erreur ? 

Voici mon code (c'est mon tout premier code en Python) : https://drive.google.com/file/d/1PfGnmTWf7fZ5JwdR2OeL7OyHUy7WU4zX/view?usp=sharing

Pourriez-vous me dire s'il y'a une erreur quelque part ?

Merci

Il faut que tu postes ton code sur le forum (ou au moins l'essentiel) parce que je suis derrière un proxy, je ne peux pas accéder à Drive.

Tu peux utiliser MathJax pour écrire des maths dans un post.

# -*-coding:Latin-1 -*

import os

import numpy as np

# Saisie des inputs par l'utilisateur

a = input ("Saisissez la valeur de a")
a = float (a)

b = input ("Saisissez la valeur de b")
b = float (b)

c = input ("Saisissez la valeur de c")
c = float (c)

L = input ("Saisissez la valeur de L")
L = float (L)

U0 = input ("Saisissez la valeur de U0")
U0 = float (U0)

UL = input ("Saisissez la valeur de UL")
UL = float (UL)

N = input ("Saisissez la valeur de N")
N = int (N)

# Calcul du pas et des constantes e, f et g

h = L/(N+1)
e = a/(h**2)+b/h
f = c - (2*a)/(h**2)- b/h
g = a/(h**2)

# Déclaration de la matrice tridiagonale A de taille N*N 

def tridiag(P, Q, R, k1=-1, k2=0, k3=1):
    return np.diag(P, k1) + np.diag(Q, k2) + np.diag(R, k3)


P = g*np.ones(N-1)
Q = e*np.ones(N)
R = f*np.ones(N-1)

A = tridiag(P,Q,R)

# Déclaration du vecteur colonne B de taille N

B = np.zeros((N,1))
B[0,0] = -g*U0
B[N-1,0] = -e*UL


# Résolution du système linéaire A*U=B

U = np.linalg.solve(A, B)

# Vérification de la solution

np.allclose(np.dot(A, U), B)

# Résolution analytique de l'équation pour D > 0

def SolAnal1 (x, a, b, L, U0, UL, D):
	r1 = (-b+np.sqrt(D))/(2*a)
	r2 = (-b-np.sqrt(D))/(2*a)
	B = (UL-U0*np.exp(r1*L))/(np.exp(r2*L)-np.exp(r1*L))
	A = U0-B
	return A*np.exp(r1*x)+B*np.exp(r2*x)

# Résolution analytique pour D = 0

def SolAnal2 (x, a, b, L, U0, UL):
	r = -b /(2*a)
	B = U0
	A = ((UL/np.exp(r*L))- U0)/L
	return (A*x+B)*np.exp(r*x)

# Discrétisation de l'intervalle [h,Nh] en N intervalles

X = np.arange(h , L, h)
print (X)

# Calcul du Delta de l'équation caractéristique et des images de chaque noeud

D = b**2-4*a*c

if D > 0 :
	Y = SolAnal1 (X, a, b, L, U0, UL, D)
else:
	Y = SolAnal2 (X, a, b, L, U0, UL)

# Traçage des solutions numérique et analytique

import matplotlib.pyplot as plt

plt.plot(X,U,"b-.", label="Solution numérique")
plt.plot(X,Y,"r-.", label="Solution analytique")
plt.plot(X,np.abs(Y-U),"y-.", label="Erreur")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.show()


os.system("pause")

Là, perso ça dépasse ce que je peux étudier en un temps raisonnable.

Je vois. Merci en tout cas

Salut,

il y a une erreur dans l'énoncé: tu as \(g\) sur la diagonale du bas, \(f\) au milieu et \(e\) en haut. Par ailleurs, utilise plutôt

B = np.zeros(N)
B[0] = -g*U0
B[N-1] = -e*UL

Sinon le code est ok.

Nozio a écrit:

Salut,

il y a une erreur dans l'énoncé: tu as \(g\) sur la diagonale du bas, \(f\) au milieu et \(e\) en haut. Par ailleurs, utilise plutôt

B = np.zeros(N)
B[0] = -g*U0
B[N-1] = -e*UL

Sinon le code est ok.

Bonjour,

Oui, t'as raison, merci :) 

En fait, j'ai fait des essais avec a = -1, b = 0, c = 1, L = 1, U0 = 2, UL = 1, N ={10, 100, 1000} et ça marche parfaitement, mais pour N = 10000 par exemple, le programme affiche un message mais la fenêtre se ferme avant que je ne puisse lire ce qui y est écrit ... Aucun graphe n'est affiché.

Autre question, pourquoi np.allclose(np.dot(A,U),B) n'affiche pas True ?

-
Edité par HanaeAttou 17 juillet 2019 à 14:56:40

Pas de problèmes de mon côté avec 10000. Cela dit, il se peut que tu n'aies tout simplement pas assez de RAM. En fait, np.diag retourne une matrice pleine. Tu devrais peut-être essayer spdiags du package scipy.sparse pour avoir une matrice stockée sous forme réduite. Pour ce qui est du allclose, il se peut que ce soit simplement de l'erreur numérique. Pour le vérifier à la main : tu affichage A*U et B, si les deux ont l'air de se superposer, c'est des problèmes de nombre flottant et tu ne peux rien y faire. Tu peux aussi vérifier la valeur de |A*U - B|.

Je vois. Est-ce que la syntaxe suivante est correcte ?

import scipy.sparse as sp

# Déclaration de la matrice tridiagonale A de taille N*N 

P = g*np.ones(N)
Q = f*np.ones(N)
R = e*np.ones(N)

data = np.array([P, Q, R])
diags = np.array([-1, 0, 1])
A = sp.spdiags(data, diags, N, N).toarray()

Voici ce que j'obtiens pour N = 10000 (toujours pour la même équation et les mêmes conditions) :

Je ne comprends pas pourquoi la solution numérique n'est plus correcte !

-
Edité par HanaeAttou 17 juillet 2019 à 16:33:41

Mmmmm, pour commencer, vire le toarray() et utilise plutôt csc_matrix() pour la convertir, ensuite utilise spsolve() pour résoudre le problème. Je te donne mon code.

import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
 
# Saisie des inputs par l'utilisateur
a = input ("Saisissez la valeur de a")
a = float (a)
b = input ("Saisissez la valeur de b")
b = float (b)
c = input ("Saisissez la valeur de c")
c = float (c)
L = input ("Saisissez la valeur de L")
L = float (L)
U0 = input ("Saisissez la valeur de U0")
U0 = float (U0)
UL = input ("Saisissez la valeur de UL")
UL = float (UL)
N = input ("Saisissez la valeur de N")
N = int (N)

# Calcul du pas et des constantes e, f et g
h = L/(N+1)
e = a/(h**2)+b/h
f = c - (2*a)/(h**2)- b/h
g = a/(h**2)

# Déclaration de la matrice tridiagonale A de taille N*N
def tridiag(P, Q, R, k1=-1, k2=0, k3=1):
    N = len(Q)
    return (sps.spdiags(P,-1,N,N) + sps.spdiags(Q,0,N,N) + sps.spdiags(R,1,N,N))

P = g*np.ones(N)
Q = f*np.ones(N)
R = e*np.ones(N)
 
A = sps.csc_matrix(tridiag(P,Q,R))

# Déclaration du vecteur colonne B de taille N
B = np.zeros(N)
B[0] = -g*U0
B[N-1] = -e*UL
 
 
# Résolution du système linéaire A*U=B
U = spsolve(A,B)

# Résolution analytique de l'équation pour D > 0
 
def SolAnal1 (x, a, b, L, U0, UL, D):
    r1 = (-b+np.sqrt(D))/(2*a)
    r2 = (-b-np.sqrt(D))/(2*a)
    B = (UL-U0*np.exp(r1*L))/(np.exp(r2*L)-np.exp(r1*L))
    A = U0-B
    return A*np.exp(r1*x)+B*np.exp(r2*x)
 
# Résolution analytique pour D = 0
 
def SolAnal2 (x, a, b, L, U0, UL):
    r = -b /(2*a)
    B = U0
    A = ((UL/np.exp(r*L))- U0)/L
    return (A*x+B)*np.exp(r*x)
 
# Discrétisation de l'intervalle [h,Nh] en N intervalles
X = np.arange(h , L, h)
 
# Calcul du Delta de l'équation caractéristique et des images de chaque noeud
D = b**2-4*a*c
if D > 0 :
    Y = SolAnal1 (X, a, b, L, U0, UL, D)
else:
    Y = SolAnal2 (X, a, b, L, U0, UL)
 
# Traçage des solutions numérique et analytique
import matplotlib.pyplot as plt
plt.plot(X,U,"b-+", label="Solution numérique")
plt.plot(X,Y,"r-.", label="Solution analytique")
plt.plot(X,np.abs(Y-U),"y-.", label="Erreur")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.show()

Merciiiiiii infiniment ça marche parfaitement maintenant.

J'ai encore deux petites questions:

- La ligne 29 N = len(Q) sert à quoi exactement ? 

- N'y aurait-il pas moyen d'inclure les 2 points de coordonnées (0,U0) et (L,UL) dans les courbes ? Je pensais à une fonction qui permettrait d'insérer U0 avant le 1er élément de U et UL après le dernier élément, ça existe ça?

-
Edité par HanaeAttou 17 juillet 2019 à 18:16:49

c'est juste pour récupérer la taille du problème. C'est une facilité qui ne coûte rien en terme de temps mais ça surcharge un peu moins. Pour ce qui est de rajouter U0 et Ul, oui c'est possible avec numpy.intert() si tu connais la position à laquelle tu veux ajouter une valeur, ou numpy.append() pour rajoute une valeur à la fin. Je te renvoie à la doc pour le détail du fonctionnement et les exemples.

J'ai essayé la syntaxe suivante :

# Insertion de U0 au début de U et de UL à la fin de U
np.insert (U,0,[U0])
np.append (U,[UL])

J'ai aussi inclus 0 et L à X :

# Discrétisation de l'intervalle [0,L] en N+1 intervalles
X = np.linspace(0.0 , L, N+2)

mais ça ne marche pas, le programme affiche les valeurs de U (U0 et UL n'ont pas été rajoutés) et X suivies d'autres lignes que je n'arrive pas à lire parce qu'il se ferme vite ...

La réponse est dans la doc Comme beaucoup de fonctions, numpy.insert() a des entrées ... et des sorties !

Ok, ça marche, merci :) 

-
Edité par HanaeAttou 19 juillet 2019 à 16:43:16

Petite question concernant sps.spadiags(P,-1,N,N) : à ce que j'ai compris, cette fonction crée une matrice de taille N*N où les éléments de P sont insérés dans la 1ère diagonale inférieure de la matrice créée. Les éléments de P sont choisis de quel indice jusqu'à quel indice ? Est-ce qu'on peut modifier ce choix ?