python - odeint 为包含离散函数的 ODE 返回错误结果

Question

我正在尝试为 ODE 建模：

我实现了：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
m = 1
k = 1
M = 0.1
b = 1
Fmax = 1
def dXdt(X,t):
    return [X[1], - b * X[1] / m - k * X[0] / m - M * np.sign(X[1]) / m + Fmax / m ]

X0 = [1, 2]
ts = np.linspace(0, 10, 200)
Xs = odeint(dXdt, X0, ts)
plt.plot(ts, Xs[:, 0])

结果：

这与我从 Modelica (OpenModelica) 得到的内容相矛盾：

model model1
//constants
  parameter Real m = 1;
  parameter Real k = 1;
  parameter Real b = 1;
  parameter Real M = 0.1;
  parameter Real Fmax = 1;
//variables
  Real x, v, a;
initial equation
  x = 1;
  v = 2;
equation
  v = der(x);
  a = der(v);
  m * a + b * v + k * x + M * sign(v) = Fmax;
end model1;

如果您能帮助我知道我的错误在哪里以及如何解决它，我将不胜感激。

标签： pythonscipyodeodeint

Answer 1

根据 的符号，您的系统有 3 个平滑的子系统或阶段x'。只要积分保持在这些平滑阶段内，步长控制器就会按预期工作。然而，在相位发生变化的那一刻，步长控制器会看到用于调整步长的量的巨大变化和振荡，因此需要将步长调低。

接下来是 , , 中的方法odeint是lsode隐式的，并且隐式方法的假设是等式的右侧再次充分可微，至少一次。如果隐式求解器不能去任何地方，你会在尖峰中观察到。

---

一种解决方案是通过继续超出边界的每个阶段来消除不连续性，并使用 ODE 求解器的事件机制来找到跨越边界的点。为此引入符号/相位选择器参数S并求解系统

m*x''+b*x'+k*x+M*s = F

将该函数e(t)=x'(t)用作事件函数。

# Define differential equation
m,b,k,M,F = 1., 1., 1., 0.1, 1.
def fun(t, x, S):
    dx = [x[1], (F-b*x[1]-k*x[0]-M*S)/m]
    return np.asarray(dx)

# Define event function and make it a terminal event
def event(t, x):
    return x[1]
event.terminal = True

t = 0
x = [1.,2.]; 
S = np.sign(u[1]);
tend = 10

由于我们需要在事件位置更改相位选择器，因此事件的方式必须是terminal. 然后循环遍历相位段并将它们组合成一个全局解决方案，就像 chthonicdaemon 对问题“如何在微分方程 (SciPy) 中使用 if 语句？ ”的回答一样。

要获得确定的行为，请确保在每个事件中都跨越相边界（如果加速度不为零（并且几乎总是非零））。

ts = []
xs = []
eps=1e-8
for _ in range(50):
    sol = solve_ivp(lambda t,u:fun(t,u,S), (t, tend), x, events=event, atol=1e-12, rtol=1e-8, max_step=0.01);
    ts.append(sol.t)
    xs.append(sol.y)
    if sol.status == 1: # Event was hit
        # New start time for integration
        t = sol.t[-1]
        # Reset initial state
        x = sol.y[:, -1].copy()
        # ensure the crossing of the phase boundary
        dx = fun(t,x,S)
        dt = -(eps*S+x[1])/dx[1]; # should be minimal
        if dt > 0: t += dt; x += dt*dx;
        # new phase parameter
        S = np.sign(x[1]);
        # stop the iteration if it stalls 
        if t-sol.t[0] <5e-12: break
    else:
        break

# We have to stitch together the separate simulation results for plotting
t=np.concatenate(ts);
x=np.concatenate(xs, axis=1);

然后绘制解决方案，例如如下所示。集成在t=4.7880468with处停止x=0.9453532。围绕这一点x'=0的加速度是x''=-0.0453532略为正x'，x''=0.15464678略为负x'和x''=0.05464678上x'=0。没有平衡的位置，也没有办法及时前行。由于强制通过边界，在数值计算中动态可以及时进行，但是 的大小越小eps，振荡的幅度越小，波长越小，因此步长越小。eps如果振荡波长变得太小，积分循环中的最后一个条件完成（对于更小）积分。