Partage
  • Partager sur Facebook
  • Partager sur Twitter

Fonction numpy.linalg.solve

Sujet résolu
    15 juillet 2019 à 16:53:35

    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 :)

    • Partager sur Facebook
    • Partager sur Twitter
      15 juillet 2019 à 17:10:24

      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.
      • Partager sur Facebook
      • Partager sur Twitter
        15 juillet 2019 à 17:47:19

        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 :)

        • Partager sur Facebook
        • Partager sur Twitter
          16 juillet 2019 à 9:34:09

          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.

          • Partager sur Facebook
          • Partager sur Twitter
            16 juillet 2019 à 15:23:38

            # -*-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")
            • Partager sur Facebook
            • Partager sur Twitter
              16 juillet 2019 à 16:03:03

              Là, perso ça dépasse ce que je peux étudier en un temps raisonnable.
              • Partager sur Facebook
              • Partager sur Twitter
                16 juillet 2019 à 16:25:32

                Je vois. Merci en tout cas :)
                • Partager sur Facebook
                • Partager sur Twitter
                  16 juillet 2019 à 17:41:05

                  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.

                  • Partager sur Facebook
                  • Partager sur Twitter

                  Avez-vous entendu parler de Julia ? Laissez-vous tenter ...

                    17 juillet 2019 à 14:48:24

                    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

                    • Partager sur Facebook
                    • Partager sur Twitter
                      17 juillet 2019 à 15:29:45

                      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|.
                      • Partager sur Facebook
                      • Partager sur Twitter

                      Avez-vous entendu parler de Julia ? Laissez-vous tenter ...

                        17 juillet 2019 à 16:17:37

                        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

                        • Partager sur Facebook
                        • Partager sur Twitter
                          17 juillet 2019 à 17:07:48

                          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()
                          • Partager sur Facebook
                          • Partager sur Twitter

                          Avez-vous entendu parler de Julia ? Laissez-vous tenter ...

                            17 juillet 2019 à 17:38:26

                            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

                            • Partager sur Facebook
                            • Partager sur Twitter
                              18 juillet 2019 à 12:50:33

                              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.
                              • Partager sur Facebook
                              • Partager sur Twitter

                              Avez-vous entendu parler de Julia ? Laissez-vous tenter ...

                                18 juillet 2019 à 15:59:48

                                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 ...


                                • Partager sur Facebook
                                • Partager sur Twitter
                                  18 juillet 2019 à 17:51:06

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

                                  Avez-vous entendu parler de Julia ? Laissez-vous tenter ...

                                    19 juillet 2019 à 16:05:45

                                    Ok, ça marche, merci :) 

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

                                    • Partager sur Facebook
                                    • Partager sur Twitter
                                      21 juillet 2019 à 20:20:11

                                      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 ?
                                      • Partager sur Facebook
                                      • Partager sur Twitter

                                      Fonction numpy.linalg.solve

                                      × Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
                                      • Editeur
                                      • Markdown