從odeint scipy使用的函數提取值python

Question

我有以下腳本來使用odeint來計算dRho。

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和Q是在功能dRho內計算，它們P作用和輸入對於Q和兩個P，Q和rho_y充當出輸入。該函數返回「out」。我可以繪製出沒有任何問題，但是，我有興趣繪製P和Q。

我已經嘗試了各種方法來實現這一點：在集成方法之後重新計算P和Q，但是這增加了腳本的運行時間。因此，由於計算是在dRho內完成的，我想知道是否以及如何從外部訪問它來繪製它。

我也嘗試將P和Q連同rho0一起作爲odeint的輸入，但是P和Q在集成中被採用，導致函數返回時出現錯誤的結果。

的簡化版本：

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

在這個簡單的例子

，我想創造ydiff額外的情節。

這裏的另一個簡單的例子：

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

我感興趣的繪製所有X new和X ynew值。

又如：

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

我想對X

Answer 1

繪製MULT值以下可能是你想要的。您可以將中間值存儲在列表中，然後繪製該列表。這也需要存儲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() 

虛線是用於相同顏色的res溶液各自y_diff值。