我用scipy.integrate.odeint来模拟双摆翻转时间.按照我的代码 struct ,odeint在固定的时间间隔内求解系统;然后我必须判断结果,看看是否发生翻转.相反,我想在每一步之后判断条件是否满足,如果满足,停止.然后,没有必要解决剩余的时间间隔.

以下是我当前的代码:

from matplotlib.colors import LogNorm, Normalize
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from tqdm import tqdm
import seaborn as sns
import numpy as np

def ode(y, t, length_1, length_2, mass_1, mass_2, gravity):
    angle_1, angle_1_d, angle_2, angle_2_d = y
    angle_1_dd = (-gravity * (2 * mass_1 + mass_2) * np.sin(angle_1) - mass_2 * gravity * np.sin(angle_1 - 2 * angle_2) - 2 * np.sin(angle_1 - angle_2) * mass_2 * (angle_2_d ** 2 * length_2 + angle_1_d ** 2 * length_1 * np.cos(angle_1 - angle_2))) / (length_1 * (2 * mass_1 + mass_2 - mass_2 * np.cos(2 * angle_1 - 2 * angle_2)))
    angle_2_dd = (2 * np.sin(angle_1 - angle_2) * (angle_1_d ** 2 * length_1 * (mass_1 + mass_2) + gravity * (mass_1 + mass_2) * np.cos(angle_1) + angle_2_d ** 2 * length_2 * mass_2 * np.cos(angle_1 - angle_2))) / (length_2 * (2 * mass_1 + mass_2 - mass_2 * np.cos(2 * angle_1 - 2 * angle_2)))

    return [angle_1_d, angle_1_dd, angle_2_d, angle_2_dd]

def double_pendulum(length_1, length_2, mass_1, mass_2, angle_1_init, angle_2_init, angle_1_d_init, angle_2_d_init, gravity, dt, num_steps):
    time_span = np.linspace(0, dt*num_steps, num_steps)
    y0 = [np.deg2rad(angle_1_init), np.deg2rad(angle_1_d_init), np.deg2rad(angle_2_init), np.deg2rad(angle_2_d_init)]
    sol = odeint(ode, y0, time_span, args=(length_1, length_2, mass_1, mass_2, gravity))

    return sol 

def flip(length_1, length_2, mass_1, mass_2, angle_1_init, angle_2_init, angle_1_d_init, angle_2_d_init, gravity, dt, num_steps):
    solution = double_pendulum(length_1, length_2, mass_1, mass_2, angle_1_init, angle_2_init, angle_1_d_init, angle_2_d_init, gravity, dt, num_steps)
    #angle_1 = solution[:, 0]
    angle_2 = solution[:, 2]
    for index, angle2 in enumerate(angle_2):
        if abs(angle2 - angle_2_init) >= 2*np.pi:
            return index*dt
    return dt*num_steps

angle_1_range = np.arange(-172, 172, 1)
angle_2_range = np.arange(-172, 172, 1)

fliptime_matrix = np.zeros((len(angle_1_range), len(angle_2_range)))

for i, angle_1 in tqdm(enumerate(angle_1_range), desc='angle_1'):
    for j, angle_2 in tqdm(enumerate(angle_2_range), desc='angle_2', leave=False):
        fliptime = flip(1, 1, 1, 1, angle_1, angle_2, 0, 0, 9.81, 0.01, 10000)
        fliptime_matrix[i, j] = fliptime

sns.heatmap(fliptime_matrix, square=True, cbar_kws={'label': 'Divergence'}, norm=LogNorm())
plt.xlabel('Angle 2 (degrees)')
plt.ylabel('Angle 1 (degrees)')
plt.title('Fliptime Heatmap')
plt.gca().invert_yaxis()
plt.show()

如何在满足翻转条件(angle2中的变化大于2*pi)后立即停止积分?

推荐答案

在初始值问题的语言中,您希望检测一个"事件".scipy.integrate.odeint并不提供接口.这是文件建议的原因之一:

对于新代码,使用scipy.integrate.solve_ivp解微分方程.

一旦您将代码转换为使用solve_ivp,您可以使用events参数在检测到事件后终止集成.事件函数看起来像:

def event(t, y):
    angle2 = y[2]
    return abs(angle2 - angle_2_init) - 2*np.pi
event.terminal = True

Python相关问答推荐

线性模型PanelOLS和statmodels OLS之间的区别

Pandas 在最近的日期合并,考虑到破产

比较两个数据帧并并排附加结果(获取性能警告)

如何使用matplotlib在Python中使用规范化数据和原始t测试值创建组合热图?

Pandas 有条件轮班操作

通过Selenium从页面获取所有H2元素

在Pandas DataFrame操作中用链接替换'方法的更有效方法

将输入管道传输到正在运行的Python脚本中

发生异常:TclMessage命令名称无效.!listbox"

组/群集按字符串中的子字符串或子字符串中的字符串轮询数据框

转换为浮点,pandas字符串列,混合千和十进制分隔符

未调用自定义JSON编码器

python—telegraph—bot send_voice发送空文件

用SymPy在Python中求解指数函数

提高算法效率的策略?

从一个df列提取单词,分配给另一个列

裁剪数字.nd数组引发-ValueError:无法将空图像写入JPEG

Seaborn散点图使用多个不同的标记而不是点

极柱内丢失类型信息""

大Pandas 中的群体交叉融合