Simulation with RK4, update ODE variable as the simulation goes

Question

Problem :

I’m currently making a python app that simulates a set of coupled ordinary differentials equations that depends on a variable, let’s call it « X ».

As for now, I’m basically simulating this set of ODE with RK4 for given time then I’m plotting the graph on an animated plot with « matplotlib animation » embedded in tkinter.

The fact is that I would like to be able to modify « X » as the equations are resolved so that the simulation can change as we modify this variable.

Context :

The solution changes from that time on. To put things in context, these are the equations for the xenon and iodine abundance as well as neutrons flow in a nuclear reactor. The variable that change is "big sigma_b" that governs the control bars.

ODE :

xenon and iodine abundance ODE :

neutrons flow ODE :

Summary :

To summarize, « X » belongs to the range [1.0, 2.0], let’s say we want to run 400 hours of simulations :

• We launch the simulation of ODE
• We update the animated plot every second (which is equal to 1 hour of simulation on the graph)
• We modify « X » with the help of a slider for instance (in fact "X" is not directly modified, there is a dynamic process that progressively move the initial "X" value toward its chosen value).
• The RK4 algorithm take that into account
• We can see that change on the updated graph
• etc…
• When the simulation is finished, we stop the plot updating

Hopefully, it’s clear enough. How could I do that ?

Answer 1

Translating the ideas of the comments into a mockup gives some code to experiment with

import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import matplotlib.animation as anim
import numpy as np
from scipy.integrate import odeint

# model is dampened oscillation with control variable u 
# representing the equilibrium position
# x'' + 0.4*x + x == u

def model(x,u): return np.array([x[1], u - 0.4*x[1] - x[0]])

# integrate over about 10 periods
N=250
t = np.linspace(0, 12*np.pi,  N+1)              # time
# arrays to hold the computed values
X = np.zeros([N+1,2])
U = np.zeros([N+1])
# initial values
X[0] = [3, 2]          
U[0] = 2;

# initialize plot
fig, ax = plt.subplots()
plt.subplots_adjust(left=0.25, bottom=0.25)
# initialize graphs
l1, = plt.plot(t[0], X[0,0], '-+b', ms=2 )
l2, = plt.plot(t[0], U[0], '-or', ms=2, lw=1 )
plt.xlim(t[0], t[-1]); plt.ylim(-1,3)

# construct slider
axcolor = 'black'
ax_U = plt.axes([0.25, 0.1, 0.65, 0.03], facecolor=axcolor)
sU = Slider(ax_U, 'U', 1.0, 2.0, valinit=U[0])


def animate(i):
    # read the slider value
    U[i] = sU.val
    # integrate over the sub-interval
    X[i] = odeint(lambda x,t: model(x,U[i]), X[i-1], [t[i-1],t[i]], atol=1e-4, rtol=1e-8)[-1];
    # set the data to plot
    l1.set_data(t[0:i+1],X[0:i+1,0]);
    l2.set_data(t[0:i+1],U[0:i+1]);
    return l1,l2

# start the animation, one should set the parameter to keep it from repeating
anim = anim.FuncAnimation(fig, animate, frames = range(1,N+1), blit=True, interval = 100 )
plt.show()