交互式matplotlib小部件沒有用scipy ode求解器更新

Question

我一直在使用交互式matplotlib小部件來可視化微分方程的解。我已經在scipy中使用了odeint函數，但是沒有設法使用ode類來更新它。我寧願使用後者，因爲它可以更好地控制使用哪個求解器。

以下代碼用於求解指數衰減的微分。 y0是衰減的幅度。當更新函數中調用solver.integrate（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()