Python scipy solve_ivp / odeint：打印和绘制 time/solution-dependent 参数
Python scipy solve_ivp / odeint: print & plot time/solution-dependent parameters
我一直在使用 scipy.integrate.solve_ivp 求解二阶 ODE 系统。
假设代码在这里可用（注意这使用 odeint，但与 solve_ivp 的方法类似）：https://scipy-cookbook.readthedocs.io/items/CoupledSpringMassSystem.html
现在，使一些参数（b1 和 b2）依赖于时间和解决方案，例如：
b1 = 0.8 * x1
if (x2 >= 1.0):
b2 = 20*x2
else:
b2 = 1.0
solve_ivp 的解决方案可以很容易地打印和绘制。
我还可以通过在定义的函数中插入 print(b1) 来打印部分 time/solution-dependent 参数。但我对这个输出感到困惑。我所要做的是打印并绘制在解决方案时间戳处解决的已定义函数中设置的 time/solution-dependent 参数：例如print(t, b1).
希望这是有道理的。任何提示都会很棒。
干杯
示例代码（使用 odeint 和 b1 & b2 作为动态参数）：
def vectorfield(w, t, p):
"""
Defines the differential equations for the coupled spring-mass system.
Arguments:
w : vector of the state variables:
w = [x1,y1,x2,y2]
t : time
p : vector of the parameters:
p = [m1,m2,k1,k2,L1,L2,b1,b2]
"""
x1, y1, x2, y2 = w
m1, m2, k1, k2, L1, L2, = p
# Friction coefficients
b1 = 0.8 * x1
if (x2 >= 1.0):
b2 = 20*x2
else:
b2 = 1.0
# Create f = (x1',y1',x2',y2'):
f = [y1,
(-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
y2,
(-b2 * y2 - k2 * (x2 - x1 - L2)) / m2]
print('b1', b1)
# print('b2', b2)
return f
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint
# Parameter values
# Masses:
m1 = 1.0
m2 = 1.5
# Spring constants
k1 = 8.0
k2 = 40.0
# Natural lengths
L1 = 0.5
L2 = 1.0
# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial velocities
x1 = 0.5
y1 = 0.0
x2 = 2.25
y2 = 0.0
# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 10.0
numpoints = 101
# Create the time samples for the output of the ODE solver.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, L1, L2]
w0 = [x1, y1, x2, y2]
# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
atol=abserr, rtol=relerr)
# Print Output
print(t)
# print(wsol)
# print(t, b1)
# print(t, b2)
# Plot the solution that was generated
from numpy import loadtxt
from pylab import figure, plot, xlabel, grid, legend, title, savefig
from matplotlib.font_manager import FontProperties
figure(1, figsize=(6, 4.5))
xlabel('t')
grid(True)
lw = 1
plot(t, wsol[:,0], 'b', linewidth=lw)
plot(t, wsol[:,1], 'g', linewidth=lw)
plot(t, wsol[:,2], 'r', linewidth=lw)
plot(t, wsol[:,3], 'c', linewidth=lw)
# plot(t, b1, 'k', linewidth=lw)
legend((r'$x_1$', r'$x_2$',r'$x_3$', r'$x_4$'), prop=FontProperties(size=16))
title('Mass Displacements for the\nCoupled Spring-Mass System')
# savefig('two_springs.png', dpi=100)
新的 ODE 求解器 API
调整您的代码以使用 solve_ivp（这是现代的 scipy API 来求解 ODE）我们可以以非侵入式的方式求解系统：
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
numpoints = 1001
t0 = np.linspace(0, stoptime, numpoints)
def func(t, x0, beta):
return vectorfield(x0, t, beta)
#sol1 = odeint(vectorfield, w0, t0, args=(p,), atol=abserr, rtol=relerr)
sol2 = solve_ivp(func, [t0[0], t0[-1]], w0, args=(p,), t_eval=t0, method='LSODA')
新的 API returns 与 odeint 的结果非常相似。 solve_ivp 返回的解决方案对象提供了有关作业执行的多个额外信息：
message: 'The solver successfully reached the end of the integration interval.'
nfev: 553
njev: 29
nlu: 29
sol: None
status: 0
success: True
t: array([ 0. , 0.01, 0.02, ..., 9.98, 9.99, 10. ])
t_events: None
y: array([[ 0.5 , 0.50149705, 0.50596959, ..., 0.60389634,
0.60370222, 0.6035087 ],
[ 0. , 0.29903815, 0.59476327, ..., -0.01944383,
-0.01937904, -0.01932493],
[ 2.25 , 2.24909287, 2.24669965, ..., 1.62461438,
1.62435635, 1.62409906],
[ 0. , -0.17257692, -0.29938513, ..., -0.02583719,
-0.02576507, -0.02569245]])
y_events: None
捕获动态参数
如果我正确理解你的问题，你也想将你的摩擦参数作为时间序列。
跟随求解器执行
这可以使用 callback function 或在评估 ODE 函数时传递额外的数据结构和存储系数值来实现。
如果您想跟踪求解器的执行，这将很有帮助，但结果将取决于求解器积分网格而不是所需的求解时间尺度。
无论如何，如果你想这样做，只需调整你的功能：
def vectorfield2(w, t, p, store):
# ...
store.append({'time': t, 'coefs': [b1, b2]})
# ...
c = []
sol1 = odeint(vectorfield2, w0, t0, args=(p, c), atol=abserr, rtol=relerr)
如您所见，结果不遵循求解时间尺度，而是求解器积分网格：
[{'time': 0.0, 'coefs': [0.4, 45.0]},
{'time': 3.331483023263848e-07, 'coefs': [0.4, 45.0]},
{'time': 3.331483023263848e-07,
'coefs': [0.4000000000026638, 44.999999999955605]},
{'time': 6.662966046527696e-07,
'coefs': [0.4000000000053275, 44.99999999991122]},
{'time': 6.662966046527696e-07,
'coefs': [0.40000000000799124, 44.99999999986683]},
{'time': 0.0001385450759485236,
'coefs': [0.4000002303394594, 44.99999616108455]}, ...]
作为状态函数的系数
另一方面，您的系数仅取决于系统状态，这意味着我们可以在 ODE 函数内部计算它（每次求解器调用它时），然后在时间相关的解决方案已知时计算它。
我建议按如下方式重写您的系统：
def coefs(x):
# Ensure signature:
if isinstance(x, (list, tuple)):
x = np.array(x)
if len(x.shape) == 1:
x = x.reshape(1, -1)
# Compute coefficients
b1 = 0.8 * x[:,0]
b2 = 20. * x[:,2]
q2 = (x[:,2] < 1.0)
b2[q2] = 1.0
return np.stack([b1, b2]).T
def system(t, w, p):
x1, y1, x2, y2 = w
m1, m2, k1, k2, L1, L2, = p
b = coefs(w)
return [
y1, (-b[:,0] * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
y2, (-b[:,1] * y2 - k2 * (x2 - x1 - L2)) / m2,
]
sol3 = solve_ivp(system, [t0[0], t0[-1]], w0, args=(p,), t_eval=t0, method='LSODA')
请注意，coefs 方法的接口旨在对简单状态（列表或向量）和时间解决方案（矩阵）进行操作。
然后，很容易显示时间的解决方案：
fig, axe = plt.subplots()
axe.plot(t0, sol3.y.T)
axe.set_title("Dynamic System: Coupled Spring Mass")
axe.set_xlabel("Time, $t$")
axe.set_ylabel("System Coordinates, $x_i(t)$")
axe.legend([r'$x_%d$' % i for i in range(sol3.y.T.shape[1])])
axe.grid()
和系数：
fig, axes = plt.subplots(2,1, sharex=True, sharey=False)
axes[0].plot(t0, coefs(sol3.y.T)[:,0])
axes[1].plot(t0, coefs(sol3.y.T)[:,1])
axes[0].set_title("Dynamic System: Coupled Spring Mass")
axes[1].set_xlabel("Time, $t$")
axes[0].set_ylabel("$b_0(t)$")
axes[1].set_ylabel("$b_1(t)$")
for i in range(2):
axes[i].grid()
这似乎更有可能是您想要实现的目标。
我一直在使用 scipy.integrate.solve_ivp 求解二阶 ODE 系统。
假设代码在这里可用（注意这使用 odeint，但与 solve_ivp 的方法类似）：https://scipy-cookbook.readthedocs.io/items/CoupledSpringMassSystem.html
现在，使一些参数（b1 和 b2）依赖于时间和解决方案，例如：
b1 = 0.8 * x1
if (x2 >= 1.0):
b2 = 20*x2
else:
b2 = 1.0
solve_ivp 的解决方案可以很容易地打印和绘制。
我还可以通过在定义的函数中插入 print(b1) 来打印部分 time/solution-dependent 参数。但我对这个输出感到困惑。我所要做的是打印并绘制在解决方案时间戳处解决的已定义函数中设置的 time/solution-dependent 参数：例如print(t, b1).
希望这是有道理的。任何提示都会很棒。
干杯
示例代码（使用 odeint 和 b1 & b2 作为动态参数）：
def vectorfield(w, t, p):
"""
Defines the differential equations for the coupled spring-mass system.
Arguments:
w : vector of the state variables:
w = [x1,y1,x2,y2]
t : time
p : vector of the parameters:
p = [m1,m2,k1,k2,L1,L2,b1,b2]
"""
x1, y1, x2, y2 = w
m1, m2, k1, k2, L1, L2, = p
# Friction coefficients
b1 = 0.8 * x1
if (x2 >= 1.0):
b2 = 20*x2
else:
b2 = 1.0
# Create f = (x1',y1',x2',y2'):
f = [y1,
(-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
y2,
(-b2 * y2 - k2 * (x2 - x1 - L2)) / m2]
print('b1', b1)
# print('b2', b2)
return f
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint
# Parameter values
# Masses:
m1 = 1.0
m2 = 1.5
# Spring constants
k1 = 8.0
k2 = 40.0
# Natural lengths
L1 = 0.5
L2 = 1.0
# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial velocities
x1 = 0.5
y1 = 0.0
x2 = 2.25
y2 = 0.0
# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 10.0
numpoints = 101
# Create the time samples for the output of the ODE solver.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, L1, L2]
w0 = [x1, y1, x2, y2]
# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
atol=abserr, rtol=relerr)
# Print Output
print(t)
# print(wsol)
# print(t, b1)
# print(t, b2)
# Plot the solution that was generated
from numpy import loadtxt
from pylab import figure, plot, xlabel, grid, legend, title, savefig
from matplotlib.font_manager import FontProperties
figure(1, figsize=(6, 4.5))
xlabel('t')
grid(True)
lw = 1
plot(t, wsol[:,0], 'b', linewidth=lw)
plot(t, wsol[:,1], 'g', linewidth=lw)
plot(t, wsol[:,2], 'r', linewidth=lw)
plot(t, wsol[:,3], 'c', linewidth=lw)
# plot(t, b1, 'k', linewidth=lw)
legend((r'$x_1$', r'$x_2$',r'$x_3$', r'$x_4$'), prop=FontProperties(size=16))
title('Mass Displacements for the\nCoupled Spring-Mass System')
# savefig('two_springs.png', dpi=100)
新的 ODE 求解器 API
调整您的代码以使用 solve_ivp（这是现代的 scipy API 来求解 ODE）我们可以以非侵入式的方式求解系统：
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
numpoints = 1001
t0 = np.linspace(0, stoptime, numpoints)
def func(t, x0, beta):
return vectorfield(x0, t, beta)
#sol1 = odeint(vectorfield, w0, t0, args=(p,), atol=abserr, rtol=relerr)
sol2 = solve_ivp(func, [t0[0], t0[-1]], w0, args=(p,), t_eval=t0, method='LSODA')
新的 API returns 与 odeint 的结果非常相似。 solve_ivp 返回的解决方案对象提供了有关作业执行的多个额外信息：
message: 'The solver successfully reached the end of the integration interval.'
nfev: 553
njev: 29
nlu: 29
sol: None
status: 0
success: True
t: array([ 0. , 0.01, 0.02, ..., 9.98, 9.99, 10. ])
t_events: None
y: array([[ 0.5 , 0.50149705, 0.50596959, ..., 0.60389634,
0.60370222, 0.6035087 ],
[ 0. , 0.29903815, 0.59476327, ..., -0.01944383,
-0.01937904, -0.01932493],
[ 2.25 , 2.24909287, 2.24669965, ..., 1.62461438,
1.62435635, 1.62409906],
[ 0. , -0.17257692, -0.29938513, ..., -0.02583719,
-0.02576507, -0.02569245]])
y_events: None
捕获动态参数
如果我正确理解你的问题，你也想将你的摩擦参数作为时间序列。
跟随求解器执行
这可以使用 callback function 或在评估 ODE 函数时传递额外的数据结构和存储系数值来实现。
如果您想跟踪求解器的执行，这将很有帮助，但结果将取决于求解器积分网格而不是所需的求解时间尺度。
无论如何，如果你想这样做，只需调整你的功能：
def vectorfield2(w, t, p, store):
# ...
store.append({'time': t, 'coefs': [b1, b2]})
# ...
c = []
sol1 = odeint(vectorfield2, w0, t0, args=(p, c), atol=abserr, rtol=relerr)
如您所见，结果不遵循求解时间尺度，而是求解器积分网格：
[{'time': 0.0, 'coefs': [0.4, 45.0]},
{'time': 3.331483023263848e-07, 'coefs': [0.4, 45.0]},
{'time': 3.331483023263848e-07,
'coefs': [0.4000000000026638, 44.999999999955605]},
{'time': 6.662966046527696e-07,
'coefs': [0.4000000000053275, 44.99999999991122]},
{'time': 6.662966046527696e-07,
'coefs': [0.40000000000799124, 44.99999999986683]},
{'time': 0.0001385450759485236,
'coefs': [0.4000002303394594, 44.99999616108455]}, ...]
作为状态函数的系数
另一方面，您的系数仅取决于系统状态，这意味着我们可以在 ODE 函数内部计算它（每次求解器调用它时），然后在时间相关的解决方案已知时计算它。
我建议按如下方式重写您的系统：
def coefs(x):
# Ensure signature:
if isinstance(x, (list, tuple)):
x = np.array(x)
if len(x.shape) == 1:
x = x.reshape(1, -1)
# Compute coefficients
b1 = 0.8 * x[:,0]
b2 = 20. * x[:,2]
q2 = (x[:,2] < 1.0)
b2[q2] = 1.0
return np.stack([b1, b2]).T
def system(t, w, p):
x1, y1, x2, y2 = w
m1, m2, k1, k2, L1, L2, = p
b = coefs(w)
return [
y1, (-b[:,0] * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
y2, (-b[:,1] * y2 - k2 * (x2 - x1 - L2)) / m2,
]
sol3 = solve_ivp(system, [t0[0], t0[-1]], w0, args=(p,), t_eval=t0, method='LSODA')
请注意，coefs 方法的接口旨在对简单状态（列表或向量）和时间解决方案（矩阵）进行操作。
然后，很容易显示时间的解决方案：
fig, axe = plt.subplots()
axe.plot(t0, sol3.y.T)
axe.set_title("Dynamic System: Coupled Spring Mass")
axe.set_xlabel("Time, $t$")
axe.set_ylabel("System Coordinates, $x_i(t)$")
axe.legend([r'$x_%d$' % i for i in range(sol3.y.T.shape[1])])
axe.grid()
和系数：
fig, axes = plt.subplots(2,1, sharex=True, sharey=False)
axes[0].plot(t0, coefs(sol3.y.T)[:,0])
axes[1].plot(t0, coefs(sol3.y.T)[:,1])
axes[0].set_title("Dynamic System: Coupled Spring Mass")
axes[1].set_xlabel("Time, $t$")
axes[0].set_ylabel("$b_0(t)$")
axes[1].set_ylabel("$b_1(t)$")
for i in range(2):
axes[i].grid()
这似乎更有可能是您想要实现的目标。