Решение дифференциальных уравнений в scilab

Решение дифференциальных уравнений в Scilab

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

Мы рассмотрим способы исследования сложных механических систем, моделируемых дифференциальными уравнениями. Рассмотрим способы конструирования управления нескольких видов для механической системы с двумя степенями свободы. А также, проведем визуализацию процесса стабилизации движения.

Дифференциальные уравнения n-го порядка

Дифференциальным уравнением n-го порядка называется соотношение вида

Решением дифференциального уравнения является функция \( x(t) \), которая обращает это уравнение в тождество. Дифференциальные уравнения имеют бесконечное множество решений.

Чтобы однозначно определить решение некоторой задачи, необходимо задавать начальные условия. Количество таких условий должно совпадать с порядком уравнения. А решение, найденное для таких условий называется решением задачи Коши.

6.2 Пример решения оду в Scilab

Для решения обыкновенных дифференциальных уравнений (ОДУ) в Scilab используется
функция y = ode([type], y0, t0, t, func).
Разберём, что обозначает каждый из параметров у этой функции.

type — необязательный строковый параметр, с помощью которого можно выбрать метод решения ОДУ. Обычно этот параметр опускается.
t0 — скаляр начальный момент отрезка интегрирования. Обычно \( t0 = 0 \).
y0 — начальные условия. Отметим, что \(y0 \) это вектор, размерность которого совпадает с порядком ОДУ. Так, для ОДУ 2-го порядка необходимо задать значения в начальный момент времени для функции и её производной, т.е. использовать запись \(y0 = [0.1, 0.3]\).
t — вектор, задающий узлы сетки, в которых ищем решение. Как правило, вектор t задается следующим образом t=t0:d:tmax, где \(t0\) — начальный момент отрезка интегрирования, \(d\) — шаг дискретизации, \(tmax\) — конечный момент отрезка интегрирования.
func — пользовательская функция, определяющая правую часть уравнения.
y — вектор решений.

Рассмотрим использование функции ode() на примере решения следующей задачи.

Найти угловую скорость \( \omega(t) \) твердого тела вокруг неподвижной оси, если заданы начальная угловая скорость тела \( \omega_0=1 \) и угловое ускорение \( \varepsilon(t) = 2+0.5sin(t) \). Найти значение \( \omega(t) \) в момент времени \( t_1=1.3 \). Построить график функций.

Рисунок 9. Сравнение графиков угловой скорости точки, по функции, найденной аналитическим и численным интегрированием. Красным кружком обозначено значение скорости точки в момент времени t=1.3c.

Угловая скорость точки может быть найдена из дифференциального уравнения

В нашем случае по условию задачи указаны следующие параметры для функции ode():
\(t0=0 \) начальный момент времени,
\(y0 = 1\) начальное условие одно, т.к. порядок уравнения равен 1, это начальная угловая скорость тела \( \omega_0=1 \)
t=t0:h:tmax — вектор, задающий узлы сетки.

Результат работы программы представлен на рис.9, исходный код на листинге 15.

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

Функция \( f = 2t — 0,5cos(t) + 1,5 \) единственное решение дифференциального уравнения, удовлетворяющее заданным начальным условиям, то есть, функция \(f\) решение задачи Коши.

6.3 Cистемы дифференциальных уравнений

Для решения систем ОДУ в Scilab используется та же функция y = ode([type], y0, t0, t, func). Однако важным требованием является запись исследуемой системы в нормальной форме Коши.

Системой дифференциальных уравнений, записанной в нормальной форме Коши называется система, где слева стоят производные фазовых переменных, а справа некоторые функции, т.е. система вида

\begin x_1′ = f_1(t,x_1. x_n) \\ x_2′ = f_2(t,x_1. x_n) \\ . \\ x_n’ = f_n(t,x_1. x_n) \end

Решением системы ОДУ является вектор x(t), который обращает это уравнение в тождество. Размерность вектора равна количеству уравнений в системе.

Стоит отметить, что дифференциальное уравнение n-ой степени может быть представлено в виде системы из n-уравнений первой степени, что позволяет решать задачу Коши для полученной системы ОДУ.

6.4 Примеры поиска решения систем ОДУ в Scilab

Решение системы ОДУ

Рассмотрим решение задачи Коши для системы уравнений

\begin x’ = cos(xy) \\ y’ = sin(x+ty) \end

на интервале \( [0; 10]\) и с начальными условиями \(x(0)=0, y(0)=0 \).

Для поиска решения данной задачи, нам необходимо привести исходную систему к нормальному виду Коши. Для этого введём новые переменные \( (z1, z2) \) и сделаем необходимые переобозначения в исходной системе:

Рассмотрим решение задачи Коши для системы уравнений

\begin z_1 = x \\ z_2 = y \end \begin z_1′ = cos(z_1 z_2) \\ z_2′ = sin(z_1 + t z_2) \end

Код программы, реализующей поиск решения системы дифференциальных уравнений представлен на листинге 16. Обратите внимание, как происходит обращение к компонентам вектора решения \( z(t) \).

Для наглядной реализации сформировано графическое решение на рис. 10.

Рисунок 10. Графическое решение задачи Коши с помощью функции ode().

Решение системы ОДУ в матричной форме

Рассмотрим решение задачи Коши для системы уравнений, заданных в матричном виде

на интервале \( [0; 10] \) и начальными условиями \( y_1(0)=1, y_2(0)=0 \).

Код программы, реализующей поиск решения системы дифференциальных уравнений представлен на листинге 17. Обратите внимание, как происходит обращение к компонентам вектора решения \(y(t)\).

Для наглядной реализации сформировано графическое решение на рис. 11.

Рисунок 11. Графическое решение задачи Коши с помощью функции ode().

Решение ДУ 2-го порядка путём сведения к системе уравнений

Продолжим знакомиться с возможностями функции ode() на примере решения задачи механики на второй закон Ньютона.

Груз находится на пружине жёсткости \( c=12 \) H/м, масса груза \( m=68.7 \) кг. Определить закон движения груза, если на него действует сила \( F=100.5sin(2t) + <2\pi \over 3>\) H.

По 2-му закону Ньютона, движение груза описывается с помощью дифференциального уравнения 2-й степени, которое имеет вид:

и начальными условиями \( x(0)= 0, ;\ \dot(0)= 0 \).

Как известно, ОДУ второго порядка сводится к системе в нормальной форме Коши, состоящей из двух уравнений первой степени. Введём новые переменные \( (z_1, z_2) \) и сделаем необходимые переобозначения в исходной системе:

\begin z_1 = x \\ z_2 = \dot \end \begin \dot_1 = z_2 \\ \dot_2 = <1 \over m>(-cz_1 + F(t)) \end

Код программы, реализующей поиск решения системы дифференциальных уравнений представлен на листинге 18. Обратите внимание, как происходит обращение к компонентам вектора решения \( y(t) \).

Для наглядной реализации сформировано графическое решение на рис. 12. Компонента \( z_1 \) представляет собой координату груза, а компонента \( z_2 \) — скорость груза.

Рисунок 12. Графики движения груза на пружине для компонент \( z_1 и z_2 \) соответственно.

Решение дифференциальных уравнений в scilab

Для решения ОДУ и системы ОДУ 1-го порядка в Scilab предусмотрена функция ode , имеющая форматы:

[y]=ode(y0,t0,t,f),

[y,w,iw]=ode(type,y0,t0,f,adams,stiff,rk,rkf,fix,rtol,adol,jac,w,iw),

которые содержат обязательные и необязательные параметры.

Первый формат функцииode содержит только обязательные параметры, к которым относятся:

y 0 – вектор начальных условий;

t 0 – начальная точка интервала интегрирования;

t – координаты узлов сетки, в которых происходит поиск решения;

f– имя внешней функция, определяющей правую часть уравнения

или системы уравнений;

y – вектор решений (выходной параметр).

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

type– строка, указывающая тип используемой программы решения, может принимать значения: » adams «, » stiff «, » rk «, » rkf «, » fix «, » discrete «или » root «;

adams– применяют при решении дифференциальных уравнений или систем методом прогноза и коррекции Адамса;

stiff– указывают при решении жестких задач;

rk– используют при решении дифференциальных уравнений или систем методом Рунге-Кутты четвертого порядка;

rkf– указывают при выборе пятиэтапного метода Рунге-Кутты четвертого порядка;

fix– тот же метод Рунге-Кутты, но с фиксированным шагом;

rtol, atoll– относительная и абсолютная погрешности вычислений, соответственно, по умолчанию rtol =0.00001, atol =0.0000001 (при использовании параметров rkfи fix– rtol =0.001, atol =0.0001);

jac– матрица, представляющая собой якобиан правой части жесткой системы дифференциальных уравнений и заданная в виде внешней функции вида j = jak ( t , y );

w, iw– векторы, предназначенные для сохранения информации о параметрах интегрирования, которые применяют для того, чтобы последующие вычисления выполнялись с теми же параметрами.

Описание этого формата достаточно подробно рассмотрено в справочной системе Scilab, где и приведены примеры его использования для различных типов ОДУ [13].

Остановимся на использовании первого формата функции ode, для чего рассмотрим решение следующего примера: найти решение ОДУ y ‘=- sin ( x * y ) на отрезке [0;5]cшагом 0.5 при начальных условиях y (0)=1.5.

Решение ОДУ, полученное с применением функции ode и график полученного решения, представлены на рис.2.5.2-1.

—> // Решение ОДУy’=-sin(x*y) —> y0 = 1.5; t0 = 0; t = 0:0.5:5; // Начальные условия —> // Загрузка и выполнение сценарияРИС2521 —> exec(‘РИС2521.sce’, 0); ans = 0. 1.5 0.5 1.3302706 1. 0.9566472 1.5 0.5574285 2. 0.2477507 2.5 0.0822207 3. 0.0208664 3.5 0.0041103 4. 0.0006303 4.5 0.0000753 5. 0.000007

Рис. 2.5.2-1. Решение ОДУ с использованием функции ode

Для решения системы обыкновенных дифференциальных уравнений в Scilab предназначена функция:

y=ode(x0,t0,t,sys),

где х0 – вектор начальных условий ОДУ;

t 0 – начальная точка интегрирования;

t– вектор значений независимой переменной;

sys– имя функции, в которой исходя из вектора значенийtвычисляется матрица решенийy;

y– матрица решений (выходной параметр), первый столбец которой

y (1)– значение функции y ( x ), а второй — y (2 )–значение производной y ‘( x ).

Решение системы ОДУ формируется в матрице y, и выводится на экран в виде таблицы

В качестве примера рассмотрим решение системы ОДУ:

c начальными условиями x(0)=0, y (0)=0 на отрезке [ 0;10] и шагом 1 .

Решение системы ОДУ в Scilab начинается с создания функции sys, описывающей систему (рис.2.5.2-2). После того, как заданы начальные условия ОДУ, производится обращение в функции ode, в результате которого формируется матрица решения y. Решение системы выводится в виде таблицы и графика. (Шапка таблицы и легенда для графика!)

—> // Решение системы ОДУ —> // Загрузка сценария РИС2522 и его выполнение —>clear —> // Начальные условия —>x0 = [0; 0]; t0 = 0; t = 0:1:10; —> —> exec(‘РИС2522.sce’, 0); ans = 0. 0. 0. 1. 0.9802401 0.533358 2. 1.4096497 0.9693978 3. 1.7429464 0.6024417 4. 2.4027415 0.2586293 5. 3.3312751 0.005761 6. 4.2504071 -0.1650347 7. 4.826261 -0.2358589 8. 5.1581587 -0.2515654 9. 5.3963437 -0.2509281 10. 5.5981318 -0.2461951

Рис.2.5.2-2. Решение системы ОДУ

Рассмотрим задачу построения временной зависимости тока i ( t ) в RC-цепи.

Построить временную зависимость тока i ( t ) в неразветвленной RC-цепи (рис. 3.3-4), если цепь включается на постоянное напряжение при нулевых начальных условиях.

Неразветвленная RC -цепь Дано: E=1 В – ЭДС источника; R=1 Ом – сопротивление; С=1 Ф – емкость; t=[0, 6τ] c – временной интервал.

Рис.2.5.2-3. Построить временную зависимость тока i ( t ) в RC-цепи

источники:

http://findout.su/8×11365.html