Méthode des différences finies en Python

[HanaChan]

Bonjour,

J'ai besoin de résoudre une équation différentielle homogène autonome de second ordre :

a d²u/dx²+b du/dx+c u=0 pour x ϵ [0 , L], les conditions de bord sont u(0) = U0 et u(L) = UL.

Voici mon code en Python :

# -*-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 = f*np.ones(N)
R = e*np.ones(N-1)

A = 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 = 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)


# 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")

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

Quel est le problème ?
-----

[Calvert]

Salut,

chez moi, ça marche au moins jusqu'à N=30000, mais c'est extrêmement lent. Du coup, dur de dire quel peut-être ton problème sans avoir les messages d'erreur.

[HanaChan]

Merci pour ta réponse Il s'agit en fait d'un problème de RAM. En fait, np.diag retourne une matrice pleine. J'ai essayé spdiags du package scipy.sparse pour avoir une matrice stockée sous forme réduite et ça marche bien maintenant.

Voici le nouveau code si jamais qqn en a besoin :

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