implémenter une équation mathématique d'intégration en utilisant odeint dans Python

Question

J'essaie de résoudre l'équation suivante en python en utilisant la fonction scipy.odeint.

Actuellement, je suis en mesure de mettre en œuvre cette forme de l'équation

en python utilisant le script suivant:

def dY(y1, x): 
    a = 0.001 
    yin = 1 
    C = 0.01 
    N = 1 
    dC = C/N 
    b1 = 0 
    return (a/dC)*(yin-y1)+b1*dC 

x = np.linspace(0,20,1000) 
y0 = 0 
res = odeint(dY, y0, x) 
plt.plot(t,res, '-') 
plt.show() 

Mon problème avec la première équation est 'je'. Je ne sais pas comment intégrer l'équation et être toujours capable de fournir les valeurs y et yi (yi-1 et yi) actuelles et précédentes. 'i' est simplement un numéro de séquence compris dans une plage de 0..100.

Edit 1:

L'équation d'origine est:

Ce que je réécrit en utilisant Y, X, a, b et C

Edit2: I édité code « Pierre de Buyl et changé la valeur de N. Heureusement, j'ai une table de validation pour valider le résultat. Malheureusement, les résultats ne sont pas égaux.

Voici ma table de validation:

et est ici la sortie numpy:

Code occasion:

def dY(y, x): 
    a = 0.001 
    yin = 1 
    C = 0.01 
    N = 3 
    dC = C/N 
    b1 = 0.01 
    y_diff = -np.copy(y) 
    y_diff[0] += yin 
    y_diff[1:] += y[:-1] 
    return (a/dC)*(y_diff)+b1*dC 

x = np.linspace(0,20,11) 
y0 = np.zeros(3) 
res = odeint(dY, y0, x) 
plt.plot(x,res, '-') 

que vous pouvez voir la les valeurs sont différentes d'un n décalage de 0,02 ..

Ai-je manqué quelque chose qui entraîne ce décalage?

Answer 1

L'équation est un « couplé » équation différentielle ordinaire (voir « Système de EDO » sur Wikipedia.

La variable est un vecteur contenant y[0], y[1], etc. Pour résoudre le ODE vous devez nourrir un vecteur . la condition initiale et la fonction dY doit retourner un vecteur aussi bien

j'ai modifié votre code pour obtenir ce résultat:

def dY(y, x): 
    a = 0.001 
    yin = 1 
    C = 0.01 
    N = 1 
    dC = C/N 
    b1 = 0 
    y_diff = -np.copy(y) 
    y_diff[0] += yin 
    y_diff[1:] += y[:-1] 
    return (a/dC)*y_diff+b1*dC 

Je WRI La partie y[i-1] - y[i] est une opération vectorielle NumPy et une case spéciale est la coordonnée y[0] (c'est le y1 dans votre notation mais les tableaux commencent à 0 en Python).

La solution, en utilisant une valeur initiale de 0 pour tous yi est

x = np.linspace(0,20,1000) 
y0 = np.zeros(4) 
res = odeint(dY, y0, x) 
plt.plot(x,res, '-')