从odeint scipy python使用的函数中提取值

如何解决从odeint scipy python使用的函数中提取值

以下可能是您想要的。您可以将中间值存储在列表中，然后再绘制该列表。那也需要存储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()

虚线是相同颜色y_diff的res溶液的相应值。

解决方法

我有以下脚本使用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的输入。该函数返回“
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],1])
plt.show()

我对绘制所有xnew和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,1000)
y0 = np.zeros(4)+1.25
res = odeint(dY,x)
dif = np.array([dY(y,x) for y in res])
print(dif)
plt.plot(x,'-')
plt.show()

我想针对x绘制多值