使用RK4进行仿真,在仿真过程中更新ODE变量 [英] Simulation with RK4, update ODE variable as the simulation goes

查看:93
本文介绍了使用RK4进行仿真,在仿真过程中更新ODE变量的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我目前正在制作一个python应用,用于模拟一组依赖于变量的耦合的常微分方程组,我们将其称为"X".

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 ».

至于现在,我基本上是在给定的时间内用 RK4 模拟这组 ODE,然后我在动画图上绘制图形,并在 tkinter 中嵌入了matplotlib 动画".

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.

事实是,我希望能够在求解方程时修改X",以便在我们修改此变量时模拟可以改变.

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.

解决方案从那时起发生了变化.综上所述,这些是氙和碘丰度以及核反应堆中中子流动的方程式.更改的变量是控制控制栏的大sigma_b".

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 :

中子流ODE :

总而言之,«X»属于[1.0,2.0]范围,也就是说我们要进行400小时的模拟:

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

  • 我们启动了 ODE 的模拟
  • 我们每秒更新一次动画图(相当于图形上的1小时仿真)
  • 例如,我们借助滑块来修改«X»(实际上,"X"并没有直接修改,有一个动态过程会逐渐将"X"的初始值移向其选择的值).
  • RK4算法已将其考虑在内
  • 我们可以在更新后的图表上看到变化
  • 等…
  • 模拟完成后,我们将停止绘图更新

希望很清楚.我怎么能这样做?

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

推荐答案

将评论的想法翻译成模型,提供一些代码进行实验

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()

这篇关于使用RK4进行仿真,在仿真过程中更新ODE变量的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆