从 odeint scipy python 使用的函数中提取值
extract values from function used by odeint scipy python
我有以下脚本来使用 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的输入。函数 returns "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()
我有兴趣绘制所有 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,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
绘制多个值
以下可能是您想要的。您可以将中间值存储在列表中,然后绘制该列表。这也需要存储 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
值。
我有以下脚本来使用 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的输入。函数 returns "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()
我有兴趣绘制所有 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,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
绘制多个值以下可能是您想要的。您可以将中间值存储在列表中,然后绘制该列表。这也需要存储 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
值。