Обыкновенные дифференциальные уравнения

ЛЕКЦИИ №9,10. ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ. ЗАДАЧА КОШИ 1. Общие замечания Обыкновенными дифференциальными уравнениями (ОДУ) называются уравнения с одной независимой переменной. Если независимых переменных больше, чем одна, то уравнение называется дифференциальным уравнением с частными производными. С помощью обыкновенных дифференциальных уравнений строятся модели движения систем взаимодействующих частиц, электротехнических процессов в электрических цепях, кинетики химических реакций, процессов заселения уровней энергии в высокотемпературных средах и многих других объектов и процессов. К задачам для обыкновенных дифференциальных уравнений сводятся некоторые задачи для уравнений в частных производных, когда многомерное уравнение позволяет провести разделение переменных (например, при вычислении энергетического спектра частиц в полях определенной симметрии). Обыкновенное дифференциальное уравнение любого порядка при помощи замены переменных может быть сведено к системе уравнений первого порядка. Рассмотрим в связи с последним пример. Дифференциальное уравнение третьего порядка a x  d 3v d 2v dv      b x  c x  d  x v  f  x  dx 3 dx 2 dx заменой переменных d 2v dv  v2 ,  v1 , 2 dx dx приводится к следующей системе дифференциальных уравнений dv  v1 , dx dv 1  v2 , dx dv a x  2  b x v 2  c x v 1  d  x v  f  x . dx В общем виде преобразование выглядит следующим образом: дифференциальное уравнение n - го порядка v ( n ) ( x )   x , v , v , v ,..., v  n 1 , заменой переменных v ( k )  vk сводятся к системе n уравнений первого порядка vk  vk 1 , 0  k  n  2, vn 1 ( x )   x ,v0 ,v1 ,v2 ,...,vn 1 , где обозначено v0  v. В соответствии с изложенным выше далее будут рассматриваться системы уравнений первого порядка: vk ( x )   k x ,v1 ,v2 ,...,vn , 1  k  n. Решение системы n -го порядка зависит от n параметров c1 ,c2 ,...,cn . .Для выделения единственного решения необходимо использование дополнительных условий для искомой функции. В зависимости от того, каким образом ставятся данные условия, различают три типа задач для обыкновенных дифференциальных уравнений: задача Коши, краевая задача и задача на собственные значения. В задаче Коши все дополнительные условия ставятся в одной точке: vk ( x0 )  vk ,0 , 1 k  n. Решение отыскивается в некотором интервале x0  x  xl . Если правые части точки x  k уравнений непрерывны в некоторой окрестности начальной ,v1,0 ,v 2 ,0 ,...,v n ,0  и удовлетворяют условию Липшица по переменным v k , то ре- шение задачи Коши существует, единственно и непрерывно зависит от координат начальной точки, т.е. задача является корректной. Условие Липшица формулируется следующим образом  k x , v1 ,l ,v 2 ,l ,...,v n ,l    k x ,v1 ,m ,v 2 ,m ,...,v n ,m    Lv1 ,l  v1 ,m  v 2 ,l  v 2 ,m ... v n ,l  v n ,m , для любых точек x , v 1 ,l ,v2 ,l ,...,vn ,l , x ,v1,m ,v2 ,m ,...,vn ,m  . 2. Методы решения Можно выделить три метода решения обыкновенных дифференциальных уравнений: точные, аналитические приближенные и численные. Точные методы предусматривают получение решения в виде комбинации элементарных функций или в виде квадратур от последних. Возможности точных методов ограничены. Приближенные методы сводятся к построению последовательности функций wn ( x ) , имеющих пределом искомую функцию v( x ) . Обрывая эту последовательность на каком-то номере k , получают приближенное решение. Наиболее универсальными методами решения являются численные. Их основной недостаток - возможность получения только частного решения. Следует иметь в виду следующее обстоятельство. Успех от применения численного метода сильно зависит от обусловленности задачи, т.е. задача должна быть хорошо обусловлена, а именно, малые изменения начальных условий должны приводить к малому изменению решения. В противном случае (слабой устойчивости) малые погрешности в начальных данных или погрешности численного метода могут приводить к большим погрешностям в решении. Далее будут рассматриваться алгоритмы решения задачи Коши на примере одного уравнения первого порядка v ( x )  ( x , v ) . Обобщение на случай системы n уравнений осуществляется заменой v( x ) на v( x ) и  ( x , v ) на v1 v( x )  v2 .  ( x,v )  , 1 2   ( x , v ) , где . n vn 2.1. Метод Пикара Данный метод является представителем приближенных методов решения рассматриваемого класса задач. Идея метода чрезвычайно проста и сводится к процедуре последовательных приближений для решения интегрального уравнения, к которому приводится исходное дифференциальное уравнение. Пусть поставлена задача Коши v x    x , v( x ) , (1) x0  x  x l . , vx0   v0 . Проинтегрируем выписанное уравнение x vx   v0   ( t , v( t )) dt . (2) x0 Процедура последовательных приближений метода Пикара реализуется согласно следующей схеме x y s x   v0   ( t , y s 1 ( t )) dt , x0 причем y0 ( t )  v0 , ( i – номер итерации). Пример. Решить методом Пикара уравнение (3) v ( x )  x 3  v 3 , v( 0 )  0 . Решение этого уравнения не выражается через элементарные функции: x4 , y 1 ( x )  0   t dt  4 x 3 x y2 ( x )  0   [ t 3  ( t4 3 x4 1 ) ] dt  ( 1  2 x9 ) , 4 4 4  13 t4 1 y 3 ( x )  0   { t [ ( 1  2 t 9 )] 3 }dt  4 4  13 x 3 x4 x 13 3 3 1   3 (1  2 x9  4 x 18  6 x 27 ) 2 4 4  13 4  22 4  13  31 4  13  40 и т.д. Видно, что при x  1 ряд быстро сходится. Метод удобен, если интегралы можно взять аналитически. Можно доказать, что метода Пикара сходится, если в некоторой ограниченной области g( x , v ) правая часть  ( x , v ) непрерывна и, кроме того, удовлетворяет условию Липшица по переменной v т.е.  ( x , v1 )   ( x , v 2 )  L v1  v 2 , где L - некоторая константа. 2.2. Методы Рунге-Кутта Данные методы являются численными. На практике применяются методы РунгеКутта, обеспечивающие построение разностных схем (методов) различного порядка точности. Наиболее употребительны схемы (методы) второго и четвертого порядков. Их мы и рассмотрим ниже. Предварительно введем некоторые понятия и определения. Сеткой на отрезке [ a ,b ] называется фиксированное множество точек этого отрезка  N . Функция, определен- ная в данных точках, называется сеточной функцией. Координаты точек x i удовлетворяют условиям a  x0  x1  x 2 ... x N 2  x N 1  x N  b . Точки xi  N являются узлами сетки. Равномерной сеткой на [ a ,b ] называется множество точек  h  { xi  a  ih } , где h  i  0,1,2,..., N , ba - шаг сетки. N При решении дифференциальных уравнений приближенным методом основным является вопрос о сходимости. Применительно к разностным методам традиционно более употребительно понятие сходимости при h  0 . Обозначим значения сеточной функции y i , значения точного решения дифференциального уравнения (5.1) в узле i - v( xi ) ( y i являются приближенными значениями v( xi ) ). Сходимость при ющее. Фиксируем точку x и строим совокупность сеток h  0 означает следу-  h таким образом, что h  0 и xi  a  ih  x (при этом i   ). Тогда считают, что численный метод сходится в точке x , если y i  v( x i )  0 при h  0 , xi  x . Метод сходится на отрезке [ a ,b ] , если он сходится в каждой точке x  [ a , b ] . Говорят, что метод имеет p -й порядок точности, если можно найти такое число p  0 , что y i  v( x i )  O( h ) при h  0 . p Введем далее понятие невязки или погрешности аппроксимции разностного уравнения, заменяющего заданное дифференциальное уравнение, на решении исходного уравнения, т.е. невязка  i представляет собой результат подстановки точного решения уравнения (1) v( x ) в разностное уравнение. Например, (1) можно заменить следующим простейшим разностным уравнением y i 1  y i  ( xi , y i )  0 , h i  0 ,1,2,...,y0  v0 . Тогда невязка определится следующим выражением i   u i 1  u i   ( xi ,u i ) . h Приближенное решение не совпадает вообще говоря с u i , поэтому невязка  i в i ой точке не равна нулю. Вводят следующее определение: численный метод аппроксимирует исходное дифференциальное уравнение, если невязка  i  0 при h  0 , и имеет p -й порядок аппроксимации, если  i  O( h p ) . Доказывается, что порядок точности численного метода решения дифференциального уравнения совпадает с порядком аппроксимации при достаточно общих предположениях. Теперь перейдем к анализу схем Рунге-Кутта. Сначала обратимся к схемам второго порядка точности. Используя формулу Тейлора, решение дифференциального уравнения (1) можно представить в виде 1 vn 1  vn  hn vn  hn2 vn ... , 2 где обозначено (6) v n  v( x n ) , vn  v( xn ) , hn  xn 1  xn . Согласно (1) vn  ( xn , vn ) , vn   x ( xn ,vn )   v ( xn ,vn )( xn ,vn ) . Далее удерживаем только выписанные члены ряда. Представим вторую производную следующим образом v n  ( v n )   ( ~x , v~ )   ( x n , v n ) x x ,~ v - пока неизвестные величины. Пусть где ~ ~ x  xn   h, v~  vn   h . , Обозначим приближенное значение решения в узле с номером n через y n (именно это решение будет получаться после того, как мы ограничим ряд членами с порядком не выше второго). Имеем y n  1  y n  hn  ( x n , y n )  1 2  ( x n   hn , y n   hn )   ( x n , y n ) hn [ ] 2 x  y n  hn [  ( x n , y n )   ( x n   hn , y n   hn )] . Введенные здесь параметры  , , и  подлежат определению. Разлагая правую часть в ряд Тейлора до линейных членов и приводя подобные члены, получим последовательно y n 1  y n  hn   ( x n , y n )    [ ( x n , y n )    x ( x n , y n ) hn   y ( x n , y n ) hn ]    y n  (    )hn( x n , y n )   hn2 [   x ( x n , y n )    y ( x n , y n )] . Условием выбора параметров    1 ,  , ,   и 1 , 2  поставим близость выражения (7) ряду (6), тогда 1 2   ( x n , y n ) . Один параметр остается свободным. Пусть это будет  1 ,   (7)  , тогда 1 1 ( x n , y n ) ,  2 2 и окончательно из (7) с учетом найденных отношений для y n 1  y n  hn ( 1   )( x n , y n )   ( x n   , и  получим h 1 hn , y n  n  ( x n , y n )]  2 2 (8) Соотношение (8) описывает однопараметрическое семейство двучленных формул Рунге-Кутта. В специальной литературе доказывается, что если  ( x , v ) непрерывна и ограниче- на вместе со своими вторыми производными, то приближенное решение схемы (8) равно2 мерно сходится к точному решению с погрешностью O(max hn ) , т.е. схема (.8) обладает вторым порядком точности. В практике расчетов используют формулы (8) при значениях параметра   1 .Рассмотрим эти варианты. Случай  1 . 2 Из (5.8) выводим y n 1  y n  h  xn , y n    xn  hn , y n  hn xn , y n , 2 (9) Применение формулы (9) сводится к следующей последовательности шагов: 1. Вычисляется грубо значение функции y n 1 (по схеме ломаных ) y n 1  y n  hn xn , y n . 2. Определяется наклон интегральной кривой в точке ( x n 1 , y n 1 ) y n1   xn 1 , y n 1 . 3. Находится среднее значение производной функции на шаге hn y 1 n 2  1  xn , y n   y n1 , 2 4. Рассчитывается, наконец, значение функции в ( n  1 )-м узле  1 , 2 y n 1  y n  hy  1 . n 2 Данная схема имеет специальное название "предиктор - корректор". Случай   1. Согласно (8) получаем y n  1  y n  hn  [ x n  hn h , y n  n  ( x n , y n )] . 2 2 Задача решается посредством следующих шагов: 1. Вычисляется значение функции в половинном узле y n 1 2  yn  hn ( x n , y n ) . 2 2. Определяется значение производной в узле n  y 1 n 2  ( x n  1 2 hn ,y 1 ) . 2 n 2 3. Находится значение функции в ( n  1 )-м узле y n 1  y n  hn y  1 . n 2 Помимо рассмотренных выше двучленных схем широкое распространение в практике расчетов имеют схемы Рунге-Кутта четвертого порядка точности. Ниже даются без вывода соответствующие формулы y n 1  y n  k 1  2k 2  2k 3  k 4  / 6 , k 1  hn  x n , y n , k 2  hn  x n  hn / 2 , y n  k 1 / 2 , k 3  hn  x n  hn / 2 , y n  k 2 / 2 , k 4  hn  x n  hn , y n  k 3 . (10) Схемы с большим числом членов практически не применяются. Пятичленные формулы обеспечивают четвертый порядок точности, шестичленные формулы имеют шестой порядок, но их вид весьма сложен. Погрешности приведенных схем Рунге-Кутта определяются максимальными значениями соответствующих производных. Оценку погрешностей легко получить для частного случая вида правой части дифференциального уравнения ( x , v )  ( x ) . В этом варианте решение уравнения может быть сведено к квадратуре и все схемы разностного решения переходят в формулы численного интегрирования. Например, схема (9) принимает вид yn1  yn  hn  xn    xn1  , 2 то есть имеет форму метода трапеций, а схема (10) переходит в схему y n 1  y n  hn  x n   4 x n  hn / 2   x n  hn  , 2 представляющую собой формулу Симпсона с шагом hn . 2 Оценки погрешности формул трапеций и Симпсона известны. Например, мажорантные оценки погрешности указанных формул даются следующими выражениями x l  x0 max ( hn2 ) max (  '' ), 12 x  x0  l max ( hn4 ) max (  IV ) 2880 RТРАП  RСИМП Из приведенных формул видно, что точность схем Рунге-Кутта достаточно высока. Выбор той или иной из приведенных схем для решения конкретной задачи определяется следующими соображениями. Если функция  ( x , v ) в правой части уравнения не- прерывна и ограничена, а также непрерывны и ограничены ее четвертые производные, то наилучший результат достигается при использовании схемы (10). В том случае, когда функция  ( x , v ) не имеет названных выше производных, предельный (четвертый) порядок схемы (10) не может быть достигнут, и целесообразным оказывается применение более простых схем . Помимо схем Рунге-Кутта практический интерес представляют многошаговые методы, которые можно описать следующей системой уравнений a0 y n  a1 y n 1 ... a m y n  m  b0  n  b1 n 1 ... bm n  m , h (11) где n  m, m  1,..., а ak ,bk - числовые коэффициенты, k  0 ,1,2,..., m , a0  0 . Согласно данному уравнению расчет начинается со значения n  m . В этом случае получается соотношение вида a0 y m  a1 y m 1 ... a m y 0  b0 m  b1 m 1 ... bm 0 , h т.е. для начала счета надо иметь m начальных значений y i , i  0 ,1,2,..., m  1. Эти значения y i приходится вычислять каким-либо другим методом, например, методом Рунге-Кутта. В необходимости использовать разные методы счета состоит неудобство многошаговых методов. Среди многошаговых методов наиболее распространен метод Адамса, схема реализации которого следует из (11) при a0  a1  1 и a k  0 для k  2,3,..., m : y n  y n 1 m   bk  n  k . h k 0 При b0  0 метод Адамса оказывается явным, а при b0  0 - неявным. 2.3. Неявные методы Рассмотрим неявную схему Эйлера для уравнения (1).. Она имеет первый порядок точности y n 1  y n   ( xn 1 , y n1 ) . h (12) Для определения неизвестного значения y n1 придется решать нелинейное уравнение. Необходимость решать уравнение – типичная ситуация для неявных методов. И в этом состоит сложность и трудоемкость их применения, учитывая, что в большинстве случаев уравнения получаются нелинейными, и возникает проблема отыскания корня, который может быть не единственный, а может и вообще не существовать при выбранном шаге. Однако неявные методы обладают свойством устойчивости, и несмотря на все трудности их применяют очень широко, т.к. они позволяют вести расчет с увеличенными шагами по сравнению с явными методами. Среди других неявных методов можно назвать метод трапеций. При решении жестких систем дифференциальных уравнений хорошо зарекомендовал себя метод Гира, который относится к чисто неявным многошаговым разностным методам, общая формула которых выглядит следующим образом: m a k 0 k y n  k  h( x n , y n ) , При m  1 и a0  1, a1  1 имеем y n  y n 1  h( x n , y n ) , т.е. неявный метод Эйлера. При m  2 и m  3 методы выглядят следующим образом 3 1 y n  2 y n 1  y n  2  h( x n , y n ) , 2 2 (13) 11 3 1 y n  3 y n 1  y n  2  y n  3  h( x n , y n ) . 6 2 3 (14) Разностное уравнение (13) имеет второй порядок точности, а (14) - третий.