14 Численные методы решения систем обыкновенных дифференциальных уравнений. Явные методы. Неявные методы (Лекции по теории оптимизации и численным методам)

Лекция 148. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ОБЫКНОВЕННЫХДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙПОСТАНОВКА ЗАДАЧИРассматривается проблема решения систем обыкновенных дифференциальныхуравнений первого порядка, связывающих независимую переменную x , неизвестныефункции y1 ( x ),..., y n ( x ) и их производные y1 ( x ),..., y n ( x ) .В случае, если уравнения разрешимы относительно производных, систему можнозаписать в нормальной форме Коши:dy1 f1 ( x , y1 ,..., y n ) ,dxdy 2 f 2 ( x , y1 ,..., y n ) ,dx.................................dy n f n ( x , y1 ,..., y n ) ,dxгде f i ( x , y1 ,..., y n ) , i  1, n , – известные функции.Решением системы называется совокупность n функций y1 ( x ),..., y n ( x ) , непрерывных на некотором интервале (a, b ) , такая, что подстановка этих функций в системуобращает все уравнения в тождества.Задача Коши для системы состоит в нахождении решения системы, удовлетворяющего начальным условиям:y1 ( x 0 )  y10 , y 2 ( x 0 )  y 20 ,..., y n ( x 0 )  y n 0 ,где y10 , y 20 ,..., y n 0 – известные числа.В векторной форме задача Коши имеет видY   F ( x ,Y ),Y (x0 )  Y 0 ,где Y  ( y1 ,..., y n )T , F ( x ,Y )  ( f1 ( x ,Y ),..., f n ( x ,Y ))T , Y 0  ( y10 ,..., y n 0 )T .Теорема (о существовании и единственности решения задачи Коши).Пусть выполнены следующие условия:а) функции f i ( x , y1 ,..., y n ) , i  1, n , определены и непрерывны в некоторой замкнутой области D , а также имеют в D ограниченные частные производные по переменным y1 ,..., y n ;б) точка ( x 0 , y10 , y 20 ,..., y n0 ) лежит внутри области D .Тогда решение задачи Коши существует и единственно.134З а м е ч а н и я.1.

Во многих практических приложениях независимая переменная обозначаетсячерез t и имеет смысл времени, поэтому задача Коши называется начальной задачей.2. Чтобы решить задачу Коши для дифференциального уравнения n -го порядка:y (n )  f ( x , y ( x ),..., y (n 1) ( x )) ,y ( x 0 )  y 0 , y ( x 0 )  y 0 ,..., y (n 1) ( x 0 )  y 0(n 1) ,где x 0  (a, b ) , y 0 , y 0 ,..., y 0(n 1) – заданные числа, ее необходимо привести к системе nуравнений первого порядка. Обозначая y1 ( x )  y ( x ), y 2 ( x )  y ( x ),..., y n ( x )  y (n 1) ( x ) ,получаемdy1 y2 ,y1 ( x 0 )  y 0 ,dxdy 2 y3 ,y 2 ( x 0 )  y 0 ,dx..............................................................dy n f ( x , y1 ,..., y n ) , y n ( x 0 )  y 0(n 1) .dx3. Чтобы упростить изложение и в силу того, что численные методы легко обобщаются на системы уравнений, в дальнейшем будем рассматривать решение задачи Коши для уравнения первого порядкаy   f ( x , y ),y ( x 0 )  y 0 , x  (a, b ) .(*)Чтобы записать формулы для решения задачи Коши необходимо заменить функцию y (x ) на вектор-функцию Y (x ) , f ( x , y ) на F ( x ,Y ) , а y 0 – на Y 0 .ПРИНЦИПЫ ФОРМИРОВАНИЯ ЧИСЛЕННЫХ МЕТОДОВЧисленное решение задачи (*) ищется в узлах сетки  n  x 0 , x1 ,..., x n , гдеhi 1  x i 1  x i , i  0, n  1 , – расстояние между соседними узлами, называемое шагоминтегрирования (параметром сетки).

Если hi 1  h  const , сетка называется равномерной (регулярной), а если hi 1  var – неравномерной (нерегулярной). В случае равномер-ной сетки узлы находятся по формуле x i  x 0  ih, i  0, n .Решение находится в виде последовательности значений yˆ0 , yˆ1 , yˆ2 ,..., yˆn , являющихся приближением значений y 0 , y ( x1 ), y ( x 2 ),..., y ( x n ) точного решения y(x ) в узлахсетки  n (рис. 1).135yŷ iy ( xi 1 )y ( x1 )y ( xi 1 )ŷ i 1y ( xi )y  y(x )ŷ i 1ŷ1ŷ 0y ( xn )y00ŷ nx0x1xi 1xixi 1xnxРис. 1Численные дискретные методы решения обыкновенных дифференциальных уравнений, позволяющие найти решение только в узлах сетки, делятся на две группы: явные инеявные.Значение ŷi 1 на (i  1) -м шаге может определяться явно:yˆi 1  ( x i  k 1 ,..., x i 1 , x i , yˆi  k 1 ,..., yˆi 1 , yˆi ) ,где (.) – некоторая функция, зависящая от конкретного метода (кроме последней рассчитанной точки ( x i , yˆi ) могут использоваться еще (k  1) предыдущих точек), илинеявно:yˆi 1  ( x i  k 1 ,..., x i 1 , x i , x i 1 , yˆi  k 1 ,..., yˆi 1 , yˆi , yˆi 1 ) ,где искомая величина ŷi 1 входит одновременно и в левую, и в правую часть.Явные и неявные методы делятся также на одношаговые и многошаговые ( k шаговые).

В одношаговых методах для расчета очередной точки ( x i 1 , yˆi 1 ) требуетсяинформация только о последней рассчитанной точке ( x i , yˆi ) . В k -шаговых методах длянахождения точки ( x i 1 , yˆi 1 ) требуется информация о k предыдущих точках.Формулы явных или неявных методов в общем случае представляют собой нелинейные уравнения относительно ŷi 1 и называются разностными схемами.Локальной ошибкой численного метода на (i  1) -м шаге называется величина i 1 (h)  yˆi 1  y ( x i 1 ) ,где y ( x i 1 ) – значение точного решения при x  x i 1 , а ŷi 1 – приближенное решение,получаемое по формулам при условии, что вместо приближенных значенийyˆi , yˆi 1 ,..., yˆi k 1 используются значения, соответствующие точному решению, т.е.y ( x i ), y ( x i 1 ),..., y ( x i k 1 ) .136Глобальной ошибкой называется величина e n (h)  yˆn  y ( x n ) , где ŷ n – значение,получаемое по формулам при i  n  1 .Глобальная ошибка определяется:а) ошибками округления и ошибками арифметических действий, обусловленнымичислом разрядов компьютера и характером выполняемых операций для расчета значенияискомой функции в очередной точке x i 1 ;б) методическими ошибками, определяемыми выбранным алгоритмом;в) переходными ошибками, обусловленными тем, что при расчете значения ŷi 1вместо точных значений y ( x i ), y ( x i 1 ),..., y ( x i  k 1 ) берутся приближенные значенияyˆi , yˆi 1 ,..., yˆi  k 1 , полученные на предыдущих шагах.Локальные ошибки «переносятся» в точку x n и формируют глобальную ошибку.Число p называется порядком (точностью) численного метода, если его глобальная ошибка есть О большое от h p , т.е.

e n (h)  O (h p ) .Пояснение. Пусть R (h) – некоторая функция переменной h (как правило, R (h) –остаточное слагаемое некоторой аппроксимационной формулы) с конечной областью определения DR на полуоси h  0 , причем h  DR . Тогда, если при некотором h  h0справедливо неравенство R (h)  ch k , где c  const , не зависящая от h , k – целое число, h0  0 , то пишут R (h)  O (h k ) и говорят, что R (h) есть «O большое от h k » приh  0.На практике в качестве характеристики точности метода часто используется велиyˆi  y ( x i ) .чина (h)  maxi  0,1,..., nМожно показать, что если локальная ошибка имеет порядок ( p  1) , т.е. i 1 (h)  O (h p 1 ) , то глобальная погрешность имеет на единицу меньший порядок, т.е.e n (h)  O (h p ) .Перейдем теперь к рассмотрению устойчивости численных методов.

Она проверяется на «тестовом примере»y    y,y (0)  1 ,где  – в общем случае комплексная константа. Дифференциальное уравнение являетсяпростейшим линейным уравнением, и для него можно получить значимые критерии устойчивости в явной форме.Метод называется ограниченно устойчивым, если существует такое число hкр.  0 ,что при использовании метода для решения тестового примера, где Re   0 , с шагом0  h  hкр. при i   глобальная ошибка ограничена. Величина hкр. называется критическим шагом.

Если h  hкр. , глобальная ошибка может неограниченно возрастать. В ограниченно устойчивых методах при задании величины шага h необходимо учитыватьзначение критического шага hкр. . Для сложных дифференциальных уравнений и систем137нахождение hкр. является самостоятельной задачей, а свойство ограниченной устойчивости предупреждает вычислителя о возможных проблемах. Поэтому на практике становится актуальной задача конструирования таких методов, которые были бы устойчивыпри любом значении шага, а его величина выбиралась бы только исходя из желаемойточности расчетов (при этом класс решаемых задач может быть ограничен).Метод называется А-устойчивым, если при его применении с любым фиксированным положительным шагом h все численные решения тестового примера с комплекснойконстантой  ( Re   0 ) стремятся к нулю при i   .А.

ЯВНЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ОБЫКНОВЕННЫХДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙА1. Явный метод ЭйлераРассмотрим проблему нахождения численного решения задачи Коши:dy f ( x , y ),dxy(x 0 )  y0 .Вводится в общем случае неравномерная сетка n  (x0 , x1,...,xi , xi 1,...,xn ) . Величина шага hi 1  x i 1  x i выражается через узловые точки. Для аппроксимации производной dy   dx x  xiиспользуем dy Ш 2,i  x i , x i 1  :   dx x  xiформулу,yi 1  yihi 1записанную O (hi 1 )надвухточечномшаблоне h i 1M 2,i  . 2Далее заменяется правая часть уравнения ее сеточным представлением, т.е.f ( x , y )  f ( x i , yi ) , а вместо y ( x ) рассматривается сеточная функция yˆi  y ( x i ) , которая определяется только в точках сетки.

Выполняется подстановка аппроксимаций производной и правой части в дифференциальное уравнение:y i 1  y ihi 1 O(hi 1 ) f ( xi , y i ) .После отбрасывания остаточных слагаемых получается явная схема Эйлера первого порядка (явный метод Эйлера):yˆi 1  yˆi  hi 1 f ( x i , yˆi ) ,i  0, n  1 , yˆ0  y 0 ;Порядок точности метода, как правило, определяется порядком аппроксимациисхем, явный метод Эйлера является ограниченно устойчивым с критическим шагом2hкр.

  (см. тестовый пример).138А2. Метод Эйлера-КошиДля аппроксимации производной применяется формула: dy   dx x  xi h2M 3,i  . 6y i  1  y i 1 O (h 2 )2hВыполняется подстановка аппроксимаций производной и правой части в дифференциальное уравнение:y i 1  y i 12h O(h 2 )  f ( x i , y i ) .После отбрасывания остаточных слагаемых получается явная схема метода Эйлера–Коши второго порядка:yˆi 1  yˆi 1  2h  f ( x i , yˆi ) , i  1, n  1 .Для начала расчетов требуется иметь две «разгонные» точки yˆ0 , yˆ1 . Первая определяется известным начальным условием yˆ0  y 0 , а вторая может быть найдена с помощью другого метода, например, по формуле: yˆ1  y 0  h1 f ( x 0 , y 0 ) .А3.

Модифицированный метод ЭйлераМодифицированный метод Эйлера второго порядка:yˆi12 yˆi hi 12f ( x i , yˆi ) ,i  0, n  1 ,hyˆi 1  yˆi  hi 1 f  x i  i 1 , yˆi  1  ,22 i  0, n  1 .Интервал устойчивости h   (2, 0) (здесь  – действительное число в тестовомпримере) модифицированного метода Эйлера совпадает с интервалом устойчивости явного метода Эйлера.А4.