2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу 1, 2  След.
 
 Matlab: как получить численное решение системы дифур?
Сообщение07.05.2013, 23:38 
Необходимо получить численное решение системы дифференциальных уравнений методом Рунге-Кутты:

$\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')

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение08.05.2013, 11:54 
Аватара пользователя
Вектор $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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение08.05.2013, 23:51 
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);

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 11:47 
Аватара пользователя
JustAMan в сообщении #721370 писал(а):
Есть ли в этом вообще какой-то смысл? :)

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 17:23 
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 значения?

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 19:10 
Вроде уже побеждаю по-тихоньку, не без Вашей помощи, конечно :)
Код функции сейчас такой:
Код:
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 вниз, асимптотически подходя к нулю :) на экспоненту похож развёрнутую слева направо :)

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 20:48 
Аватара пользователя
Так, ну коды верные (с точностью до нужного Вам интервала интегрирования и начальных данных и непонятного $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 открыт.

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 21:12 
Спасибо за объяснение по ode45 и по моим уравнениям! :)

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

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

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

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 21:22 
Аватара пользователя
Вы видимо не правильно записали систему диффуров в своем первом сообщении. Напишите сюда правильную систему. Исходя из нее уже следует смотреть на код в Матлабе.

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 21:39 
Вроде бы как раз именно такая... :)

$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. Но в коде ошибки у меня не было, поэтому пока синусоид ещё нет :)

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 21:51 
Аватара пользователя
Ну тогда, как это обычно делается, введите новые переменные $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$

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 21:59 
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'?

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 22:05 
Аватара пользователя
JustAMan в сообщении #721690 писал(а):
По сути, у меня ведь dy(3) и dy(4) это и есть v1', v2'?

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

 
 
 
 Re: Matlab: как получить численное решение системы дифур?
Сообщение09.05.2013, 22:13 
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) );

 
 
 [ Сообщений: 21 ]  На страницу 1, 2  След.


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group