如何模拟弹跳球？ odeint 不适用？答案

【问题标题】：How to simulate bouncing ball? odeint not applicable?如何模拟弹跳球？ odeint 不适用？
【发布时间】：2015-02-20 22:32:12
【问题描述】：

这是一个模拟最简单的落球的代码：

%pylab
from scipy.integrate import odeint
ts = linspace(0, 1)
def f(X, t):
    dx0 = X[1]
    dx1 = -9.8
    return [dx0, dx1]
X = odeint(f, [2, 0], ts)
plot(ts, X[:, 0])

但是球在 y=0 处弹跳怎么样？

我知道，一般来说，碰撞是物理模拟的一个困难部分。但是，我想知道是否真的无法模拟这个简单的系统，希望使用 odeint。

【问题讨论】：

您的f(X,t) 独立于t，这很有意义。检查docs.scipy.org/doc/scipy/reference/tutorial/… 以获取类似示例

标签： python scipy ode

【解决方案1】：

最简单的方法就是添加一个很大的力，将粒子踢出禁区。这需要一个能够处理“僵硬”问题的 ODE 求解器，但幸运的是 odeint 使用 lsoda 自动切换到适合此类问题的模式。

例如，

F(x) = { 0              if x > 0
       { big_number     if x < 0

如果阶跃函数稍微四舍五入，数值求解器会稍微容易一些，例如F(x) = big_number / (1 + exp(x/width))：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

big_number = 10000.0
width = 0.0001

ts = np.linspace(0, 10, 2000)
def f(X, t):
    dx0 = X[1]
    dx1 = -9.8
    dx1 += big_number / (1 + np.exp(X[0]/width))
    return [dx0, dx1]

with np.errstate(over='ignore'):
    # don't print overflow warning messages from exp(),
    # and limit the step size so that the solver doesn't step too deep
    # into the forbidden region
    X = odeint(f, [2, 0], ts, hmax=0.01)
plt.plot(ts, X[:, 0])
plt.show()

【讨论】：