如何模拟弹跳球? odeint不适用? [英] How to simulate bouncing ball? odeint not applicable?

问题描述

以下是模拟最简单的落球的代码:

Here is a code simulating the simplest falling ball:

%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处弹跳的球怎么样?

But how about a ball bouncing at y=0?

我知道，通常来说，碰撞是物理模拟中的一个困难部分.但是，我想知道是否真的不可能用odeint模拟这个简单的系统.

I know that, in general, collision is a difficult part of physics simulations. However, I wonder if it's really not possible to simulate this simple system, hopefully with odeint.

推荐答案

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

The simplest way is just to add a big force that kicks the particle out from the forbidden region. This requires an ODE solver that is able to handle "stiff" problems, but luckily odeint uses lsoda which automatically switches to a mode suitable for such problems.

例如，

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

如果将步长函数四舍五入，则数值求解器的时间会稍微容易一些，例如，F(x) = big_number / (1 + exp(x/width)):

The numerical solver will have slightly easier time if the step function is rounded a bit, for example, 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()

这篇关于如何模拟弹跳球? odeint不适用?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！