Дифференциальные уравнения высших порядков matlab

Решение ОДУ в Matlab

Доброго времени суток! Сегодня мы поговорим о решении ОДУ (обыкновенных дифференциальных уравнений) в Matlab. Перед тем как мы начнём обсуждать данную тему, советую вам ознакомиться с темой: Численное дифференцирование в Matlab, чтобы лучше понимать теоретическую составляющую решения ОДУ.

Обыкновенные дифференциальные уравнения

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

• Задача Коши
• Краевая задача
• Задача на собственные значения

Кратко расскажу о их сути:

Задача Коши предполагает дополнительные условия в виде значения функции в определённой точке.
Краевая задача подразумевает поиск решения на заданном отрезке с краевыми (граничными) условиями в концах интервала или на границе области.
Задача на собственные значения — помимо искомых функций и их производных, в уравнение входят дополнительное несколько неизвестных параметров, которые являются собственными значениями.

Методы решения дифференциальных уравнений

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

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

Решение обыкновенных дифференциальных уравнений в Matlab можно реализовать «своими ручками», прописав алгоритм по разным схемам. Но также в Matlab есть встроенные функции, выполняющие все стандартные задачи.

Метод Рунге-Кутта первого порядка

Методы Рунге-Кутта представляют собой разложения в ряд Тейлора и от количества использованных элементов ряда зависит порядок этого метода. Следовательно, помимо Рунге-Кутта первого порядка, вы сможете увидеть методы других порядков. Иногда их называют другими именами.

Например, Метод Рунге-Кутта первого порядка, также известен как Метод Эйлера или Метод ломаных. Информацию о его математическом и графическом представлении советую поискать в гугл. Мы же поговорим о том, как Метод Рунге-Кутта первого порядка реализуется в Matlab для решения ОДУ. Например:

Решить и привести график ошибки уравнения y’ = y*x методом Рунге-Кутта первого порядка (Методом Эйлера, Методом ломаных).

Погрешность Метода Рунге-Кутта 1 порядка

» data-medium-file=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?fit=300%2C236&ssl=1″ data-large-file=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?fit=622%2C489&ssl=1″ loading=»lazy» src=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/%D0%A0%D1%83%D0%BD%D0%B3%D0%B5-1-%D0%BF%D0%BE%D0%B3%D1%80%D0%B5%D1%88%D0%BD%D0%BE%D1%81%D1%82%D1%8C.png?resize=622%2C489″ alt=»Погрешность метода 1 порядка» width=»622″ height=»489″ srcset=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?w=629&ssl=1 629w, https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?resize=300%2C236&ssl=1 300w» sizes=»(max-width: 622px) 100vw, 622px» data-recalc-dims=»1″ />
На данном графике показана зависимость величины ошибки от шага.

Метод Рунге-Кутта второго порядка

Также известен как Метод Эйлера-Коши. Как видите, во второй части уравнения происходит обращения к следующему шагу. Но как тогда быть, если нам ещё не известен следующий шаг? Всё просто. Метод Рунге-Кутта второго порядка — это всё тот же метод первого порядка, однако, на половине шага происходит нахождение «первичного» решения, а затем происходит его уточнение. Это позволяет поднять порядок скорости сходимости до двух.

Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта второго порядка.

По сравнению с Рунге-Куттом первого порядка изначальная ошибка уже гораздо меньше.

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

Метод Рунге-Кутта четвёртого порядка

Метод Рунге-Кутта четвёртого порядка считается самым распространённым. Тем не менее, работает он аналогично второму и третьему порядку.

Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта четвёртого порядка.

Как видите, на последней картинке размерность ошибки на столько мала, что пришлось воспользоваться loglog() для лучшей видимости.

Решение ОДУ в Matlab стандартными средствами

Стоит отметить, что мы с вами разобрали только один самый известный метод решения ОДУ с разными порядками. Однако, методов очень много.

Для решения дифференциальных уравнений и систем в MATLAB предусмотрены следующие функции:

ode45 (f, interval, X0, [options])
ode23 (f, interval, X0, [options])
ode113 (f, interval, X0, [options])
ode15s (f, interval, X0, [options])
ode23s (f, interval, X0, [options])
ode23t (f, interval, X0, [options])
ode23tb (f, interval, X0, [options])

Входными параметрами этих функций являются:

• f — вектор-функция для вычисления правой части уравнения системы уравнений;
• interval — массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы;
• Х0 — вектор начальных условий системы дифференциальных уравнений;
• options — параметры управления ходом решения дифференциального уравнения или системы.

Все функции возвращают:

• массив Т — координаты узлов сетки, в которых ищется решение;
• матрицу X, i-й столбец которой является значением вектор-функции решения в узле Тi.

В функции ode45 реализован метод Рунге-Кутта 4-5 порядка точности, в функции ode23 также реализован метод Рунге-Кутта, но 2-3 порядка, а функция ode113 реализует метод Адамса.

Для решения жёстких систем предназначены функция ode15s, в которой реализован метод Гира, и функция ode23s, реализующая метод Розенброка. Для получения более точного решения жёсткой системы лучше использовать функцию ode15s. Для решения системы с небольшим числом жёсткости можно использовать функцию ode23t, а для грубой оценки подобных систем служит функция ode23tb.

Символьное решение обыкновенных дифференциальных уравнений произвольного порядка осуществляет функция dsolve r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…‘, ‘v’)
Пример использования:

На этом мы закончим. Если остались вопросы, задавайте их в комментариях. Также вы можете скачать исходники чтобы лучше понять тему: «Решение ОДУ в Matlab».

Дифференциальные уравнения высших порядков matlab

Дифференциальные уравнения и системы уравнений

Для решения дифференциальных уравнений и систем в MATLAB предусмотрены следующие функции ode45(f, interval, X0 [, options]), ode23(f, interval, X0 [, options]), ode113(f, interval, X0 [, options]), odel5s(f, interval, X0 [, options]), ode23s(f, interval, X0 [, options]), ode23t (f, interval, X0 [,options]) и ode23tb(f, interval, X0 [, options]).

Входными параметрами этих функций являются:

• f — вектор-функция для вычисления правой части уравнения системы уравнений
• interval — массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы;
• Х0 — вектор начальных условий системы дифференциальных систем
• options — параметры управления ходом решения дифференциального уравнения или системы.

Все функции возвращают:

• массив Т — координаты узлов сетки, в которых ищется решение;
• матрицу X, i-й столбец которой является значением вектор-функции решения в узле Тi

В функции ode45 реализован метод Рунге-Кутта 4-5 порядка точности, в функции ode23 также реализован метод Рунге-Кутта, но 2-3 порядка, а функция ode113 реализует метод Адамса.

В М-файле с именем pr 7. m пишем:

Потом в командном окне вызываем функцию ode113:

ode113(@pr7,[0 20],0) %Метод Адамса: @ pr 7 – ссылка на М-функцию, [0 20]- интервалы интегрирования,0 — условие: y(0)=0

Результатом будет график:

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

В М-файле с именем pr 8. m пишем:

Потом в командном окне вызываем функцию ode 45:

Библиотека функций MATLAB’а

Важным средством, позволяющим ускорить разработку программ, являются библиотеки стандартных программ. Есть такая библиотека и в MATLAB’е.

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

Что еще почитать

Как вы понимаете, на одном занятии невозможно освоить принципы отладки и тестирования в полной мере. Если вы намерены совершенствовать свое мастерство программиста, вам придется ознакомиться с более серьезными работами по этому вопросу, чем эти краткие заметки.

Надо сказать, что такой литературы, которая учила бы принципам программирования, а не синтаксическим правилам какого-либо языка программирования не так уж много. «Программистов часто учат программировать и редко – отладке» (ван Тассел).

Могу порекомендовать следующие книги, в которых вопросы отладки и тестирования подробно рассмотрены:

1. Ван Тассел Д. Стиль, разработка, эффективность, отладка и испытания программ. 1981
2. Керниган, Пайк. Практика программирования. 2004
3. Иванова Г.С. Технология программирования. 2002

Задания

Написать программу решения дифференциальных уравнений по методу Эйлера. При отладке и тестировании вести учет ошибок. Получить решение этой же задачи, используя библиотечную функцию MATLAB’а. Сравнить полученные результаты

№	Задача Коши для уравнения первого порядка
Задача Коши для уравнения высшего порядка

Справочная информация

Обыкновенное дифференциальное уравнение 1-го порядка

связывает независимую переменную t,искомую функцию y и ее производную. Решение дифференциальногоуравнениязаключается в отысканиифункции y = y(t), обращающей это уравнение в тождество на конечном или бесконечном интервале (a, b). Различают общее и частные решения дифференциального уравнения. Общее решение имеет вид y = y(t, C), где C — произвольная постоянная интегрирования. Его графическим отображением является семейство кривых (см. рис.1), называемых интегральными. Каждая интегральная кривая является отображением частного решения, соответствующего своему значению постоянной C. Для выделения частного решения из множества общего решения необходимо задать начальное условие

.

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

Для решения задачи Коши существует множество методов, которые делятся на одношаговые и многошаговые. Все они позволяют получить искомое решение в виде таблично заданной функции, в той или иной мере согласующееся с истинным частным решением (см. рис.2). Эти группы методов различаются объемом информации, которая используется для вычисления координат очередной точки табличной функции. Одношаговые методы используют значения функции и ее производной только в одной предыдущей точке, в то время как многошаговые — в нескольких. К одношаговым методам решения задачи Коши относятся метод Эйлера, модифицированный метод Эйлера, методы Рунге-Кутта и другие.

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

,

или в матричной форме следующим образом

, ,

, , .

Система дифференциальных уравнений связывает независимую переменную x, искомые функции y₁, y₂, . y_n и их первые производные. В данном случае решение задачи Коши заключается в отыскании функций y₁ = y₁(x), y₂ = y₂(x). y_n = y_n(x), обращающих каждое уравнение системы в тождество на конечном или бесконечном интервале (a, b) и удовлетворяющих начальным условиям.

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

,

, , , … , ,

, , , … , ,

,

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

Для решения такой задачи Коши используются те же методы, что для обыкновенных дифференциальных уравнений 1-го порядка. Это обуславливается тем, что матричная форма записи задачи Коши для нормальной системы полностью совпадает с ее формулировкой для этих уравнений. Единственным отличием является то, что вместо функций y(t) и f(t,y) используются вектор-функции y и f, состоящие из n функций y₁(t), y₂(t). y_n(t) и f₁(t,y₁. y_n), f₂(t,y₁. y_n). f_n(t,y₁. y_n), соответственно. При этом расчетные схемы методов и оценки их погрешностей сохраняются.

Он является старейшим методом решения задачи Коши и заключается в последовательном применении следующих формул

,

,

,

геометрическая интерпретация которых представлена на рис.3. В точке t₀ вычисляется значение производной dy/dt через f(t,y), которое определяет угол наклона касательной к графику точного решения задачи Коши. Следующая точка численного решения определяется как точка на этой касательной с абсциссой t₁ = t₀ + h. В компактном виде эти соотношения записываются следующим образом

,

.

Метод Эйлера — самый простой метод численного интегрирования. Он относится к методам первого порядка точности, поскольку его решение совпадает с истинным только в том случае, когда истинное решение является линейной функцией y = a₁ + a₂t. Его погрешность e_k на каждом шаге пропорциональна величине h 2 . Это обусловлено тем, что в качестве направления, определяющего положение следующей точки численного решения, используется касательная в крайней левой точке каждого отрезка [t_k,t_k+1]. Из рис.3 видно, что для получения более точного численного решения надо использовать некоторое промежуточное направление между направлениями касательных в крайних точках рассматриваемого отрезка.

Соотношения метода Эйлера для нормальной системы в матричной форме имеют вид

,

.

или в развернутой форме

, .

Геометрическая интерпретация работы метода Эйлера решения задачи Коши для нормальной системы идентична его геометрической интерпретации для дифференциальных уравнений 1-го порядка. Однако в данном случае движение осуществляется вдоль некоторойгиперкривойв (n+1)-мерном пространстве переменных t, y₁, y₂. y_n.

Программное обеспечение

Ниже приведен текст подпрограммы на алгоритмическом языке FORTRAN, реализующей метод Эйлера:

Subroutine Eu(n, dt, k, Fun, t, Y)

В качестве параметров в подпрограмме используются:

n — порядок системы дифференциальных уравнений;

dt — длина отрезка интегрирования уравнений;

k — количество шагов на длине отрезка интегрирования;

Fun — имя внешней подпрограммы типа Subroutine, с помощью которой вычисляются значения вектора правой части f(t, y)системы;

t — значение аргумента системы в начале отрезка интегрирования при обращении к подпрограммам и в конце этого отрезка (t + dt)после отработки подпрограмм;

Y — массив (n элементов) значений вектора решения y(t)системы уравнений в начале отрезка интегрирования при обращении к подпрограммам и в конце этого отрезка после отработки подпрограмм.

Подпрограмма Fun оформляется в виде

Subroutine Fun(n, t, Y, F)

Dimension Y(n), F(n)

Здесь F — массив значений вектора правой части f(t, y) системы уравнений.

Subroutine Eu(n, dt, k, Fun, t, Y)

Call Fun(n, t, Y, F1)

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

· [t,Y] = ode23(fun, tspan, y0) – решает задачу Коши для системы обыкновенных дифференциальных уравнений с использованием комбинации методов Рунге-Кутта относительно невысокого порядка (2-го и 3-го);

· [t,Y] = ode45(fun, tspan, y0) — решает задачу Коши для системы обыкновенных дифференциальных уравнений с использованием комбинации методов Рунге-Кутта более высокого порядка (4-го и 5-го).

Входные аргументы указанных функций:

fun — функция для вычисления вектора правой части системы, интерфейс которой должен иметь вид:

где t – независимый параметр системы, Y- вектор значений неизвестных функций, а выходной параметр F – вектор вычисленных компонент правой части;

tspan — вектор, определяющий интервал интегрирования системы. Если этот вектор состоит из двух компонент, то их значения задают начало и конец интервала. Если требуется, чтобы в качестве решения были выданы значения искомых функций в конкретно заданных точках интервала, то число компонент должно быть не менее трех (можно задавать в виде имени массива или при помощи конструктора массивов, например в виде [2:0.1:5]);

y0 — вектор начальных значений.

t — вектор, состоящий из значений независимой переменной, которым соответствуют вычисленные значения решения, помещенные в массив Y (если вектор tspan состоит лишь из двух компонент, то количество значений автоматически определяется функцией с целью наиболее точного отображения на графике);

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

Параметр fun, являющийся фактическим параметром, указывающим предварительно написанную функцию для вычисления правой части системы, задается в виде:

@fname , где fname — имя m-функции.

Точность получаемого решения по умолчанию регулируется условием обеспечения относительной погрешности не более 10 -3 или абсолютной погрешности не более 10 -6 . Если требуется задать иные требования точности, то функции следует вызывать с дополнительным аргументом options:

[t,Y] = ode23(fun, tspan, y0, options) или

[t,Y] = ode45(fun, tspan, y0, options) .

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

где namei — название управляющего параметра, а valuei — новое присваиваемое ему значение. Среди управляющих параметров имеются RelTol и AbsTol , задающие относительную и абсолютную погрешности решения систем. Так, например, если предполагается, что вызываемая функция должна обеспечить относительную погрешность решения не более 10 -5 или абсолютную не более 10 -8 , то следует выполнить оператор

а затем использовать переменную options в качестве фактического аргумента при решении системы дифференциальных уравнений.