Scipy: Deux façons d'implémenter une équation différentielle: deux solutions différentes:

Question

J'ai essayé de résoudre une équation différentielle pour ma thèse en chimie et là je suis tombé sur une question concernant le solveur d'équation différentielle "odeint" de scipy.

D'abord j'ai implémenté le différentiel par la fonction CIDNP_1 (CIDNP est un phénomène chimique, ce qui explique les variables inhabituelles) selon l'exemple sur le site scipy. Mais la solution est loin de la bonne direction.

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.integrate 

R0 = 5e+5 
kt = 5e5/R0 
beta = 3/R0 

def CIDNP_1(y, t): 
    dP_dt, dQ_dt = y 

    def R(t): 
     return R0/(1 + kt*R0*t) 

    dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 
    dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 
    return [dP_dt, dQ_dt] 


def CIDNP_2(y, t):  
    dP_dt, dQ_dt = y 

    def R(t): 
     return R0/(1 + kt*R0*t) 

    return [-kt*dP_dt*R(t) - kt*beta*(R(t))**2, \ 
      +kt*dP_dt*R(t) + kt*beta*(R(t))**2] 

y0 = [-1, +1] 
t = np.linspace(1e-9, 100e-6, 1e3) 
sol_1 = scipy.integrate.odeint(CIDNP_1, y0, t) 
sol_2 = scipy.integrate.odeint(CIDNP_2, y0, t) 

Puis j'ai changé ma solution à CIDNP_2 ce qui a donné le résultat correct, mais à mes yeux il n'y a pas de différence dans la mise en œuvre que les variables dP_dt et dQ_dt ne sont pas modifiés dans la mise en œuvre CIDNP_1.

Alors comme tout le monde peut me donner un indice pourquoi la mise en œuvre CIDNP_1 donne le mauvais résultat, je serais très chanceux car au moins les deux dernières heures n'étaient pas complètement perdues.

Cordialement,

Jakob

Answer 1

Dans votre première version que vous effectuez les mises à jour pas en même temps, que vous exécutez la ligne à

dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 

pas simulatnously; par conséquent, vous utilisez le dP_dt déjà mis à jour pour mettre à jour dQ_dt. Ceci est une implémentation défectueuse du système ODE. Votre deuxième approche est meilleure, car elle n'a pas ce genre d'erreur. Vous devez soit renvoyer directement les valeurs mises à jour, soit enregistrer les nouvelles valeurs calculées de dP_dt dans une autre variable avant de calculer le nouveau dQ_dt.

new_dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 
new_dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 

return [new_dP_dt, new_dQ_dt] 

cela résoudrait votre problème.

Answer 2

En CIDNP_1 vous modifiez la valeur de dP_dt avant d'utiliser la nouvelle valeur pour le calcul dQ_dt:

dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 # changed dP_dt! 
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 # use the changed value! 

Dans CIDNP_2 vous les calculer simultanément, à savoir dQ_dt est calculée avec la valeur d'origine de dP_dt, pas celui changé . Vous pouvez penser comme

a = -kt*dP_dt*R(t) - kt*beta*(R(t))**2  # no overwriting - 
b = +kt*dP_dt*R(t) + kt*beta*(R(t))**2  # uses original value of dP_dt 
return [a, b] 

Vous pouvez aussi accélérer probablement les choses un peu comme

def CIDNP_3(y, t): 
    dP_dt, dQ_dt = y 
    R_t = R0/(1 + kt * R0 * t) 
    res = kt * R_t * (dP_dt + beta * R_t) 
    return [-res, res]