odeint对于包含降级功能的ODE返回错误结果 [英] odeint returns wrong results for an ODE including descrete function

问题描述

我正在尝试对ODE建模:

I'm trying to model the 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)获得的内容相矛盾:

which contradicts what I get from 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;

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

I would appreciate if you could help me know where is my mistake and how I can solve it.

推荐答案

根据x'的符号，您的系统具有3个平滑子系统或阶段.只要积分停留在这些平稳阶段内，步长控制器就会按预期工作.但是，在相位变化的时刻，步长控制器会看到巨大的变化和其用于调整步长的数量的波动，从而需要减小步长.

Your system has 3 smooth sub-systems or phases according to the sign of x'. As long as the integration stays inside these smooth phases, the step size controller works as expected. At the moment that the phase changes however, the step size controller sees huge changes and oscillations in the quantities that it uses to adapt the step size, requiring to regulate the step size down.

接下来是odeint中的方法，lsode，是隐式的，并且隐式方法的假设是等式的右侧再次足够微分，至少一次.隐式求解器无法到达任何地方，您会在峰值中观察到.

Next comes that the method in odeint, lsode, is implicit, and that the assumption for the implicit method is that the right side of the equation is again sufficiently differentiable, at least once. Failing that the implicit solver can go anywhere, which you observe in the spike.

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

One solution is to eliminate the discontinuities by continuing each phase beyond their boundaries and use the event mechanism of the ODE solver to find the points where the boundary is crossed. To that end introduce a sign/phase selector parameter S and solve the system

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

使用功能e(t)=x'(t)作为事件功能.

using the function e(t)=x'(t) as event function.

# 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)?".

As we need to change the phase selector at the event location, the modus of the event has to be terminal. Then loop through the phase segments and combine them to a global solution like in the answer of chthonicdaemon to the question "How to use if statement in a differential equation (SciPy)?".

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

To get a definitive behavior make sure that at each event the phase boundary is crossed at each event (if the acceleration is non-zero (and it is almost always non-zero)).

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.7880468和x=0.9453532之间.在x'=0上的该点附近，加速度为x''=-0.0453532，对于稍微为正的x'，x''=0.15464678为x''=0.15464678，对于稍微为负的x'和x''=0.05464678，对于x'=0.没有平衡的立场，也没有及时前进的方法.由于强制穿过边界，因此在数值计算中动力学可以及时进行，但是eps的大小越小，该振荡的幅度越大，波长越小，因此步距也越小.如果振荡波长变得太小，积分环路中的最后一个条件将结束积分(对于eps要小得多).

Then plot the solutions, for instance like below. The integration stalls at t=4.7880468 with x=0.9453532. Around this point on x'=0 the acceleration is x''=-0.0453532 for slightly positive x' and x''=0.15464678 for slightly negative x' and x''=0.05464678 on x'=0. There is no equilibrium position and no way to proceed forward in time. Due to the forcing to pass through the boundary, in the numerical computation the dynamic can proceed in time, however the smaller the size of eps, the amplitude of that oscillation, the smaller the wave length and thus the step size. The last condition in the integration loop finishes (for much smaller eps) the integration if the oscillation wave length becomes too small.

这篇关于odeint对于包含降级功能的ODE返回错误结果的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！