如何模拟弹跳球? odeint 不适用?
How to simulate bouncing ball? odeint not applicable?
下面是模拟最简单的落球的代码:
%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。
最简单的方法就是加一个大的力把粒子踢出禁区。这需要能够处理 "stiff" 问题的 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()
下面是模拟最简单的落球的代码:
%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。
最简单的方法就是加一个大的力把粒子踢出禁区。这需要能够处理 "stiff" 问题的 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()