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) ?
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.
# -*-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")
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 ?
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|.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
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 !
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()
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
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?
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.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
# 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 ...
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 ?
Fonction numpy.linalg.solve
× Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
× Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...