Résolution d'une équation aux dérivées partielles en Python

[invitea45e89ab]

Bonjour,

J'ai besoin de résoudre une équation aux dérivées partielle par la méthode des différences finies : https://drive.google.com/file/d/1QQt...ew?usp=sharing

Je choisis comme a(x) la fonction carré.

Voici mon code :

import os
import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
from scipy import misc
*
# Saisie des inputs par l'utilisateur
c = input ("Saisissez la valeur de c")
c = float (c)
e = input ("Saisissez la valeur de epsilon")
e = float (e)
L = input ("Saisissez la longeur de l'intervalle [0,L]")
L = float (L)
T = input ("Saisissez la longeur de l'intervalle [0,T]")
T = float (T)
U0t = input ("Saisissez la valeur de u(0,t)")
U0t = float (U0t)
ULt = input ("Saisissez la valeur de u(L,t)")
ULt = float (ULt)
Ux0 = input ("Saisissez la valeur de u(x,0)")
Ux0 = float (Ux0)
N = input ("Saisissez le nombre des noeuds de l'espace")
N = int (N)
M = input ("Saisissez le nombre des noeuds de du temps")
M = int (M)
*
# Calcul du pas de l'espace
h = L/(N+1)
*
# Discrétisation de l'intervalle [0,L] selon le pas h
X = np.linspace (0.0,L,N+2)
*
*
# Calcul du pas du temps
t = T/(M+1)
*
# Discrétisation de l'intervalle [0,T] selon le pas t
Y = np.linspace (0.0,T,M+2)
*
*
# Fonction a(x)
def a(x):
****return(x**2)
*
# Remplissage 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 = [- misc.derivative(a,(i)*h)/(2*h) + a((i)*h)/h**2 for i* in range(2,N+2)]
Q = [c - 2*a((i+1)*h)/h**2 for i in range(N)]
R = [misc.derivative(a,(i)*h)/(2*h) + a((i)*h)/h**2 for i in range(N)]
*
A = sps.csc_matrix(tridiag(P,Q,R))
print (A)
*
# Remplissage du vecteur colonne B de taille N
B = np.zeros(N)
B[0] = (- misc.derivative(a,h)/(2*h) + a(h)/h**2)*U0t
B[N-1] = (misc.derivative(a,N*h)/(2*h) + a(N*h)/h**2)*ULt
print (B)
*
# Initialisation de Uj à l'instant tj=0
U = Ux0*np.ones(N)
print(U)
*
# Déclaration de la matrice identité de taille N*N
I = np.eye(N)
*
# Nouvelle A
A = (t/e)*A + I
*
# Nouveau B
B = (t/e)*B
*
# Calcul de Uj+1 à partir de Uj
*
for j in range (1,M+2):
****U = np.dot(A,U)+ B
****print(U)
*
os.system("pause")

Le problème réside dans la boucle For (ou plus précisément dans le calcul de Uj+1 à partir de Uj).

Lorsque j'ai fait un essai pour c = -1, epsilon = 1, longueur de l'intervalle [0,L] = 1, longueur de l'intervalle [0,T] = 1, u(0,t)= 2, u(L,t) = 1, u(x,0) = 1.5, Nombre des nœuds de l'espace N = 2, Nombre des nœuds du temps M = 1, le programme demande tous les inputs puis se ferme rapidement.

Lorsque j'ai fait un autre essai (toujours les mêmes inputs) mais sans boucle (je voulais juste calculer U1 à l'instant tj = t1), le programme marche bien et affiche la valeur correcte de U1. Mais lorsque je retape le même code pour calculer U2 à partir de U1, le programme demande juste les inputs puis se ferme !

Je ne comprends pas non plus pourquoi la boucle n'affiche même pas la première itération sachant qu'elle est correcte.

Merci de votre aide
-----

[invite6c250b59]

Est-ce qu'il n'y a pas un pb avec les t vs T et L vs h? (cf ci-dessous une tentative pour une version un peu plus lisible)

import os
import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
from scipy.misc import derivatives as drv

t1, t2, t3 = "Saisissez la valeur de ", "Longueur de l'intervalle ", "Nombre de noeuds "
c = float(input(t1 + "c"))
e = float(input(t1 + "epsilon")                                 
L = float(input(t2 + "[0,L]?")
T = float(input(t2 + "[0,T]?")
U0t = float(input(t1 + "u(0,t)")
ULt = float(input(t1 + "u(L,t)")
Ux0 = float(input(t1 + "u(x,0)")
N = int(input(t3 + "de l'espace")
M =  int(input(t3 + "du temps")

h = L/(N+1)                                                
X = np.linspace (0.0,L,N+2)                         # ICI L ou h?
t = T/(M+1)                                                
Y = np.linspace (0.0,T,M+2)                         # ICI T ou t?

def a(x):                                                     
     return(x**2)

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

d1, d2, P, Q, R = 2*h, h**2, [], [], []
for i in range(N): 
    P.append( a( (i+2)*h ) / d2 - drv(a, (i+2)*h) / d1 )
    Q.append( c - 2 * a( (i+1)*h ) / d2 )
    R.append( a( i*h ) / d2 + drv(a, i*h) / d1 )

A = sps.csc_matrix(tridiag(P,Q,R))
print (A)

B = np.zeros(N)                                         
B[0] = (- drv(a,h)/d1 + a(h)/d2)*U0t
B[N-1] = (drv(a,N*h)/d1 + a(N*h)/d2)*ULt
print (B)

U = Ux0*np.ones(N)                                   # Initialisation de Uj à l'instant tj=0
print(U)

A, B = (t/e)*A + np.eye(N), (t/e)*B                                         # update A, B

for j in range (1,M+2):                                # Uj+1 from Uj
    U = np.dot(A,U)+ B
    print(U)

os.system("pause")

[invitea45e89ab]

Merci pour votre réponse En fait, non, le problème n'est pas là.

[invitea45e89ab]

Je viens de découvrir l'erreur: il s'agissait d'une incompatibilité entre les tailles de A, U et B.

J'ai modifié le code au niveau de la boucle For comme suivant et ça marche bien maintenant:

# Calcul de Uj+1 à partir de Uj

for j in range (1,M+2):
	U = np.dot(A,U)+ B
	print(U)
	U = np.reshape (U,(N,1))
	B = np.reshape (B,(N,1))

[A voir en vidéo sur Futura]