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

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

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

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

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

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

(1.5.1-1)

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

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

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

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

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

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

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

Рис. 1.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), (1.5.2-1)

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

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

Рис. 1.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). Для этого проинтегрируем производную, заданную (1.5.2-1) на интервале [x₀; x₁], по формуле Ньютона – Лейбница:

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

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

Полученное выражение имеет наглядную геометрическую интерпретацию (рис. 1.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]. График решения представляет собой ломаную линию, которая называется ломаной Эйлера. При этом общая формула для определения очередного значения функции имеет вид:

(1.5.2-2)

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

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

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

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

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

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

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

(1.5.3-1)

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

(1.5.3-2)

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

(1.5.3-3)

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

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

(1.5.3-4)

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

Рис. 1.5.3-1

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

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

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

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

(1.5.3-5)

Величина n в (1.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.

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

(1.5.3-6)

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

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

Пример 1.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.

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

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

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

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

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

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

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

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

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

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

, .

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

(1.5.4-1)

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