交互式matplotlib小部件不更新scipy代码求解器

Question

我一直在使用交互式matplotlib小部件来可视化微分方程的解。我已经得到了它与odeint函数在scipy工作，但没有设法让它与ode类更新。我宁愿使用后者，因为它可以更好地控制使用哪个求解器。

下面的代码用于求解一个指数衰减的微分。y0是衰减的振幅。当在更新函数中调用solve . integration (t1)时，代码停止工作。我不知道为什么会这样。

from scipy.integrate import ode
# solve the system dy/dt = f(t, y)
def f(t, y):
    return -y / 10
# Array to save results to
def solout(t, y):
    sol.append([t, *y])
solver = ode(f).set_integrator('dopri5')
solver.set_solout(solout)
# Initial conditions
y0 = [1]  # Initial amplitude
t0 = 0      # Start time
t1 = 20    # End time
fig = plt.figure(figsize=(10, 6))
fig.subplots_adjust(left=0.25, bottom=0.4)
ax = plt.subplot(111)
# solve the DEs
solver.set_initial_value(y0, t0)
sol = []
solver.integrate(t1)
sol = np.array(sol)
t = sol[:, 0]
y = sol[:, 1]
l, = plt.plot(t, y, lw=2, color='red')
plt.axis([0, 20, 0, 1.1])
plt.xlabel('Time (ms)')
plt.ylabel('n1(t)')
plt.grid()
axcolor = 'lightgoldenrodyellow'
axn1 = plt.axes([0.25, 0.1, 0.65, 0.03], axisbg=axcolor)
sn1 = Slider(axn1, 'y(0)', 0, 1.0, valinit=1)
def update(val):
    y0 = [sn1.val]
    solver.set_initial_value(y0, t0)
    sol = []
    solver.integrate(t1)
    sol = np.array(sol)
    t = sol[:, 0]
    y = sol[:, 1]
    l.set_data(t, y)
    plt.draw()
sn1.on_changed(update)

Answer 1

我想把计算和绘图分开总是明智的。因此，首先尝试在给定初始条件下求解ODE。一旦成功，尝试绘制和交互的东西。

在您的例子中，我们将构建一个函数来解决ODE，然后在绘图更新中使用具有不同初始条件的该函数。

import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import numpy as np
from scipy.integrate import ode

# solve the system dy/dt = f(t, y)
def f(t, y):
    a = np.zeros((1,1))
    a[0] = -y / 10.
    return a
#define a function to solve the ODE with initial conditions
def solve(t0, t1, y0, steps=210):
    solver.set_initial_value([y0], t0)
    dt = (t1 - t0) / (steps - 1)
    solver.set_initial_value([y0], t0)
    t = np.zeros((steps, 1))
    Y = np.zeros((steps, 1))
    t[0] = t0
    Y[0] = y0
    k = 1
    while solver.successful() and k < steps:
        solver.integrate(solver.t + dt)
        t[k] = solver.t
        Y[k] = solver.y[0]
        k += 1
    return t, Y
# set the ODE integrator   
solver = ode(f).set_integrator("dopri5")  
# Initial conditions
y0 = 1.  # Initial amplitude
t0 = 0.      # Start time
t1 = 20.    # End time
#solve once for given initial amplitude
t, Y = solve(t0, t1, y0)

fig = plt.figure(figsize=(10, 6))
fig.subplots_adjust(left=0.25, bottom=0.4)
ax = plt.subplot(111)
l, = plt.plot(t, Y, lw=2, color='red')
plt.axis([0, 20, 0, 1.1])
plt.xlabel('Time (ms)')
plt.ylabel('n1(t)')
plt.grid()
axn1 = plt.axes([0.25, 0.1, 0.65, 0.03], axisbg='#e4e4e4')
sn1 = Slider(axn1, 'y(0)', 0, 1.0, valinit=1)
def update(val):
    #solve again each time 
    t, Y = solve(t0, t1, sn1.val)
    l.set_data(t, Y)
    plt.draw()
sn1.on_changed(update)
plt.show()