extrait les valeurs de la fonction utilisée par odeint scipy python

Question

J'ai le script suivant pour calculer dRho en utilisant odeint.

P_r = 10e5 
rho_r = 900 
L = 750 
H = 10 
W = 150 
A = H * W 
V = A * L 
fi = 0.17 

k = 1.2e-13 
c = 12.8e-9 
mu = 2e-3 

N = 50 
dV = V/N 
dx = L/N 

P_in = P_r 
rho_in = rho_r 

P_w = 1e5  
rho_w = rho_r* np.exp(c*(P_w-P_r)) 

# init initial case 
P = np.empty(N+1)*10e5 
Q = np.ones(N+1) 
out = np.empty(N+1) 

P[0] = P_w 
Q[0] = 0 
out[0] = 0 

def dRho(rho_y, t, N): 

    P[1:N] = P_r + (1/c) * np.log(rho_y[1:N]/rho_r) 
    P[N] = P_r + (1/c) * np.log(rho_y[N]/rho_r) 


    Q[1:N] = (-A*k/mu)*((P[1-1:N-1] - P[1:N])/dx) 
    Q[N] = (-A*k/mu)*((P[N]-P_r)/dx) 


    out[1:N] = ((Q[1+1:N+1]*rho_y[1+1:N+1] - Q[1:N]*rho_y[1:N])/dV) 
    out[N] = 0 

    return out 

t0 = np.linspace(0,1e9, int(1e9/200)) 
rho0 = np.ones(N+1)*900 
ti = time.time() 
solve = odeint(dRho, rho0, t0, args=(N,)) 
plt.plot(t0,solve[:,1:len(rho0)], '-', label='dRho') 
plt.legend(loc='upper right') 
plt.show() 

P et Q sont calculées dans la fonction dRho, ils P agit et d'entrée pour les deux Q et P, Q et agir comme entrée pour rho_y out. La fonction renvoie "out". Je peux tracer sans aucun problème, cependant, je suis intéressé à tracer P et Q ainsi.

J'ai essayé différentes approches pour y parvenir comme: recalculer P et Q après la méthode d'intégration mais cela a augmenté le temps d'exécution du script. Donc, puisque le calcul est fait dans dHho je me demandais si et comment je pouvais y accéder de l'extérieur pour le tracer.

J'ai également essayé d'ajouter P et Q ensemble avec rho0 comme entrée pour odeint mais les deux P et Q ont été pris dans l'intégration qui a entraîné un mauvais résultat quand retourné par la fonction.

Une version simplifiée:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.integrate import odeint 
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] 
    print(y) 
    return (a/dC)*y_diff+b1*dC 

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

dans cet exemple simplifié, je voudrais créer une parcelle supplémentaire de ydiff.

Voici un autre cas simple:

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

def func(z,t): 
    x, y=z 
    xnew = x*2 
    print(xnew) 
    ynew = y*0.5 
#  print y 
    return [x, y]  

z0=[1,3] 
t = np.linspace(0,10) 
xx=odeint(func, z0, t) 
plt.plot(t, xx[:,0],t,xx[:,1]) 
plt.show() 

Je suis intéressé à tracer toutes les valeurs xnew et YNouveau.

Un autre exemple:

xarr = np.ones(4) 
def dY(y, x): 
    a = 0.001 
    yin = 1 
    C = 0.01 
    N = 1 
    dC = C/N 
    b1 = 0 
    xarr[0] = 0.25 
    xarr[1:] = 2 
    mult = xarr*2 
    out = mult * y 
    print(mult) 
    return out 

x = np.linspace(0,20,1000) 
y0 = np.zeros(4)+1.25 
res = odeint(dY, y0, x) 
dif = np.array([dY(y,x) for y in res]) 
print(dif) 
plt.plot(x,res, '-') 
plt.show() 

Je voudrais tracer les valeurs mult contre x

Answer 1

Le suivant pourrait être ce que vous voulez. Vous pouvez stocker les valeurs intermédiaires dans une liste et ensuite tracer cette liste. Cela nécessiterait également de stocker les valeurs x.

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

xs = [] 
yd = [] 

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] 
    xs.append(x) 
    yd.append(y_diff) 
    return (a/dC)*y_diff+b1*dC 

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

plt.plot(x,res, '-') 

plt.gca().set_prop_cycle(plt.rcParams['axes.prop_cycle']) 
plt.plot(np.array(xs),np.array(yd), '-.') 

plt.show() 

lignes en pointillé sont les valeurs respectives pour les y_diffres solutions de la même couleur.