Типы задач для обыкновенных дифференциальных уравнений

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

В дифференциальное уравнение n-го порядка в качестве известных входят функции y(x) и ее первые n производных по аргументу x:

(5.1)

 

Из теории ОДУ известно. Что уравнение (5.1) эквивалентно системе уравнений первого порядка

(5.2)

где k =1,…, n.

 

Уравнение (5.1.) и эквивалентная ему система (5.2) имеют бесконечное множество решений. Единственные решения выделяют с помощью дополнительных условий, которым должны удовлетворять искомые решения. В зависимости от вида таких условий рассматривают три типа задач, для которых доказано существование и единственность решений.

Первый тип – это задачи Коши, или задачи с начальными условиями. Для таких задач кроме исходного решения (5.1) в некоторой точке должны быть заданны начальные условия, т.е. значение функции y(x) и ее производных:

 

 

Для системы ОДУ типа (5.2) начальные условия задаются в виде

 

(5.3)

 

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

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

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

Большинство методов решения ОДУ основано на задаче Коши, ал­горитмы для которой рассматриваются в настоящей главе.

 

Метод Эйлера

Систему ОДУ (5.2) часто удается представить в каноническом ви­де, в так называемом виде Коши:

(5.4)

где k = 1,2,..., n.

При формулировке задачи Коши система (5.4) дополняется началь­ными условиями (5.3). Для простоты рассмотрим задачу Коши для од­ного уравнения типа (5.4), а затем полученные алгоритмы обобщим на систему п уравнений:

(5.5)

В окрестности точки х₀ функцию у(х) разложим в ряд Тейлора

(5.6)

который можно применить для приближенного определения искомой функции у(х). В точке при малых значениях h можно ограни­читься двумя членами ряда (5.6), тогда

y(x₀ + h) = y₀ + hy'(x₀) + O(h²), (5.7)

где 0(h²) — бесконечно малая величина порядка h². Заменим произ­водную у'(), входящую в формулу (5.7), на правую часть уравнения (5.5):

(5.8)

Теперь приближенное решение в точке = + h можно вновь рас­сматривать как начальное условие, и по формуле (5.8) найти значение искомой функции в следующей точке х₂ = + h. В результате полу­чен простейший алгоритм решения задачи Коши, который называется методом Эйлера или методом ломаных. Последнее название связано с геометрической интерпретацией процесса (рис. 5.1); искомую функцию у(х) мы заменяем ломаной линией, представляющей собой отрезки ка­сательных к этой функции в узлах , ,...

 

Рис.5.1.Метод Эйлера

 

Формула (5.8) может быть получена и из других соображений. За­меним производную в левой части уравнения (5.5) приближенным ко­нечно-разностным отношением

Нетрудно видеть эквивалентность последнего выражения с алгорит­мом Эйлера (5.8).

На каждом шаге метода Эйлера решение у(х) определяется с по­грешностью за счет отбрасывания членов ряда Тейлора, пропорцио­нальных h в степени выше первой. Это означает, что метод Эйлера имеет второй порядок локальной погрешности. Глобальная погреш­ность, как показано в [19], имеет первый порядок и при постоянном шаге ft для оценки погрешности применима первая формула Рунге (4.14):

где — приближенное решение дифференциального уравнения в точке х, полученное с шагом h; — приближенное решение того же уравнения с шагом kh; p — порядок метода.

Фомула (5.9) позволяет опытным путем определять шаг h, обеспе­чивающий требуемую точность решения у{х). Так же, как и при вы­числении определенных интегралов, можно осуществить автоматиче­ское изменение шага в процессе интегрирования дифференциального уравнения.

Для уточнения решения применима вторая формула Рунге (4.15):

(5.10)

Формула Эйлера (5.8) обобщается для системы ОДУ, записанных в форме Коши (5.4) с начальными условиями (5.3)

(5.11)

 

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

Для уменьшения погрешности метода интегрирования ОДУ, исполь­зующего разложение искомого решения вряд Тейлора (5.6), необходи­мо учитывать большее количество членов ряда. Однако при этом воз­никает необходимость аппроксимации производных от правых частей ОДУ. Оснавная идея методов Рунге-Кутта заключается в том, что про­изводные аппроксимируются через значения функции f(x, у) в точках на интервале которые выбираются из условия наибольшей близости алгоритма к ряду Тейлора. В зависимости от старшей сте­пени h, которой учитываются члены ряда, построены вычислительные схемы Рунге-Кутта разных порядков точности.

Так, например, для второго порядка получено однопараметрическое семейство схем вида [2, 1]:

(5.12)

где 0 < £ 1 — свободный параметр,

.

Локальная погрешность схемы (5.12) имеет третий порядок, гло-* бальная — второй, т.е. решение ОДУ, полученное по этой схеме, рав­номерно сходится к точному решению с погрешностью 0(h²).

Для параметра £ чаще всего используют значение £ = 0,5 и £ = 1. В первом случае формула (5.12) приобретает вид

, (5.13)

геометрическая интерпретация которой представлена на рис.5.2,а [2, 1 ]. Вначале вычисляется приближенное решение ОДУ в точке х₀+h по формуле Эйлера Затем определяется наклон интеграль­ной кривой в найденной точке и после нахождения сред­него наклона на шаге h находится уточненное значение .. Схемы подобного типа называют "прогноз — коррекция", что подра­зумевает грубое вычисление решения по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении инте­гральной кривой.

 

Рис. 5.2. Метод Рунге-Кутта второго порядка: а) £ = О,5; б) £ = 1

 

С целью экономии памяти при программировании алгоритма (5.13), обобщенного для решения системы ОДУ, изменим его запись с учетом того, что :

(5.14)

где к — номер решения для системы ОДУ. Теперь не придется держать в памяти ЭВМ массив начальных значений , его можно не сохра­нять после вычисления значений эйлеровских приближений .Хотя по сравнению с методом Эйлера схема (5.14) требует дополнительного массива для запоминания значений

Во втором случае при £ = 1 от формулы (5.12) переходим к схеме

(5.15)

геометрический смысл которой отражает рис.5.2,6. Здесь при прогнозе определяется методом Эйлера решение в точке

,

а после вычисления наклона касательной к интегральной кривой в сред­ней точке решение корректируется по этому наклону.

Формула (5.15) обобщается для решения системы ОДУ аналогично схеме (5.14). По сравнению с программой метода Эйлера для сохране­ния начальных значений придется ввести дополнительный массив.

Схему (5.14) можно получить из метода Эйлера с помощью первой и второй формул Рунге без использования общего соотношения (5.12) для методов второго порядка.

Для построения вычислительных схем методов Рунге-Кутта дру­гих порядков в тейлоровском разложении искомого решения у(х) учи-, тываются члены, содержащие степени шага h до m включительно, где m — порядок схемы. После аппроксимации производных правой ча­сти ОДУ f(x, у) можно получить семейство схем Рунге-Кутта различ­ных порядков [20, 21, 1], наиболее используемые из которых третьего и четвертого порядка приведены ниже. Они обобщаются для решения систем ОДУ, записанных в форме Коши.

Одна из схем метода Рунге-Кутта третьего порядка имеет вид

, (5.16)

где

.

 

Также используют другую схему

(5.17)

где

 

Схемы (5.16) и (5.17) на каждом шаге вычислений требуют нахо­ждения правой части ОДУ в трех точках. Локальная погрешность схем имеет четвертый порядок, глобальной — третий.

Схема метода Рунге-Кутта четвертого порядка имеет вид

 

, (5.18)

 

 

Известна также схема метода Рунге-Кутта четвертого порядка:

(5.19)

.

 

 

Наиболее часто в вычислительной практике используют схему ме­тода Рунге-Кутта четвертого порядка, поскольку она имеет рад пре­имуществ при программной реализации

(5.20)

 

 

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

(5.21)

i = 1,2,..., n — номер уравнения в системе ОДУ из п уравнений.

Вычислительная схема метода Рунге-Кутта пятого порядка имеет вид

(5.22)

где

 

Схема (5.22) на каждом шаге требует вычисления правой части ОДУ в пяти точках. Локальная погрешность этих схем имеет шестой поря­док, глобальная — пятый.

 

Метод Рунге-Кутта-Мерсона

Мерсон предложил модификацию метода Рунге-Кутта четвертого порядка, позволяющую оценивать погрешность на каждом шаге и при­нимать решение об изменении шага. Схему Мерсона [22] с помощью эквивалентных преобразований приведем к виду, удобному для про­граммирования:

(5.23)

Схема Мерсона требует на каждом шаге вычислять правую часть ОДУ в пяти точках, но за счет только одного дополнительного коэффи­циента ki по сравнению с классической схемой Рунге-Кутта на каждом шаге можно определить погрешность решения R по формуле

(5.24)

Для автоматического выбора шага интегрирования рекомендуется следующий критерий. Если абсолютное значение величины R, вычис­ленное по формуле (5.24), больше допустимой погрешности е, то шаг h уменьшается в два раза и вычисления по схеме (5.23) повторяются с точки - Шаг h можно удвоить при выполнении условия:

32|R | < . (5.25)

Автоматический выбор шага позволяет значительно сократить вре­мя решения ОДУ. Схема (5.23) обобщается для решения системы ОДУ аналогично классической схеме Рунге-Кутта.