odeint возвращает неправильные результаты для ODE, включая функцию descrete - PullRequest
Новая
       112

odeint возвращает неправильные результаты для ODE, включая функцию descrete

0 голосов
/ 19 февраля 2019

Я пытаюсь смоделировать ODE:

enter image description here

Я реализовал: 

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])

результат:

enter image description here

, что противоречит тому, что я получаю от 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;

enter image description here

Буду признателен, если вы поможете мне узнать, где моя ошибка и как я могу ее исправить.

1 Ответ

0 голосов
/ 21 февраля 2019

Ваша система имеет 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 на вопрос « Как использовать оператор if в дифференциальном уравнении (SciPy)? ».

Чтобы получить определенное поведение, убедитесь, что в каждом событии фазовая граница пересекается при каждом событии (если ускорение не равно нулю (и оно почти всегда не равно нулю)). 

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' и x''=0.05464678 на x'=0.Нет равновесной позиции и нет возможности двигаться вперед во времени.Из-за принудительного прохождения через границу в численном расчете динамика может продолжаться во времени, однако чем меньше размер eps, амплитуда этого колебания, тем меньше длина волны и, следовательно, размер шага.Последнее условие в контуре интегрирования завершает (для гораздо меньшего eps) интегрирование, если длина волны колебаний становится слишком малой.

enter image description here

...