Почему численное решение ОДУ отходит от неустойчивого равновесия?

15

Я хочу смоделировать поведение системы, подобной двойному маятнику. Система представляет собой робот-манипулятор с 2 степенями свободы, который не приводится в действие и, следовательно, будет вести себя в основном как двойной маятник, на который действует сила тяжести. Единственное основное отличие двойного маятника состоит в том, что он состоит из двух твердых тел с массой и инерционными свойствами в их центрах масс.

По сути, я запрограммировал ode45под Matlab решить систему ODE следующего типа:

[10000M110M1200100M120M22][Икс˙1Икс˙2Икс˙3Икс˙4]знак равно[Икс2-В1-грамм1Икс4-В2-грамм2]

где - угол первого тела относительно горизонтали, - угловая скорость первого тела; - угол второго тела относительно первого тела, а - угловая скорость второго тела. Все коэффициенты указаны в следующем коде, в и функциях , которые я создал.Икс1Икс2Икс3Икс4rhsfMass

clear all
opts= odeset('Mass',@fMass,'MStateDependence','strong','MassSingular','no','OutputFcn',@odeplot);
sol = ode45(@(t,x) rhs(t,x),[0 5],[pi/2 0 0 0],opts);

function F=rhs(t,x)
    m=[1 1];
    l=0.5;
    a=[0.25 0.25];
    g=9.81;
    c1=cos(x(1));
    s2=sin(x(3));
    c12=cos(x(1)+x(3));
    n1=m(2)*a(2)*l;
    V1=-n1*s2*x(4)^2-2*n1*s2*x(2)*x(4);
    V2=n1*s2*x(2)^2;
    G1=m(1)*a(1)*g*c1+m(2)*g*(l*c1+a(2)*c12);
    G2=m(2)*g*a(2)*c12;

    F(1)=x(2);
    F(2)=-V1-G1;
    F(3)=x(4);
    F(4)=-V2-G2;
    F=F';     
end

function M=fMass(t,x)
    m=[1 1];
    l=0.5;
    Izz=[0.11 0.11];
    a=[0.25 0.25];
    c2=cos(x(3));
    n1=m(2)*a(2)*l;
    M11=m(1)*a(1)^2+Izz(1)+m(2)*(a(2)^2+l^2)+2*n1*c2+Izz(2);
    M12=m(2)*a(2)^2+n1*c2+Izz(2);
    M22=m(2)*a(2)^2+Izz(2);
    M=[1 0 0 0;0 M11 0 M12;0 0 1 0;0 M12 0 M22];
end

Обратите внимание, как я установил начальное условие Икс1 (угол первого тела относительно горизонтали), чтобы система начинала в полностью вертикальном положении. Таким образом, поскольку действует только гравитация, очевидным результатом является то, что система вообще не должна двигаться из этой позиции.

ПРИМЕЧАНИЕ: на всех графиках ниже я построил решения и относительно времени.Икс1Икс3

ode45

Когда я запускаю симуляцию в течение 6 секунд ode45, я получаю ожидаемое решение без каких-либо проблем, система остается на месте и не двигается:

введите описание изображения здесь

Однако, когда я запускаю симуляцию в течение 10 секунд, система начинает двигаться неоправданно:

введите описание изображения здесь

ode23

Затем я запустил симуляцию, ode23чтобы проверить, не исчезла ли проблема. Я получаю такое же поведение, только на этот раз расхождение начинается через 1 секунду:

введите описание изображения здесь

ode15s

Затем я запустил симуляцию, ode15sчтобы увидеть, сохраняется ли проблема, и нет, система кажется стабильной даже в течение 100 секунд:

введите описание изображения здесь

Опять же, ode15sэто только первый порядок и обратите внимание, что есть только несколько шагов интеграции. Так что я запустил еще одну симуляцию в ode15sтечение 10 секунд, но с MaxStepразмером для повышения точности, и, к сожалению, это приводит к тому же результату, что и для обоих, и для .0,01ode45ode23

введите описание изображения здесь

Обычно очевидным результатом этих симуляций будет то, что система останется в своем первоначальном положении, поскольку ничто не мешает ей. Почему происходит это расхождение? Это как-то связано с тем, что системы такого типа хаотичны по своей природе? Это нормальное поведение для odeфункций в Matlab?

jrojasqu
источник
Кроме уравнений, я думаю, что схема также очень поможет понять вопрос.
Никогуаро
Если вы считаете это целесообразным, вы можете принять один из ответов (есть зеленая кнопка).
Ertxiem - восстановить Монику
Вы не говорите, но вы, кажется, заговор x1и x3. (Вставьте сухой комментарий о графиках без легенд и описаний.) Попробуйте построить логарифмы (абсолютные значения) x2и x4.
Эрик Тауэрс

Ответы:

15

Я думаю, что два основных момента уже были высказаны Брайаном и Эртциемом: ваше начальное значение является неустойчивым равновесием, а тот факт, что ваши численные вычисления никогда не бывают точными, обеспечивает небольшое возмущение, которое вызовет нестабильность.

Чтобы немного подробнее описать, как это работает, рассмотрим вашу проблему в форме общей проблемы начальных значений.

Y˙(T)знак равноM-1е(T,Y(T))
где иY(T)знак равно(Икс1(T),Икс2(T),Икс3(T),Икс4(T))

f(t,y(t))=[x2V1G1x4V2G2]

В точной арифметике вы должны иметь и, следовательно, и ничто не изменится: ваша система остается в равновесии , Тем не менее, арифметика в компьютере не является точной, а это означает, что по разным причинам ваша правая часть не точно равна нулю, но равна некоторому который почти, но не совсем равен нулю.f(0,y0)=0y˙(0)=0f~

В своем коде вы можете проверить это, вычисляя

norm(rhs(0,[pi/2 0 0 0]))

что дает 6.191e-16 - так почти, но не точно ноль. Как это влияет на динамику вашей системы?

Согласно некоторым предположениям, эффект того, что не является точно нулем, является тем же, если вы не начнете с начального значения вы прописываете, но со значения, которое немного отличается, давайте вызовем это .fy0y~0

Кроме того, за очень короткое время решение вашей системы выглядит как решение линеаризованной системы

Y˙(T)знак равное(0,Y0)+е'(0,Y0)(Y(T)-Y0)знак равное'(0,Y0)(Y(T)-Y0)

где - это якобиан вашей функции или вашего кода. Поскольку является константой, мы можем преобразовать это в уравнение для , где говорит, насколько далеко от первоначального значения мы находимся:е'еrhsY0d(T)знак равноY(T)-Y0d

d˙(T)знак равное'(0,Y0)d(T),

Я не мог потрудиться вычислить якобиан вручную, поэтому я использовал автоматическое дифференцирование, чтобы получить хорошее приближение:

Jзнак равное'(0,Y0)знак равно[01009,8102,4525000012,452502,45250]

так что ваше уравнение становится

d˙(T)знак равноJd(T),d(0)знак равноY~0-Y0

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

Jзнак равноQDQ-1

где - диагональная матрица, содержащая собственные значения и - ортогональные матрицы, представляющие преобразования координат. Затем мы можем преобразовать уравнение для в уравнение для которое читаетDJQdе(T)знак равноQ-1d(T)

е˙(T)знак равноDе(T),е(0)знак равноQ-1d0,

Поскольку является диагональю, это фактически четыре независимых уравненияD

е˙я(T)знак равноλяея(T),ея(0)знак равноя-й компонент Q-1d0

с . Если вы вычислите собственные значения, вы обнаружите, что наибольшее значение равно . Следовательно,язнак равно1,2,3,4λ1знак равно3,2485

е1(T)знак равное1(0)е3,2485T,

Теперь, если бы арифметика на вашем компьютере была точной, вы бы имели , поэтому и, следовательно, и ничего не произойдет. Но поскольку это не так, у вас есть небольшое, но конечное которое экспоненциально усиливается. Отсюда и быстрое отклонение от равновесия в вашем решении.d(0)знак равно0е(0)знак равноQ-1d(0)знак равно0е1(0)знак равно0е1(0)

Даниил
источник
16

Обратите внимание, что представлен в формате двойной точности таким образом, что он точно не равен . Это только с точностью до 15 цифр. Таким образом, вы начинаете чуть-чуть от положения равновесия. Поскольку равновесие неустойчиво, оно в конечном итоге начнет двигаться.π/2π/2

Брайан Борхерс
источник
4
Если вы внимательно следите за переменными состояния (глядя на значения, напечатанные в научной записи), вы сможете увидеть начальное очень медленное движение от равновесия.
Брайан
Это имеет смысл, и действительно, когда я запускаю систему в вертикальном положении вниз (являясь точкой стабильного равновесия), система вообще не движется, по крайней мере, для имитации 1000 секунд, которую я считаю очень длительным периодом времени. ,
jrojasqu
6
Икс1sin(0)cos(0)sin(pi/2)cos(pi/2)rhs(t,[0,0,0 0] == [0,0,0,0]
π/2
1
θзнак равно0 10-16
4

Посмотрите на составляющие сил, рассчитанных в ваших функциях.

π

10-16

aзнак равно1,0aзнак равноa+10-16

alephzero
источник
4

Первоначальное предположение состояло в том, что исходное положение находилось в устойчивом равновесии (т. Е. Минимуме потенциальной энергии) с нулевой кинетической энергией, и система начала отходить от равновесия.
Поскольку физически это не может произойти (если мы рассмотрим классическую механику), мне пришло в голову две вещи:

  1. π/2-π/2

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

Ertxiem - восстановить Монику
источник
1
π
2
Кстати, просто для удовольствия, если вы хотите сохранить систему в нестабильном вертикальном положении, вы можете изменить начало координат, чтобы угол, равный нулю, был направлен вверх.
Ertxiem - восстановить Монику
@Ertxiem другой вариант заключается в том, чтобы ввести небольшое трение в штырьки, что приведет к числовым ошибкам.
Свавил
грех(π)
Поскольку физически это не может произойти - учитывая понимание того, что мы находимся в нестабильном равновесии, я несколько оспариваю это. Физические системы (без слишком большого трения) не остаются в неустойчивых равновесиях. В целом, если вы моделируете реальные системы, вам следует избегать того, чтобы они случайно застряли в неустойчивом равновесии (как бы там ни было) - это особенность, а не ошибка. (Есть некоторые редкие исключения из этого, такие как неинфицированное состояние в иммунологии, которое является нестабильным равновесием, которое можно поддерживать.)
Wrzlprmft
0

Вам следует поискать больше о двойных маятниках: это то, что мы называем «хаотическими системами». Несмотря на то, что они ведут себя по простым правилам, начиная с немного разных начальных условий, решения расходятся довольно быстро. Выполнение численного моделирования для систем такого типа не легко. Посмотрите следующее видео, чтобы получить более полное представление о проблеме.

Для простого или двойного маятника вы можете написать формулу для полной энергии системы. Предполагая, что силами трения пренебрегают, эта полная энергия сохраняется аналитической системой. Численно это совсем другая проблема.

Прежде чем попробовать двойной маятник, попробуйте простой маятник. Вы заметите, что для методов Рунге-Кутты малого порядка энергия системы будет расти в численном моделировании вместо того, чтобы оставаться постоянной (это то, что происходит в ваших симуляциях: вы получаете движение из ничего). Чтобы предотвратить это, можно использовать методы RK более высокого порядка (ode45 имеет порядок 4; RK порядка 8 будет работать лучше). Существуют и другие методы, называемые «симплектическими методами», которые разработаны таким образом, чтобы численное моделирование сохраняло энергию. В общем, вы должны остановить симуляцию, как только энергия значительно возрастет по сравнению с вашей инициализацией.

И попытайтесь понять простой маятник, прежде чем перейти к двойному.

Бени Богосел
источник
2
Это не вопрос хаотичности системы. Вы можете иметь неустойчивое равновесие и в не хаотических системах, например, один маятник находится «на голове», и он будет демонстрировать то же поведение, что и в вопросе.
Даниэль
1
Также неверно, что энергия увеличивается для РКМ малого порядка: неявный Эйлер имеет первый порядок и демонстрирует совершенно противоположное поведение.
Даниэль
@BeniBogosel Вы упоминаете симплектические методы, которые привлекли мое внимание, потому что, очевидно, в моем примере энергия не сохраняется. Однако не могли бы вы указать конкретный симплектический метод, который мог бы быть реализован здесь?
jrojasqu
@jrojasqu, почему вы говорите, что энергия не сохраняется в вашей системе?
Ertxiem - восстановить Монику
ode45π