Matlab: как получить численное решение системы дифур?

[JustAMan]

Необходимо получить численное решение системы дифференциальных уравнений методом Рунге-Кутты:

$\dot x_1 = \frac{-k_1}{m_1} \cdot x_1 + \frac{k_2}{m_1} \cdot (x_2-x_1)$
$\dot x_2 = \frac{-k_2}{m_2 }\cdot x_1 \cdot (x_2-x_1)$
$\dot x_3 = \dot x_1$
$\dot x_4 = \dot x_2$

Я пытался реализовать код на Матлабе, но не так он работает, как ожидал я. Ожидал синусоиды, получаю почти прямую линию :) Вот мой код:

% Система дифференциальных уравнений:
function dy = oscillator(t, x)
global k1 k2 m1 m2
dy = zeros(4, 1); % 4 equations
dy(1) = x(3);
dy(2) = x(4);
dy(3) = -k1/m1 * x(1) + k2/m1 * (x(2) - x(1));
dy(4) = -k2/m2 * (x(2) - x(1));

% Использование:
>> k1 = 2; k2 = 3; m1 = 2; m2 = 3;
>> [T,Y] = ode45(@oscillator,[0 5],[1 0]);
>> plot(T,Y(:,1),'-o')

Может кто видит ошибку, что не так с этим кодом?

[ShMaxG]

Вектор $dy$ должен выдавать значение правой части диф. уравнений. Почему у Вас

dy(1) = x(3);
dy(2) = x(4);

Это же означает, что $\dot x_1 = x_3$ и $\dot x_2 = x_4$. А еще вектор $x$ должен иметь 4 компоненты, а в коде начальное значение при $t=0$ Вы подставляете вектор $(1,0)^T$ только двух элементов.

-- Ср май 08, 2013 13:00:07 --

И вообще, приведенный код работать не будет, потому что дано неправильное начальное значение иксов, и значения констант не пройдут внутрь функции, для этого в Command Window следует еще в самом начале строчку добавить

global k1 k2 m1 m2

[JustAMan]

ShMaxG в сообщении #721080 писал(а):

Почему у Вас

dy(1) = x(3);
dy(2) = x(4);

Это же означает, что $\dot x_1 = x_3$ и $\dot x_2 = x_4$.

А каким образом сослаться на производные? Что-то вроде этого?

dy(1) = dy(3);
dy(2) = dy(4);

Есть ли в этом вообще какой-то смысл? :)

[ShMaxG]

JustAMan в сообщении #721370 писал(а):

Есть ли в этом вообще какой-то смысл? :)

Будет, если вы сначала определите 3-ю и 4-ю компоненты $dy$, а потом определите так 1-ю и 2-ю.
Но не надо так делать. Пишите все как есть. Задайте первую компоненту $dy$ через $x(1)$ и $x(2)$, затем задайте вторую, тоже через $x(1)$ и $x(2)$, просто по первым двум формулам из Вашего первого сообщения, ничего лишнего. А затем напишите, что

dy(3) = dy(1); dy(4) = dy(2);

ведь ровно это значат третье и четвертое равенство из Вашего поста. Вот, собственно, и все дела.

[JustAMan]

ShMaxG
Спасибо за ответ! Исправляю по-тихоньку. Но что-то 4 компоненты, вместо двух у меня не принимаются:

>> [T,Y] = ode45(@oscillator,[0:0.1:4],[1 0 0 0]);
Error using odearguments (line 93)
OSCILLATOR returns a vector of length 2, but the length of initial
conditions vector is 4. The vector returned by OSCILLATOR and the initial
conditions vector must have the same number of elements.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
 

Это функция должна соответственно принять 4 значения?

[ShMaxG]

Во-первых, попробуйте вторым аргументом в ode45 подавать не [0:0.1:4], а [0 4], ибо как разбивать интервал времени пусть сам решает (да и вообще, туда по-определению этой функции должны подаваться именно границы интервала). Если проблема не исчезнет, напишите сюда еще раз код функции oscillator.

[JustAMan]

Вроде уже побеждаю по-тихоньку, не без Вашей помощи, конечно :)
Код функции сейчас такой:

function dy = oscillator1( t, x )
global k1 k2 m1 m2;
dy = zeros( 4, 1 ); % 4 equations
dy(1) = -k1 / m1 * x(1) + k2 / m1 * ( x(2) - x(1) );
dy(2) = -k2 / m2 * ( x(2) - x(1) );
dy(3) = dy(1);
dy(4) = dy(2);

Я не понимаю откуда вообще параметры t и x беруться в функции oscillator1? Откуда они передаются туда? Я не понимаю где именно они передаются и сколько именно переменных передаётся в эту функцию... :)

В Command Window у меня:

clear all;
global k1 k2 m1 m2;
k1 = 2; k2 = 3; m1 = 2; m2 = 3;
t = 1;
[T,Y] = ode45(@oscillator1,[0 4],[1 1 0 0]);
plot(T,Y(:,1),'-o')

Но синусоиды, по-прежнему, не получается... Получается какой-то спадающий график с y=1 вниз, асимптотически подходя к нулю :) на экспоненту похож развёрнутую слева направо :)

[ShMaxG]

Так, ну коды верные (с точностью до нужного Вам интервала интегрирования и начальных данных и непонятного $t=1$). На первый взгляд синусы не понятно откуда появятся, ибо там отчетливо видно, что первое уравнение имеет вид $\dot{x}_1 = -A x_1 + \ldots$ и $A>0$, поэтому ожидать убывающую экспоненту -- это нормально. Правда, на месте троеточия какая-то фигня, поэтому это будет не совсем экспонента. Я не очень верю, что даже при каких-то других начальных данных будут синусы, но не знаю наверняка.
Но мои экстрасенсорные способности шепчут мне, что Вам были изначально даны 2 диффура второго порядка, которые, для решения на машине, нужно привести сначала к 4 диффурам первого порядка. А такой подход действительно дает синусы. Проверьте, какую задачу Вы решаете?

JustAMan в сообщении #721635 писал(а):

Я не понимаю откуда вообще параметры t и x беруться в функции oscillator1? Откуда они передаются туда? Я не понимаю где именно они передаются и сколько именно переменных передаётся в эту функцию... :)

Они передаются туда из функции ode45. Функция ode45 реализует популярный интегратор диффуров 1-го порядка по вложенным формулам Рунге-Кутты 5(4) Дормана-Принса. Интегратору требуется вычислять значение правой части диффуров в некоторых точках. Вот он смотрит в каких и подставляет нужные значения в функцию oscillator. Вы можете сами это проследить, благо исходный код функции ode45 открыт.

[JustAMan]

Спасибо за объяснение по ode45 и по моим уравнениям! :)

ShMaxG в сообщении #721665 писал(а):

Но мои экстрасенсорные способности шепчут мне, что Вам были изначально даны 2 диффура второго порядка, которые, для решения на машине, нужно привести сначала к 4 диффурам первого порядка. А такой подход действительно дает синусы. Проверьте, какую задачу Вы решаете?

Да, именно так :) Даны два дифференциальных уравнения второго порядка (система уравнений колебания гармонического осциллятора). После чего я привёл их к уравнениям первого порядка, получилось 4 штуки. Как Вы уже заметили, метод Рунге-Кутты принимает уравнения первого порядка, поэтому привёл к первому порядку. Я неправильно записал эту систему в Матлабе?

Там: k1, k2 - жёсткости пружин; m1, m2 - массы колеблющихся элементов; x1, x2 - координаты отклонения от положения равновесия.

t = 1, там не нужно на самом деле, убрал. Пытался отлаживать :)

[ShMaxG]

Вы видимо не правильно записали систему диффуров в своем первом сообщении. Напишите сюда правильную систему. Исходя из нее уже следует смотреть на код в Матлабе.

[JustAMan]

Вроде бы как раз именно такая... :)

$x_1'' = \frac{-k_1}{m_1} \cdot x_1 + \frac{k_2}{m_1} \cdot (x_2-x_1)$
$x_2'' = \frac{-k_2}{m_2 } \cdot (x_2-x_1)$

PS Упс, одна ошибка, всё-таки, была в начальном сообщении - во втором уравнении не было домножения на x1. Но в коде ошибки у меня не было, поэтому пока синусоид ещё нет :)

[ShMaxG]

Ну тогда, как это обычно делается, введите новые переменные $v_1$ и $v_2$ и перепишите свои уравнения по типу

\begin{gather*} \begin{split} &\dot{x}_1 = v_1, \\ &\dot{x}_2 = v_2, \\ &\dot{v}_1 = \ldots, \\ &\dot{v}_2 = \ldots . \end{split} \end{gather*}

Тогда матлабовский икс $x_{\mathrm {matlab}}$ будет вектором $(x_1,x_2,v_1,v_2)^T$

Думаю, Вы уже можете все запрограммировать.

[JustAMan]

ShMaxG в сообщении #721686 писал(а):

Ну тогда, как это обычно делается, введите новые переменные $v_1$ и $v_2$ и перепишите свои уравнения по типу

\begin{gather*} \begin{split} &\dot{x}_1 = v_1, \\ &\dot{x}_2 = v_2, \\ &\dot{v}_1 = \ldots, \\ &\dot{v}_2 = \ldots . \end{split} \end{gather*}

Тогда матлабовский икс $x_{\mathrm {matlab}}$ будет вектором $(x_1,x_2,v_1,v_2)^T$

Думаю, Вы уже можете все запрограммировать.

Что-то всё равно не пойму как это будет выглядеть в Матлабовской функции :) Попробовал переписать на такое:

dy(1) = x(1);
dy(2) = x(2);
dy(3) = -( k1 / m1 ) * x(1) + k2 / m1 * ( x(2) - x(1) );
dy(4) = -( k2 / m2 ) * ( x(2) - x(1) );

но тоже видимо не то... То ли я в переменных совсем запутался... По сути, у меня ведь dy(3) и dy(4) это и есть v1', v2'?

[ShMaxG]

JustAMan в сообщении #721690 писал(а):

По сути, у меня ведь dy(3) и dy(4) это и есть v1', v2'?

Да. Но $x(1)$ и $x(2)$ -- это $x_1$ и $x_2$! Так что поправьте первые две строчки.

[JustAMan]

ShMaxG в сообщении #721695 писал(а):

JustAMan в сообщении #721690 писал(а):

По сути, у меня ведь dy(3) и dy(4) это и есть v1', v2'?

Да. Но $x(1)$ и $x(2)$ -- это $x_1$ и $x_2$! Так что поправьте первые две строчки.

А что получиться в итоге? Опять это же? :)))

dy(1) = dy(3);
dy(2) = dy(4);
dy(3) = -( k1 / m1 ) * x(1) + k2 / m1 * ( x(2) - x(1) );
dy(4) = -( k2 / m2 ) * ( x(2) - x(1) );