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

1 ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ 15.1. Постановка задачи В классическом анализе разработано немало приемов нахождения решений дифференциальных уравнений через элементарные (или специальные) функции. Между тем весьма часто при решении практических задач эти методы оказываются либо совсем неприменимыми, либо их решение связывается с недопустимыми затратами усилий и времени. По этой причине для решения задач практики созданы методы приближенного решения дифференциальных уравнений. В зависимости от формы представления решения, эти методы подразделяются на три основные группы. 1. Аналитические методы, применение которых дает решение дифференциального уравнения в виде аналитического выражения. 2. Графические методы, дающие приближенное решение в виде графика. 3. Численные методы, когда искомая функция получается в виде таблицы. Остановимся на численных методах. Решим простейшее обыкновенное дифференциальное уравнение первого порядка, разрешенное относительно первой производной у= f(х,у). (15.1) Общее решение уравнения (15.1) имеет вид у = у(х,С ), (15.2) где С - произвольная постоянная. Геометрически общее решение (15.2) представляет собой семейство интегральных кривых, т.е. совокупность линий, соответствующих различным значениям постоянной С. Интегральные кривые обладают тем свойством, что в каждой точке М(х,у) наклон касательной удовлетворяет условию dy tg  = = f(x,y) dx Таким образом, дифференциальное уравнение в разрешенном относительно производной виде устанавливает явную связь между координа- 2 тами точки М(х,у) и угловым коэффициентом касательной к интегральной кривой в этой точке. Если функция f(х,у) определена на некоторой области D плоскости, то каждой точке М(х,у)D соответствует некоторое направление, угловой коэффициент которого равен f(х,у). Указывая это направление единичным вектором, проходящим через точку М, мы получим на D поле направлений. Интегральные кривые (15.2) суть кривые, для которых упомянутые направления являются направлениями касательных. Следовательно, с геометрической точки зрения задача решения дифференциального уравнения (15.1) заключается в нахождении кривых, направление касательных к которым в каждой точке совпадает с направлением поля. Конечно, в данном случае интегральные кривые принадлежат области D. Если задать точку M(х,у), через которую должна проходить интегральная кривая, то тем самым из бесконечного семейства интегральных кривых, в простейшем случае, выделяется некоторая определенная интегральная кривая, которая соответствует частному решению дифференциального уравнения (15.1). Аналитически это требование сводится к так называемому начальному условию: у = у при х = х. Задача Коши (начальная задача): найти решение уравнения у= f(х,у) в виде функции у(х), удовлетворяющей начальному условию: у(х) = у, (15.3) 3 т.е. принимающей при х = х заданное значение у = у. Геометрически задача Коши формулируется так: найти интегральную кривую дифференциаль- ного уравнения у= f(х,у), проходящую через заданную точку M(х,у). Замечание: Дифференциальные уравнения являются математическим аппаратом, с помощью которого мы можем изучать процессы, протекающие в природе. Если условия задачи полностью определяют процесс, то он должен протекать однозначно, т.е. решение дифференциального уравнения, дающее закон протекания процесса, должно быть единственным. Общее решение дифференциального уравнения содержит произвольные постоянные и, следовательно, не дает определенного ответа на поставленный вопрос. Поэтому при решении конкретных задач, кроме дифференциального уравнения, нужны еще дополнительные условия. В простейшем случае, это начальные условия, и мы приходим к задаче Коши. В курсе математического анализа рассматривалась теорема о существовании и единственности решения задачи Коши. Теорема: Пусть функция f(х,у) определена и непрерывна в области D,   D =  x  x  a , y  y  b  и имеет в ней ограниченную производную fу(х,у) , т.е. при (х,у)D f  N. y 1 b   Тогда на отрезке x  x  l, где l < min  a, , , M = max f(х,у) N M   (х,у)D существует и притом единственное решение уравнения у= f(х,у), удовлетворяющее начальному условию у(х) = у. При этом выполняется неравенство y(х)  y  b , при x  x  l . На рисунке в плоскости х0у изображен прямоугольник D и при  надлежащий ему прямоугольник D 1 =  x  x  a , y  y  b  4 Теорема утверждает, что если на прямоугольнике D функция f(х,у) f непрерывна и имеет ограниченную производную , то через точку y (х,у) проходит единственная интегральная кривая у = у(х ), определенная для всех значений х  [ x  l, x  l ]. Она полностью принадлежит прямоугольнику D 1. Подчеркнем, что теорема гарантирует существование определенного отрезка L = [ x  l, x  l ] , на котором заведомо существует и единственное решение у = у(х ) уравнения (1), проходящее через точку (х,у). В дальнейшем мы будем искать приближенное решение именно на этом отрезке, потому что нельзя ручаться, что указанное решение определено вне L. 15.2. Метод Эйлера Классы уравнений, для которых разработаны методы получения точных решений, сравнительно узки и охватывают только малую часть возникающих на практике задач. Общего же метода для нахождения точного решения произвольного дифференциального уравнения первого порядка не существует. Поэтому важное значение приобретают приближенные методы решений дифференциальных уравнений. Мы рассмотрим простейший из них, так называемый метод Эйлера для решения задачи Коши. Пусть дано уравнение у= f(х,у) с начальным условием у(х) = у. Найдем таблицу значений функции на отрезке [a,b] с шагом табулирования h. Будем предполагать, что на отрезке [a,b], где ищется решение, выполнены все условия, обеспечивающие существование и единственность решения задачи Коши. Выбрав достаточно малый шаг h, построим, начиная с точки х= a, систему равностоящих точек хi = х+ ih, ( i = 0, 1, 2, ...). При достаточно малом h производную в каждой точке х0, х1, ... , хn1 заменим на отношение приращения функции к приращению аргумента: y y y хi i+1 i h Подставим приближенное значение y хi в уравнение (15.1): yi+1  yi = f(хi,уi) h Выразим yi+1 : yi+1 = yi + h· f(хi,уi) Вместе с начальным условием у(х) = у получены расчетные фор- 5 мулы метода Эйлера для вычисления приближенных значений функции у(х) в точках х0, х1, ... , хn  у= у(х)  (15.4) yi+1 = yi + h· f(хi,уi) , ( i = 0, 1, ... , n 1)  Расчетные формулы метода Эйлера (15.4) можно получить и из иных соображений. Проинтегрируем обе части уравнения (15.1) на отрезке [х i , х i +1] х i+1 х i+1 уdx =   f(х,у) dx хi хi х i+1 y i+1  y i =  f(х,у) dx хi х i+1 y i+1 = y i +  f(х,у) dx (15.5) хi Вычислим приближенное значение определенного интеграла х i+1  f(х,у(х)) dx хi по методу левых прямоугольников. Получим, что х i+1  f(х,у(х)) dx = hf(х i,у i) хi Следовательно, yi+1 = yi + h· f(хi,уi). Геометрический смысл метода Эйлера Пусть на заданном отрезке [a,b] требуется найти решение дифференциального уравнения первого порядка у= f(х,у), удовлетворяющее начальному условию: у(х) = у. Геометрически это значит, что для дифференциального уравнения (1) нужно построить интегральную кривую у = у(х ), проходящую через точку M(х,у). Из геометрического смысла производной получаем, что в каждой точке М(х,у) интегральной кривой ее наклон ( т.е. угловой коэффициент касательной ) удовлетворяет условию dy k = tg  = = f(x,y). dx Напомним расчетные формулы метода Эйлера: 6  у= у(х)  yi+1 = yi + h· f(хi,уi)  ( i = 0, 1, ... , n 1)  xi+1 = xi + h Каждое последующее значение функции, начиная с у1 , получается путем добавления к предыдущему значению приращения yi = h· f(хi,уi) = h· tgi , зависящее от наклона интегральной кривой в точке Mi(хi,уi). у= у(х)  yi = h· f(хi,уi) = h· tgi   ( i = 0, 1, ... , n 1) yi+1 = yi + yi  xi+1 = xi + h Таким образом, по методу Эйлера интегральная кривая M0M1M2...Mn, где Mi(хi,у(xi)) ( i = 0,1,,n) заменяется ломаной M0N1N2...Nn с вершинами Ni(хi,уi) ( i = 0,1,,n) называемой ломаной Эйлера. Первое звено ломаной касается искомой интегральной кривой в точке M0(х0,у0). Все последующие звенья имеют такие же наклоны, что и касательные к интегральной кривой в точках Mi(хi,у(xi)) ( i = 1, 2, ... , n 1 ). M0N1N2...Nnграфик кусочно-линейной функции уn(x), определяемой на [х0,b] с помощью равенств (3) при разбиении отрезка на n частей. 7 Можно доказать1, что если существует единственное решение y(x) уравнения (15.1), удовлетворяющее начальному условию и определенное на отрезке [х0,b], то последовательность ломаных Эйлера {уn(x) }сходится на [х0,b] к истинному решению задачи Коши при n . Думаем, что теперь уже очевидна возможность вывода формул метода Эйлера и из чисто геометрических соображений. Физический смысл метода Эйлера С механической точки зрения мы непрерывный процесс, описываемый дифференциальным уравнением (15.1), заменяем импульсным процессом, протекающим с постоянной скоростью на элементарных промежутках [х i , х i +1] ( i = 0, 1, ... , n 1 ), скорость которого меняется скачками при переходе к последующему промежутку. Недостатки метода Эйлера: 1. Малая точность при значительном шаге h, большой объем работы при малом шаге. 2. Систематическое накопление ошибок, т.к. при вычислении значений на последующих отрезках исходные данные не являются точными и содержат погрешности, зависящие от неточности предшествующих вычислений. Оценка точности Метод Эйлера обладает малой точностью (в пределах одной-двух значащих цифр), к тому же погрешность каждого шага систематически возрастает. Наиболее приемлемым для практики методом оценки точности h является способ двойного счета - с шагом h и с шагом . Если у i - при2 ближенное значение решения, полученное при расчете с шагом h, у i* h улучшенное значение, полученное при шаге , и у(xi) - точное значение 2 решения, то абсолютную погрешность определяют из приближенного равенства. y i  у(xi)  y i  у i* , ( i = 1, 2, ... , n ) Совпадение десятичных знаков в полученных двумя способами результатах дает естественные основания считать их верными. 15.3. Модифицированный метод Эйлера Модифицированный метод Эйлера при практически том же объеме 1 Доказательство см., например, в книге: Петровский И.Г. Лекции по теории обыкновенных дифференциальных уравнений.- М.: Наука, 1970. 8 вычислительной работы дает погрешность порядка h2, вместо h в обычном методе Эйлера, что достигается с помощью очень простого приема. Возвратимся к формуле (15.5). При получении из нее формулы метода Эйлера мы полагали подынтегральную функцию f(х,у(х)) постоянной, равной ее значению f(хi,уi) в левом конце отрезка [х i , х i+1]. Более точное х i+1 значение получится, если при вычислении интеграла  f(х,у(х)) dx примехi нить метод центральных прямоугольников, полагая f(х,у(х)) равной значеh нию в центральной точке отрезка х i+  = х i + 2  х i+1  f(х,у(х)) dx = hf( х i+  , у i+  )   хi yi+1 = yi + hf( х i+  , у i+  )   Но для вычислений по этой формуле необходимо знать значение у i+  = у(х i+  ). Это значение можно найти приближенно с помощью   h обычного метода Эйлера с шагом . 2 Таким образом, получим следующую последовательность вычислений h у i+  = yi + f( х i, уi) 2  h y i+1 = yi + hf(х i + , у i+  ). 2  Геометрический смысл модифицированного метода Эйлера Модификации метода Эйлера применяются для того, чтобы более точно определить направление перехода из точки (хi,уi) в точку (хi+1,уi+1). В модифицированном методе Эйлера сначала делаем половинный шаг по методу Эйлера, вычисляя промежуточные значения h у i+  = yi +  f( х i, уi ) 2  h х i+  = х i + 2     Затем в найденной средней точке Мi+  ( х i+  , у i+  ) определяем    наклон интегральной кривой 9 k = у i+  = f( х i+  , у i+  ) = tg i+      По этому наклону определяем приращение функции на целом шаге yi = h· f( х i+  , у i+  ) = h· tg i+     И окончательно yi+1 = yi + yi Таким образом, интегральная кривая M0M1M2Mn, где Mi(хi,у(xi)) ( i = 0,1, ... ,n) заменяется ломаной M0N1N2...Nn. Каждое звено ломаной на отрезке [хi , х i +1] (i=0,1,…,n–1) имеет то же направление, что и касательh ная к интегральной кривой в средней точке х i+  = х i + отрезка. 2  Обосновать модифицированный метод Эйлера можно с помощью идеи графического построения решения дифференциального уравнения. Пусть для уравнения у= f(х,у) с начальным условием у(х) = у выполнены на отрезке [х0,b] все условия существования и единственности решения. Выберем шаг h для разбиения отрезка системой точек хi = х+ ih , ( i = 0,1,2, ... ,n) настолько малым, что для всех х между хi и х i +1 значения функции у мало отличались от линейной функции. На каждом участке [хi , хi+1] интегральная кривая заменяется отрезком прямой 10 линии, параллельной касательной, проведенной к графику решения задачи h Коши в средней точке х i+  = х i + отрезка. 2  15.4. Метод Эйлера-Коши Другой модификацией метода Эйлера является метод Эйлера-Коши. Его расчетные формулы легко получить, применив при вычислении интеграла х i+1  f(х,у(х)) dx хi в уравнении (15.5) метод трапеций. х i+1 f(х i,у i) + f(х i+1,у i+1)  f(х,у(х)) dx = h  2 хi Следовательно, f(х i,у i) + f(х i+1,у i+1) у i+1 = у i + h  2 В этой формуле в правой части участвует у i+1 . Можно сначала вычислить “грубое приближение” решения по методу Эйлера y*i+1 = yi + h· f(хi,уi), а затем уточнить f(х i,у i) + f(х i+1, y*i+1) у i+1 = у i + h  2 Таким образом, мы получаем расчетные формулы метода ЭйлераКоши: 11 у= у(х) y*i+1 = yi + h· f(хi,уi) f(х i,у i) + f(х i+1, y*i+1) у i+1 = у i + h  2 xi+1 = xi + h   ( i = 0, 1, ... , n 1)  Геометрический смысл метода Эйлера-Коши В этой модификации направление перехода из точки (хi,уi) в точку (хi+1,уi+1) определяется другим способом. На каждом отрезке [х i , х i+1] находим наклон интегральной кривой в точке (хi,уi), также грубо вычисляем по методу Эйлера значение y*i+1 = yi + h· f(хi,уi), а затем определяем наклон интегральной кривой в новой точке (хi+1, y*i+1) y*i+1 = f( х i+1 , y*i+1 ) = tg i+1 . Затем находим среднее арифметическое наклонов в точках (хi,уi) и * (х i+1, y i+1) и по среднему наклону определяем приращение функции на f(х i,у i) + f(х i+1, y*i+1) tg i + tg i+1 всем отрезке у i = h  =h , уточняя 2 2 значение у i+1 у i+1 = у i + у i 12 Таким образом, интегральная кривая M0M1M2...Mn, где Mi(хi,у(xi)) ( i = 0, 1, ... , n ) заменяется ломаной M0N1N2...Nn. Каждое звено ломаной на отрезке [хi , х i+1] ( i = 0, 1, ... , n 1 ) имеет усредненное направление между направлениями касательных к интегральной кривой на концах отрезка. 15.5. Метод Рунге-Кутта Метод Эйлера и метод Эйлера-Коши отыскания приближенного решения задачи Коши являются представителями метода Рунге-Кутта 1-го и 2-го порядка. Расчетные формулы метода Эйлера представим в виде: yi1  yi  r1 , где r1  h  f xi , yi  . А расчетные формулы метода Эйлера-Коши – в виде: 1 yi 1  yi  r1  r2 , где 2 r1  h  f  xi , yi , r2  h  f  xi  h, yi  r1 . Чем выше порядок формул Рунге-Кутта, тем меньше погрешность получаемого приближенного решения задачи Коши. На практике соблюдают некоторый компромисс между высоким порядком формул и громоздкостью вычислений. Одна из самых популярных формул семейства РунгеКутта – 4-го порядка, которую чаще называют просто формулой РунгеКутта: 13 r1  h  f  xi , yi , h r  r2  h  f  xi  , yi  1 , 2 2  h r   r3  h  f  xi  , yi  2 , 2 2  r4  h  f  xi  h, yi  r3 , yi 1  yi  1 r1  2r2  2r3  r4 . 6 Геометрический смысл метода Рунге-Кутта Из точки xi , yi  сдвигаются в направлении, определенном углом h r  1 tg1  f xi , yi  . На этом направлении выбирают точку  xi  , yi  1  . 2 2  Из точки xi , yi  сдвигаются в направлении, определенном углом h r h r     2 tg 2  f  xi  , yi  1  , и на нем выбирают точку  xi  , yi  2  . Из 2 2 2 2   точки сдвигаются в направлении под углом xi , yi  h r    3 tg 3  f  xi  , yi  2  и выбирают точку xi  h, yi  r3  , в которой 2 2  угол наклона касательной к интегральной кривой равен  4 tg 4  f xi  h, yi  r3  . Четыре направления усредняются, причем  2 и  3 берутся удвоенными. На этом усредненном направлении и выбирается точка приближенного решения xi1 , yi1   xi  h, yi  yi . 15.6. Решение систем обыкновенных дифференциальных уравнений первого порядка Изученные методы легко распространяются на решение задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка, разрешенных относительно производных. Рассмотрим случай двух уравнений:  y  f1 x, y, z    z  f 2 x, y, z   yx0   y0   z x0   z0 14 Будем считать, что в окрестности начальных данных, точки x0 , y0 , z0 , выполняется условие существования и единственности решения задачи Коши, в том числе: f1 , f 2 – непрерывно дифференцируемы, существуют их частные производные 1-го и 2-го порядков. Численное решение будем искать в виде двух таблично заданных функций y  yx  , z  zx  на множестве точек xi  x0  ih, i  1,2,, n . Для поиска значений функций можно использовать любой метод, в частности метод Эйлера:  y0  y  x0 , z0  z  x0   y  y  h  f x , y , z   i 1 i 1 i i i   zi 1  zi  h  f 2  xi , yi , zi   xi 1  xi  h i  0,1, n  1. 15.7. Численное решение обыкновенных дифференциальных уравнений высших порядков Основная идея методов численного решения обыкновенных дифференциальных уравнений n-го порядка (n>1) заключается в замене этого уравнения на систему n обыкновенных дифференциальных уравнений первого порядка. Для примера рассмотрим обыкновенное дифференциальное уравнение 3-го порядка, разрешенное относительно старшей производной y  f x, y, y, y . Требуется отыскать решение уравнения – функцию y  yx  , удовлетворяющую начальным условиям yx0   y0 , yx0   z0 , yx0   u0 . Произведем замену: yx   zx  , yx   ux  . Исходное уравнение эквивалентно системе с заданными начальными условиями  y  x0   y0  y  z    z  u  z  x0   z0 . u  f  x, y, z, u  u  x   u   0 Решением данной задачи Коши являются три функции y  yx  , z  zx  , u  ux  , среди которых первая и является искомым решением исходного дифференциального уравнения. Численное решение также ищется в табличном виде на множестве точек xi  x0  ih, i  1,2,, n . Для поиска значений функции можно использовать любой метод, в 15 частности метод Эйлера:  y0  y  x0 , z0  y x0 , u0  y x0  y  y  h  z i i  i 1  zi 1  zi  h  ui u  u  h  f  x , y , z , u  i i i i i  i 1  xi 1  xi  h i  0,1, n  1.