Тема 6. Численное решение дифференциальных уравнений 6.1. Постановка задачи Коши Известно, что обыкновенное дифференциальное уравнение первого порядка имеет вид: y' (t) = f(t, y(t)). (6.1) Решением уравнения (6.1) является дифференцируемая функция y(t), которая при подста- новке в уравнение (6.1) обращает его в тождество. На рис. 6.1 приведен график решения диффе- ренциального уравнения (6.1). График решения дифференциального уравнения называется инте- гральной кривой. Рис. 6.1 Производную y'(t) в каждой точке (t, y) можно геометрически интерпретировать как тангенс угла наклона касательной к графику решения, проходящего через эту точку, т е.: k = tg = f(t, y). Уравнение (6.1) определяет целое семейство решений. Чтобы выделить одно решение, за- дают начальное условие: y(t 0 ) = y 0 , (6.2) где t 0 – некоторое заданное значение аргумента t, а y 0 – начальное значение функции. Задача Коши заключается в отыскании функции y = y(t), удовлетворяющей уравнению (6.1) и начальному условию (6.2). Обычно определяют решение задачи Коши на отрезке, расположенном справа от начального значения t 0 , т. е. для t [t 0 , T]. Разрешимость задачи Коши определяет следующая теорема. Теорема 6.1. Пусть функция f(t, y) определена и непрерывна при t 0 t T, - < y < и удовлетворяет условию Липшица: | f(t, y 1 ) – f(t, y 2 )| L| y 1 – y 2 |, где L некоторая постоянная, а y 1 , y 2 – произвольные значения.
Тогда для каждого начального значения y 0 существует единственное решение y(t) задачи Коши для t [t 0 , T]. Даже для простых дифференциальных уравнений первого порядка не всегда удается полу- чить аналитическое решение. Поэтому большое значение имеют численные методы решения. Чис- ленные методы позволяют определить приближенные значения искомого решения y(t) на некото- рой выбранной сетке значений аргумента t i , (i = 0, 1, …). Точки t i называются узлами сетки, а вели- чина h i = t i+1 – t i – шагом сетки. Часто рассматривают равномерные сетки, для которых шаг h i по- стоянен, h i = h = n t T 0 . При этом решение получается в виде таблицы, в которой каждому узлу сетки t i соответствуют приближенные значения функции y(t) в узлах сетки y i y(t i ). Численные методы не позволяют найти решение в общем виде, зато они применимы к ши- рокому классу дифференциальных уравнений. Сходимость численных методов решения задачи Коши. Пустьy(t) – решение задачи Коши. Назовем глобальной погрешностью (или просто погрешностью) численного метода функ- цию i = y(t i ) – y i , заданную в узлах сетки t i . В качестве абсолютной погрешности примем величину R = n i 0 max | y(t i ) – y i | Численный метод решения задачи Коши называется сходящимся, если для него R 0 при h 0. Говорят, что метод имеет p-ый порядок точности, если для погрешности справедлива оценка R Ch p , p > 0, C – константа, C 0. 6.2. Метод Эйлера Простейшим методом решения задачи Коши является метод Эйлера. Будем решать задачу Коши y' (t) = f(t, y(t)). y(t 0 ) = y 0 , на отрезке [t 0 , T]. Выберем шаг h = n t T 0 , и построим сетку с системой узлов t i = t 0 + ih, i = 0, 1, …, n. В методе Эйлера вычисляются приближенные значения функции y(t) в узлах сетки : y i y(t i ). Заменив производную y' (t) конечными разностями на отрезках [t i , t i+1 ], i = 0, 1, …, n – 1, по- лучим приближенное равенство: h y y i i 1 = f(t i , y i ), i = 0, 1, …, n – 1, которое можно переписать так: y i+1 = y i + h f(t i , y i ), i = 0, 1, …, n – 1. (6.3) Формулы (6.3) и начальное условие (6.2) являются расчетными формулами метода Эйле- ра. Геометрическая интерпретация одного шага метода Эйлера заключается в том, что решение на отрезке [t i , t i+1 ] заменяется касательной y = y' (t i )( t - t i ), проведенной в точке (t i , y(t i )) к интеграль- ной кривой, проходящей через эту точку. После выполнения n шагов неизвестная интегральная кривая заменяется ломаной линией (ломаной Эйлера). Оценка погрешности. Для оценки погрешности метода Эйлера воспользуемся следующей теоремой. Теорема 6.2. Пусть функция f удовлетворяет условиям: дy дf K, dt df = f дy дf дt дf L. (6.4) Тогда для метода Эйлера справедлива следующая оценка погрешности:
R = ni 0 max | y( ti) – yi| KLenLl2 2 = KLehl2 2 , где l – длина отрезка [ t0 , T]. Мы видим , что метод Эйлера имеет первый порядок точности. Оценка погрешности метода Эйлера часто бывает затруднительна, так как требует вычис- ления производных функции f( t, y( t)). Грубую оценку погрешности дает правило Рунге (правило двойного пересчета), которое используется для различных одношаговых методов, имеющих p-ый порядок точности. Правило Рунге заключается в следующем. Пусть y2 / hi – приближения, получен- ные с шагом 2 h, а yhi – приближения, полученные с шагом h. Тогда справедливо приближенное равенство: | y2 / hi- y( ti)| 1 2 1 p| y2 / hi- yhi| . (6.5) Таким образом, чтобы оценить погрешность одношагового метода с шагом 2 h, нужно найти то же решение с шагом h и вычислить величину, стоящую справа в формуле (6.5), т е. R 1 2 1 p| y2 / hi- yhi| (6.6) Так как метод Эйлера имеет первый порядок точности, т. е. p = 1, то приближенное равен- ство (6.6) примет вид R | y2 / hi- yhi| (6.7) Используя правило Рунге, можно построить процедуру приближенного вычисления реше- ния задачи Коши с заданной точностью . Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y2 / hi, i = 0, 1, …, n. Вычисления прекращаются тогда, когда будет выполнено условие: R 1 2 1 p| y2 / hi- yhi| < . (6.8) Для метода Эйлера условие (6.8) примет вид R | y2 / hi- yhi| < (6.9) Приближенным решением будут значения y2 / hi, i = 0, 1, …, n. Пример 6.1. Найдем решение на отрезке [0, 1] следующей задачи Коши: y' ( t) = y – yt2 , (6.10) y(0) = 1. Возьмем шаг h = 0.2. Тогда n = 2 0 0 1 = 5. В соответствии с (6.3) получим расчетную формулу метода Эйлера: yi+1 = yi + 0.2 iiiyty2 , y0 = 1, i = 0, 1, 2, 3, 4, 5. Решение представим в виде таблицы 6.1: Таблица 6.1 i 0 1 2 3 4 5 ti 0 0.2 0.4 0.6 0.8 1.0 yi 1.0000 1.2000 1.3733 1.5294 1. 6786 1.8237 Уравнение (6.10) есть уравнение Бернулли. Его решение можно найти в явном виде: y = 1 2 t . (6.11) Для сравнения точного и приближенного решений представим точное решение (6.11) в виде таблицы 6.2: Таблица 6.2 i 0 1 2 3 4 5 t i 0 0.2 0.4 0.6 0.8 1.0 y(t i ) 1.0000 1.1832 1.3416 1.4832 1. 6124 1.7320 Из таблицы видно, что погрешность составляет R = 5 0 max k | y(t i ) – y i | = 0.0917. 6.3. Модифицированные методы Эйлера Первый модифицированный метод Эйлера. Суть этого метода состоит в следующем. Сначала вычисляются вспомогательные значения искомой функции y 2 1 i в точках t 2 1 i = t i + 2 h с помощью формулы: y 2 1 i = y i + 2 h f i = y i + 2 h f(t i , y i ). Затем находится значение правой части уравнения (6.1) в средней точке f 2 1 i = f(t 2 1 i , y 2 1 i ) и затем полагается y i+1 = y i + h f 2 1 i , i = 0, 1, …, n – 1. (6.12) Формулы (6.12) являются расчетными формулами первого модифицированного метода Эй- лера. Первый модифицированный метод Эйлера является одношаговым методом со вторым по- рядком точности Второй модифицированный метод Эйлера – Коши. Суть этого метода состоит в следую- щем. Сначала вычисляются вспомогательные значения 1
i y = y i + h f(t i , y i ). (6.13) Затем приближения искомого решения находятся по формуле: y i+1 = y i + 2 h [f(t i , y i ) + f(t i+1 , 1
i y )], i = 0, 1, …, n – 1. (6.14) Формулы (6.14) являются расчетными формулами второго модифицированного метода Эйлера – Коши. Второй модифицированный метод Эйлера – Коши, так же, как и первый, является одноша- говым методом со вторым порядком точности. Оценка погрешности. Приближенная оценка погрешности модифицированных методов Эйлера осуществляется как и для простого метода Эйлера с использованием правила Рунге (см. предыдущий раздел 6.2). Так как оба модифицированных метода Эйлера имеют второй порядок точности, т. е. p = 2, то оценка погрешности (6.6) примет вид R 3 1 |y 2 / h i - y h i |. (6.15)
Используя правило Рунге, можно построить процедуру приближенного вычисления реше- ния задачи Коши модифицированными методами Эйлера с заданной точностью . Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y 2 / h i , i = 0, 1, …, n. Вычисления прекращаются тогда, когда будет выполнено условие: R 3 1 |y 2 / h i - y h i | < . (6.16) Приближенным решением будут значения y 2 / h i , i = 0, 1, …, n. Пример 6.2. Применим первый модифицированный метод Эйлера для решения задачи Коши y' (t) = y – y t 2 , y(0) = 1, рассмотренной ранее в примере 6.1. Возьмем шаг h = 0.2. Тогда n = 2 0 0 1 = 5. В соответствии с (6.3) получим расчетную формулу первого модифицированного метода Эйлера: y i+1 = y i + h f 2 1 i = y i + 0.2 f 2 1 i , где f 2 1 i = f(t 2 1 i , y 2 1 i ) = y 2 1 i – 2 1 2 1 2 i i y t , t 2 1 i = t i + 2 h = t i + 0.1, y 2 1 i = y i + 2 h f(t i , y i ) = y i +0.1 i i i y t y 2 , t 0 = 0, y 0 = 1, i = 0, 1, …, 4. Решение представим в виде таблицы 6.3: Таблица 6.3 i t i y i 2 h f(t i , y i ) t 2 1 i y 2 1 i h f 2 1 i 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1.0 1 1.1836 1.3426 1.4850 1.6152 1.7362 0.1 0.0850 0.0747 0.0677 0.0625 0.1 0.3 0.5 0.7 0.9 1.1 1.2682 1.4173 1.5527 1.6777 0.1836 0.1590 0.1424 0.1302 0.1210 Третий столбец таблицы 6.3 содержит приближенное решение y i , i = 0, 1, …, 5. Сравним полученное приближенное решение с точным решением (6.11), представленном в таблице 6.2. Виднм, что погрешность составляет R = 5 0 max k | y(t i ) – y i | = 0.0042. Пример 6.3. Применим второй модифицированный метод Эйлера – Коши для решения задачи Коши
y' (t) = y – y t 2 , y(0) = 1, рассмотренной ранее в примерах 6.1 и 6.2. Так же, как и ранее, зададим шаг h = 0.2. Тогда n = 2 0 0 1 = 5. В соответствии с (6.14) получим расчетную формулу метода Эйлера – Коши: y i+1 = y i + 2 h [f(t i , y i ) + f(t i+1 , 1
i y )] = y i + 0.1[f(t i , y i ) + f(t i+1 , 1
i y )], где f(t i , y i ) = y i – i i y t 2 1
i y = y i + h f(t i , y i ) = y i + 0.1 i i i y t y 2 t 0 = 0, y 0 = 1, i = 0, 1, …, 4. Решение представим в виде таблицы 6.4: Таблица 6.4 i t i y i 2 h f(t i , y i ) t i+1 1
i y f(t i+1 , 1
i y ) 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1.0 1 1.1867 1.3484 1.4938 1.6272 1.7542 0.1 0.0850 0.0755 0.0690 0.0645 0.2 0.4 0.6 0.8 1.0 1.2 1.3566 1.4993 1.6180 1.7569 0.867 0.767 0.699 0.651 0.618 Таблица 6.4 заполняется последовательно по строкам, сначала первая строка, затем вторая и т. д. Третий столбец таблицы 6.4 содержит приближенное решение y i , i = 0, 1, …, 5. Сравним полученное приближенное решение с точным решением (6.11), представленном в таблице 6.2. Видим, что погрешность составляет R = 5 0 max k | y(t i ) – y i | = 0.0222. 6.4. Метод Рунге – Кутта Метод Рунге – Кутта является одним из наиболее употребительных методов высокой точно- сти. Метод Эйлера можно рассматривать как простейший вариант метода Рунге – Кутта. Рассмотрим задачу Коши для дифференциального уравнения y' (t) = f(t, y(t)) с начальным условием y(t 0 ) = y 0. Как и в методе Эйлера, выберем шаг h = n t T 0 и построим сетку с системой узлов t i = t 0 + ih, i = 0, 1, …, n. Обозначим через y i приближенное значение искомого решения в точке t i Приведем расчетные формулы метода Рунге – Кутта четвертого порядка точности: y i+1 = y i + 6 1 h(k 1 i + 2k 2 i + 2k 3 i + k 4 i ), k 1 i = f(t i , y i ),
k2 i = f( ti + 2 h, yi + 2 hk1 i), (6.17) k3 i= f( ti + 2 h, yi + 2 hk2 i), k4 i = f( ti + h, yi + hk3 i), i = 0, 1, …, n. |