Тема 6.5. Методы решения обыкновенных дифференциальных уравнений

Стр 1 из 2Следующая ⇒

Тема 6.5. Методы решения обыкновенных дифференциальных уравнений

Постановка задачи

Метод Эйлера

Методы Рунге-Кутты

Решение ОДУ n-го порядка

Сравнение методов решения ОДУ

Технология решения ОДУ средствами математических пакетов

Технология решения ОДУ средствами MathCad

Технология решения ОДУ в среде MatLab

6.2.7. Тестовые задания по теме «Методы решения ОДУ»

Постановка задачи

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

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

Из курса математического анализа известно, что обыкновенным называется такое дифференциальное уравнение от одной переменной, которое содержит одну или несколько производных от искомой функции y(x). В общем виде ОДУ можно представить следующим образом:

(6.5.1-1)

где х – независимая переменная, а n – порядок ОДУ.

Численные методы позволяют решить только ОДУ 1 -го порядка, поэтому в дальнейшем будем рассматривать только такие уравнения. Следует отметить, что ОДУ n-го порядка можно привести к системе из n уравнений 1 -го порядка, и при решении системы применить те же методы.

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

Вторая форма записи называется ОДУ, разрешенным относительно старшей производной.

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

Общее решение ОДУ содержит n произвольных постоянных С₁, С₂, ..., С_n и имеет следующий вид:

Общее решение ОДУ первого порядка содержит одну произвольную постоянную
y = j(x, C ) и описывает множество функций, удовлетворяющих уравнению y¢ = f(x, y )
(рис. 6.5.1-1).

Рис. 6.5.1-1

Если произвольная постоянная принимает конкретное значение С=С₀, то из общего решения ОДУ, в соответствии с теоремой Коши, получаем частное решение y = j(x, C₀), поскольку через каждую точку (x₀, y₀) в области допустимых значений проходит только одна интегральная кривая.

Теорема Коши для ОДУ 1 -го порядка звучит так:

Если в ОДУ функция y¢ = f(x, y) и ее частная производная f¢ (x, y ) определены и непрерывны в некоторой области G изменения переменных x и y, то для всякой внутренней точки ( x₀, y₀ ) этой области данное уравнение имеет единственное решение.

Значения x₀, y₀ называются начальными условиями. Для ОДУ 2 -го порядка, общее решение которого имеет две произвольные постоянные, в качестве начальных значений выступают x₀, y₀= j(x₀) и y₀¢ =φ ¢ (x₀).

При решении ОДУ точным решением является аналитическое выражение функции
y = j(x), а результатом решения ОДУ численными методами является таблица значений
y = j(x ) на некотором множестве значений аргумента х. Поэтому при постановке задачи численного решения ОДУ наряду с начальными условиями x₀, y₀ необходимо задать область решения - отрезок [a; b] и шаг изменения аргумента h.

Таким образом, численное решение ОДУ представляет собой таблицу значений искомой функции для заданной последовательности аргументов, x_i₊₁=x_i+h, i=0, 1, …, n, где

h = x_i₊₁-x_i называется шагом интегрирования.

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

а в многошаговых – нескольких, например

Начинать решение задачи Коши многошаговыми методами нельзя, поэтому начинают решение, используя всегда одношаговые методы.

Основная идея решения ОДУ одношаговыми методами сводится к разложению искомого решения y(x) в ряд Тейлора в окрестности текущей точки и его усечению. Число оставшихся членов ряда определяет порядок и, следовательно, точность метода.

Рассмотрим наиболее распространенные одношаговые численные методы решения ОДУ.

Метод Эйлера

Пусть дано уравнение

y¢ =f(x, y), (6.5.2-1)

с начальными условиями x₀, y₀ = y(x₀). Требуется найти решение данного уравнения на отрезке [a; b] (обычно x₀=а) в n точках.

На рис. 6.5.2-1 график искомой функции y(x) проходит через точку А(x₀, y₀), заданную начальными условиями.

Рис. 6.5.2-1

Разобьем отрезок на n равных частей и получим последовательность x₀, x₁, …, x_n, где x_i=x₀+i∙ h (i=0, 1, …, n), а h = (b-a)/n – шаг интегрирования.

Найдем y_i = y(x_i). Для этого проинтегрируем производную, заданную (6.5.2-1) на интервале [x₀; x₁], по формуле Ньютона – Лейбница:

Отсюда значение искомой функции в точке x₁

Примем допущение, что на интервале [x₀; x₁]производная исходной функции постоянна и равна своему значению в точке А(x₀, y₀). Тогда по формуле прямоугольников

Полученное выражение имеет наглядную геометрическую интерпретацию (рис. 6.5.2-1). Поскольку значение производной f’(x₀, y₀) = tga, то в прямоугольном треугольнике ABD Dy₀=h× tga, и, следовательно, y₁= y₀+Dy₀= y₀+h× f¢ (x₀, y₀). Таким образом, y₁ может быть найдено геометрически в результате замены искомой кривой y(x) касательной, проведенной в точке А.

Продолжая этот процесс и принимая подынтегральную функцию f(x) на соответствующем участке [x_i, x_i₊₁] постоянной и равной ее значению в начале отрезка, получим решение дифференциального уравнения в виде значений искомой функции y(x) на отрезке [a; b]. График решения представляет собой ломаную линию, которая называется ломаной Эйлера. При этом общая формула для определения очередного значения функции имеет вид:

(6.5.2-2)

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

Погрешность метода Эйлера связана с величиной шага интегрирования отношением e₁ =C₁h², где C₁ – произвольная постоянная.

Пример 6.5.2-1. Решить методом Эйлера ОДУ y¢ = 2x/y с начальными условиями x₀ = 1 и y₀= 1 на отрезке [1; 1.4] с шагом h = 0.2.

Методы Рунге-Кутты

Методы Рунге-Кутты – это группа методов, широко применяемых на практике для решения ОДУ. В этих методах при вычислении значения искомой функции в очередной точке х_i₊₁ используется информация о предыдущей точке х_i, y_i. Методы различаются объемом вычислений и точностью результата.

Порядок метода Рунге-Кутты определяется кратностью вычисления значения производной искомой функции f(x, y ) на каждом шаге. В соответствии с этим метод Эйлера является методом Рунге-Кутты первого порядка, поскольку для получения очередного значения y_i₊₁ функция f(x) вычисляется один раз в предыдущей точке х_i, y_i. В методах Рунге-Кутты более высоких порядков для вычисления очередного значения искомой функции в точке х_i₊₁ значение правой части уравнения y’= f(x, y ) вычисляется несколько раз, количество которых и определяет порядок метода.

Метод Рунге-Кутты 2-го порядка (Усовершенствованный метод Эйлера). Вычисление значения искомой функции в точке х_i₊₁ проводится в два этапа. Сначала вычисляют вспомогательную величину по методу Эйлера:

(6.5.3-1)

Затем значение производной искомой функции в точке (x_i₊₁, y_i₊₁) используется для вычисления окончательного значения функции:

(6.5.3-2)

Подставляя (6.5.3-1) в (6.5.3-2), окончательно получим расчетную формулу метода Рунге-Кутты 2 -го порядка:

(6.5.3-3)

Этот метод также называют методом прогноза и коррекций. Сначала находят грубое приближение по методу Эйлера (прогноз), а затем уточненное значение y_i₊₁ (коррекция).

В общем виде формулу (6.5.3-3) можно представить как

(6.5.3-4)

Метод Рунге-Кутты второго порядка имеет наглядную геометрическую интерпретацию (рис. 6.5.3-1). Построение проводится следующим образом: определяется пересечением перпендикуляра, восстановленного из точки x_i₊₁ c касательной L₁, проведенной к кривой y(x) в предыдущей точке (х_i, y_i). Затем в точке проводится прямая L₂ с тангенсом угла наклона, равным . Прямую проводят через точку под углом, тангенс которого находим усреднением значений тангенсов углов наклона L₁ и L₂. Прямая L проводится параллельно через точку (х_i, y_i). Ее пересечение с перпендикуляром, восстановленным из точки х_i₊₁_, и дает уточненное значение y_i₊₁.

Рис. 6.5.3-1

Погрешность метода Рунге-Кутты второго порядка связана с величиной шага интегрирования отношением e₂ =C₂h³, где C₂– произвольная постоянная.

Пример 6.5.3-1. Решить методом Рунге-Кутты второго порядка ОДУ y¢ = 2x/y с начальными условиями x₀ = 1 и y₀= 1 на отрезке [1; 1.4] и шагом h = 0.2.

Проводя дальнейшее обобщение формул Рунге-Кутты, для решения ОДУ первого порядка можно записать следующее:

где Ф – линейная функция аргументов x, y, h и f(x, y), которая может быть представлена как

(6.5.3-5)

Величина n в (6.5.3-4) определяется порядком метода, а коэффициентам a₂, a₃, …, a_n, Р₁, Р₂, …, P_n подбирают такие значения, которые обеспечивают минимальную погрешность. Так, для метода Рунге-Кутты четвертого порядка (n=4) получена расчетная формула при следующих коэффициентах: a₂= a₃=1/2, a₄=1, P₁ = P₄=1/6, P₂ = P₃ =2/6.

Подставив значения коэффициентов в (6.5.3-4), имеем

(6.5.3-6)

Геометрическая интерпретация этого метода очень сложна и потому не приводится.

Погрешность метода Рунге-Кутты четвертого порядка значительно меньше методов первого и второго порядков и пропорциональна величине h (e₄ =C₄h⁵).

Пример 6.5.3-2. Решить методом Рунге-Кутты четвертого порядка ОДУ y¢ = 2x/y с начальными условиями x₀ = 1 и y₀= 1 на отрезке [1; 1.4] с шагом h = 0, 2.

Сведем в таблицу результаты решения уравнения y¢ =2x/y методами Рунге-Кутты, соответственно, первого (y¹_i), второго (y²_i) и четвертого (y⁴_i) порядков и сравним с результатами, полученными точным методом (y_i).

х_i	y¹_i	y²_i	y⁴_i	y_i
1.2 1.4	1.4 1.74286	1.3714 1.7091	1.37115 1.7089	1.37113 1.7088

На практике для обеспечения требуемой точности (при использовании любого приближенного метода решения ОДУ ) применяется автоматический выбор шага методом двойного просчета. При этом в каждой точке х_i по формуле, соответствующей выбранному методу, производится расчет y_i с шагом h (y_i⁽^h⁾) и с шагом h/2 (y_i⁽^h^/2)). Цель двойного просчета состоит в том, чтобы для каждой точки численного решения эти значения отличались на величину, не превышающую заданной погрешности e. В этом случае общая формула для оценки погрешности решения ОДУ методами Рунге-Кутты имеет следующий вид:

где p – порядок метода Рунге-Кутты. Эта формула называется также правилом Рунге.

Если | y_i⁽^h⁾) - y_i⁽^h^/2)|< e, то шаг для следующей точки выбирается равным h, иначе шаг уменьшается вдвое и продолжается уточнение y _i в точке х_i.

Схемы алгоритмов интегрирования ОДУ методом Рунге-Кутты с автоматическим выбором шага приведены на рис. 6.5.3-2 и рис. 6.5.3-3.

Рис. 6.5.3-2. Схема алгоритма процедуры-функции решения ОДУ в очередной точке

Рис. 6.5.3-3. Схема алгоритма интегрирования ОДУ методом Рунге-Кутты с

автоматическим выбором шага

Решение ОДУ n-го порядка

Методы, рассмотренные выше, позволяют найти численное решение ОДУ только первого порядка. Однако они применимы и к уравнениям n -го порядка. Для этого ОДУ n -го порядка предварительно приводится к системе n уравнений первого порядка.

Пусть, например, требуется решить ОДУ второго порядка

с начальными условиями , , .

Обозначим z=y’. В результате подстановки в исходное уравнение получим систему двух уравнений первого порядка

с двумя неизвестными функциями и и начальными условиями

, .

В общем виде система уравнений может быть представлена в виде

(6.5.4-1)

Решением системы (6.5.4-1) являются две функции и , из которых - решение исходного уравнения второго порядка. Выбрав, например, метод Эйлера, приближенное решение системы (6.5.4-1) можно найти с помощью двух рекуррентных формул:

Пример 6.5.4-1. Дано обыкновенное дифференциальное уравнение второго порядка при начальных условиях , , на отрезке [0; 0.4] с шагом .

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

с начальными условиями , , .

Применим метод Эйлера для решения системы ОДУ

и т.д.

xi	yi	zi

0.2	1.4	1.6
0.4	1.72	1.808

В общем виде ОДУ n-го порядка

Введем следующие обозначения:

…

В результате этих подстановок перейдем к системе n ОДУ первого порядка:

(6.5.4-2)

Решением системы (6.5.4-2) являются функции

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

Окончательным решением ОДУ n-го порядка, согласно определению, служит функция , вычисленная на заданном множестве точек [a; b].

Средствами MathCad

При решении ОДУ его следует привести к нормальной форме (к виду разрешенному относительно производной исходного ОДУ )

Для ОДУ с разделяющимися переменными исходное уравнение можно привести к виду , тогда выражение задает решение задачи Коши с начальными условиями как функцию y от переменной х.

Пример 6.5.6-1. Решить ОДУ вида .

Найдем частное решение данного ОДУ с использованием средств Mathcad, сначала методом разделения переменных, а затем с использованием функции odesolve(x, xk, n), где х – имя переменной, относительно которой решается уравнение, xk – конец интервала интегрирования, n – количество шагов, на которых вычисляется решение ОДУ. Результаты подтверждают правильность преобразований.

произвольная постоянная

Аналитическое решение ОДУ

Численное решение ОДУ

Аналитическое выражение для решений ОДУ удается получить достаточно редко, поэтому широкое распространение при решении ОДУ получили численные методы.

Пример 6.5.6-2. Решить ОДУ на отрезке [0; 3] методами Рунге-Кутты с постоянным шагом h=0, 6.

В приведенном ниже документе решение, полученное методом Эйлера, обозначено как y1, методом Рунге-Кутты 2-го порядка – y2, а 4-го порядка – y4.

Метод Эйлера

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

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

В Mathcad нет средств символьного решения ОДУ, но достаточно широко представлены методы численного решения задачи Коши. Для этого предназначена, например, функция rkfixed(y, x₀, x_end, N, D), где y – первоначально равно y₀, x₀ и x_end – начальное и конечное значения аргумента, N – количество проводимых вычислений решения, а переменной D(x, y) должно быть присвоено выражение для вычисления правой части уравнения. Результатом вычислений функции rkfixe( ) служит матрица, в первом столбце которой содержатся координаты узлов x₀ … x_end, а во втором – значения приближенного решения в соответствующих узлах. В функции rkfixed( ) вместо метода Рунге-Кутты используется метод Булирша-Штера. Ниже приведены решения и их графическая иллюстрация, полученные с шагом 0.6 и 0.15.

Пример 6.5.6-3. Решить ОДУ у’= на отрезке [0; 3] методами Рунге-Кутты с постоянным шагом h=0.6.

Решение ОДУ 2-го порядка вида у”=F(x, y, z), где z=y’ также может быть получено методом Рунге-Кутты 4-го порядка. Ниже приведены формулы для решения ОДУ:

Система Mathcad имеет специальную встроенную функцию для решения дифференциальных уравнений. Она имеет вид: Odesolve ( x, b [, steps ] ).

Для решения задачи Коши необходимы так называемые начальные условия и указание конца интервала. Эти данные вместе с самим уравнением записываются в блок функции Given, и лишь затем применяется сама функция odesolve( ). Функция имеет ряд особенностей. Если указано число шагов step, то решение выполняется с фиксированным шагом, иначе - адаптивным методом.

Пример 6.5.6-4. Решить дифференциальное уравнение.

Задано дифференциальное уравнение

Заданы начальные условия

Задано решение дифференциального уравнения

Вычисление производной от b(a)

График решения заданного дифференциального уравнения b(a) и производной от функции решения -c(a)

В среде MatLab

В MatLab имеется ряд функций для решения задачи Коши. Одна из них - ode45 - использует метод Рунге-Кутты четвёртого-пятого порядка точности с автоматическим выбором размера шага.

Обращение к ode45 выполняется следующим образом

[T, Y] = ode45(ODEFUN, TSPAN, Y0).

Входные параметры процедуры ode45( ):

ODEFUN - имя функции (в виде строчной переменной), задающей правую часть системы дифференциальных уравнений. Уравнения должны быть записаны в нормальной форме

u' = ODEFUN(T, Y),

где u' - вычисляемый вектор-столбец производных;
TSPAN= [T0 TFINAL] - вектор, задающий интервал изменения независимой переменной; T0 - начальная точка;

TFINAL - последнее значение аргумента, при котором завершается расчёт;
Y0 - вектор начальных значений зависимых переменных.

Выходные параметры ode45( ):

T - вектор, содержащий отсчёты аргумента в точках решения;

Y - массив, содержащий вычисленные значения u и u' в точках, соответствующих отсчетам независимой переменной в T.

Требования к точности и другие параметры численного решения задаются в MatLab по умолчанию. Изменить эти настройки позволяет дополнительный аргумент OPTIONS.

Пример 6.5.6-5. Решить дифференциальное уравнение второго порядка,

описывающее колебания в нелинейной системе (например, автогенераторе), известное как уравнение Ван-дер-Поля:

u" + (u² - b) u' + u = 0,

где u, u' и u" - функция, её первая и вторая производные по времени t, b - произвольный параметр. Параметр b определяет потери в системе, а слагаемое u² в скобках - её нелинейные свойства (например, рабочие характеристики активного элемента автогенератора).

Данное дифференциальное уравнение второго порядка можно привести к системе двух дифференциальных уравнений первого порядка

y1' = ( 1 – (y2)² ) y1 - y2,
y2' = y1,

где y1=u' и y2=u, а параметр b принят равным единице.

Для выполнения расчётов воспользуемся внешней функцией vdp1( ), записанной в m-файле vdp1.m. Эта функция, вычисляющая правую часть уравнения Ван-дер-Поля, имеет следующий вид:

-Пример_6_5_6_5

function dydt = vdp1(t, y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)]

Для получения решения дифференциального уравнения, описываемого функцией vdp1( ) на интервале 0 < t < 20, при начальных условиях y1 = u' = 2 и y2 = u = 0 следует ввести в MatLab следующую строку:

Пример 6.5.6-5

[t, y]=ode45(@vdp1, [0 20], [2 0]).

На рисунке представлены рассчитанные функцией ode45 и построенные с помощью команды plot зависимости y2(t) ≡ u(t) и y1(t) ≡ u'(t).

График на рисунке показывает, что с течением времени колебания в системе нарастают. Однако в дальнейшем рост амплитуды колебаний замедляется, и система переходит в устойчивое состояние, описываемое представленными фазовыми кривыми y2 (y1). Кроме рассмотренной функции ode45( ) MatLab содержит ряд других функций для решения задачи Коши - ode23( ), ode113( ), ode15s( ), ode23s( ), ode23t( ), ode23tb( ).

Пример 6.5.6-6. Решить систему уравнений

с начальными условиями

Пример 6.5.6-6

Function F=FF(t, x) F=[-4*x(1)-13*x(2)+exp(t); x(1)]; end %Формирование вектора начальных условий X0=[1; -1]; interval=[0.25 2]; [T, X]=ode45(@FF, interval, X0); plot(T, X(:, 1), ': ', T, X(:, 2), '-'); legend('y', 'x-Решение'); grid on;