Численные методы решения задачи Коши для систем обыкновенных дифференциальных уравнений.

Лекция 6. Стр. 1 Лекция 6 Численные методы решения задачи Коши для систем обыкновенных дифференциальных уравнений Метод Эйлера, методы Рунге-Кутты. Многошаговые методы Адамса Основные понятия и обозначения. Обыкновенным дифференциальным уравнением первого порядка называется уравнение вида F x, y( x), y( x)  0 , где F — известная функция трех переменных, определенная в области G  R 3 , x — независимая переменная, x  a, b, y (x ) — неизвестная функция, y (x ) — ее производная. Функция y  y (x) называется решением дифференциального уравнения, если она непрерывно дифференцируема на a, b  и F x, y( x), y( x)  0 при всех x  a, b . График решения дифференциального уравнения называют интегральной кривой дифференциального уравнения. В дальнейшем будем рассматривать обыкновенные дифференциальные уравнения в нормальной форме, т.е. уравнения, разрешенные относительно производной.: y   f ( x, y ( x )) .(1) Если дифференциальное уравнение первого порядка y   f ( x, y ) , x, y   G имеет в области G решение, то, вообще говоря, таких решений бесконечно много, они могут быть заданы в виде y  y ( x, C ) , где C — произвольная константа; x, y( x, C)  G и y( x, C)  f x, y( x, C) . Однако, если требуется найти решение, удовлетворяющее начальному условию y ( x0 )  y0 , то при определенных условиях такая задача имеет единственное решение. Задача об отыскании решения дифференциального уравнения, удовлетворяющего заданному начальному условию, называется задачей Коши. Начальные условия для обыкновенных дифференциальных уравнений называют условиями Коши Справедлива следующая теорема существования и единственности Теорема. Если функция f ( x, y ) и ее частная производная f ( x , y ) непрерывны в области y G , x0 , y0   G , то на некотором интервале x0  h, x0  h существует единственное решение уравнения y   f ( x, y ) , удовлетворяющее начальному условию y ( x0 )  y0 . Для дифференциального уравнения y   f ( x, y ) теорема существования и единственности имеет простую геометрическую интерпретацию: через каждую точку x0 , y0   G проходит только одна интегральная кривая y  y( x, C0 ) семейства y  y ( x, C ) такая, что y ( x0 , C 0 )  y 0 . Обыкновенным дифференциальным уравнением n -го порядка называется уравнение вида F x, y( x), y ( x), y ( x), ..., y ( n) ( x)  0 , где F — известная функция n  2 переменных,   определенная в области D  R , x — независимая переменная x  a, b, y (x ) — неизвестная функция, n — порядок уравнения. Функция y  y (x) называется решением n 2 дифференциального уравнения, если она n раз непрерывно дифференцируема на a, b  и   F x, y, y , y , ..., y ( n )  0 для всех x  a, b. В дальнейшем будем рассматривать обыкновенные дифференциальные уравнения в нормальной форме: y ( n)  f ( x, y, y , ..., y ( n1) ) . Дифференциальное уравнение n -го порядка имеет, вообще говоря, бесконечное множество решений. Однако задача об отыскании решения, удовлетворяющего начальному условию y( x0 )  y0 , y( x0 )  y01 , ..., y ( n 1) ( x0 )  y 0,n 1 , при определенных ограничениях на правую часть уравнения имеет единственное решение. Справедлива следующая теорема. Лекция 6. Стр. 2 Теорема. Если функция f ( x, y) , y   y1 , y 2 , ..., y n  , и ее частные производные f ( x, y) , i  1, 2, ..., n непрерывны в области D , x, y    x, y1 , y 2 , ..., y n   D  R n1 , то на yi некотором интервале x0  h, x0  h существует единственное решение уравнения y ( n)  f ( x, y, y , ..., y ( n1) ) , удовлетворяющее начальным условиям: y( x0 )  y0 , y( x0 )  y01 , ..., y ( n 1) ( x0 )  y 0,n 1 , x0 , y 0   x0 , y0 , y01 , ..., y0,n1  D . Задача Коши для дифференциального уравнения n -го порядка y ( n)  f ( x, y, y , ..., y ( n1) ) , y( x0 )  y0 , y( x0 )  y01 , ..., y ( n 1) ( x0 )  y 0,n 1 может быть сведена к задаче Коши для системы n -го порядка ( системы n дифференциальных уравнений 1-го порядка), которая в векторной форме имеет вид    Y ) , Y ( x0 )  Y 0 , где Y   y1 ( x), y 2 ( x), ..., y n ( x)  ,   F ( x, Y )   y2 , y3 , ..., yn , f ( x, y1 ,..., yn  .  Y  F ( x, Обусловленность задачи Коши. В силу приведенных выше теорем существования и единственности решения задачи Коши и известной из курса дифференциальных уравнений непрерывной зависимости решения задачи Коши от начальных данных, эта задача корректна. Однако применять численные методы решения задачи имеет смысл только к хорошо обусловленным задачам. В инженерных и научных расчетах достаточно часто возникают плохо обусловленные задачи. Пример плохо обусловленной задачи Коши. Рассмотрим задачу Коши y '  y  x , y (0)  1 , x  0, 100 . Легко видеть, что общее решение задачи y ( x, c)  1  x  c  exp( x) , а решение задачи Коши y ( x)  1  x и y (100 )  101 . Однако, если начальные данные содержат малую погрешность, y (0)  1.000001 , то y( x, c)  1  x  10 6  exp( x) и y(100 )  2.7  10 37 . Здесь важно понимать, что плохая обусловленность задачи проявляется на относительно длинном промежутке интегрирования уравнения. Модельная задача. Полезно рассмотреть модельную задачу y'  y, y(0)  y0   0 , решение которой y  y (t ) имеет вид y   y0   0   expt  . Погрешность , вызванная возмущением начальных условий  0 , определяется формулой  0 expt  . Очевидно, что задача хорошо обусловлена при Re   0 и плохо обусловлена при Re   0 . В этой задаче видно как распространяется возмущение начальных данных. Задача Коши на конечном отрезке. Рассмотрим процесс распространения погрешностей ,внесенных в начальные значения. Пусть y0 * - возмущенное начальное значение,  0  y0  y0 * - его погрешность, а y * ( x) - решение соответствующей задачи Коши ( y*)( x)  f ( x, y * ( x)), (2)   y * ( x0 )  y 0 * Вычтем из уравнения (1) уравнение (2) и воспользуемся формулой конечных приращений Лагранжа:  f ( x, y( x))  f ( x, y * ( x))  f y ( x, ~ y ( x))( y( x)  y * ( x)), где ~y ( x) некоторое промежуточное значение. В результате получим, что погрешность  ( x)  y ( x)  y * ( x) удовлетворяет дифференциальному уравнению:  ( x)   ( x) ( x),  ( x)  f y  ( x, ~y ( x)) и начальному условию  ( x0 )   0 . Решение этой задачи Коши выражается формулой Лекция 6. Стр. 3  x   x0   ( x)   0 exp    ( )d  .Таким образом,  x   x  x0   x0     y ( )) d  играет в задаче Коши роль величина C ( x )  exp   ( ) d   exp  f y ( , ~   коэффициента роста ошибки. Заметим, что знак производной f y оказывает существенное  влияние на поведение погрешности  (x ) . Если f y  0 , то величина C (x) , а вместе с ней и модуль погрешности монотонно возрастают. При этом соответствующие интегральные кривые  расходятся. ( см. рис.) Иначе ведет себя погрешность в случае f y  0 . Здесь C (x) и  (x ) с ростом x монотонно убывают, а соответствующие интегральные кривые сближаются. Ошибка, внесенная в начальное значение, имеет тенденцию к затуханию. В случае, когда  производная f y незнакопостоянна, поведение погрешности может быть более сложным. Постановка задачи численного решения задачи Коши. Аналитическое выражение для решений дифференциальных уравнений, за исключением линейных дифференциальных уравнений с постоянными коэффициентами, удается получить достаточно редко. Численное решение задачи Коши состоит в построении таблицы приближенных значений y1 , y2 , ..., y N решения y (x ) в узлах сетки x1 , x2 , ..., x N , y( xi )  yi . Если xk  x0  kh , k  1, 2, ..., N , то сетка называется равномерной, h — шаг сетки. Как уже говорилось, задача Коши для любого дифференциального уравнения n -го порядка сводится к задаче Коши для системы n дифференциальных уравнений 1-го порядка. Численное решение задачи Коши для системы дифференциальных уравнений состоит в построении таблицы приближенных значений yi ,1 , yi , 2 , ..., yi , N компонент yi (x) вектора решения y(x) в точках x1 , x2 , ..., x N . В дальнейшем будем заниматься численным решением задачи Коши для дифференциального уравнения первого порядка, записанного в нормальной форме. Используемые при этом подходы без труда переносятся на системы дифференциальных уравнений. Классификация методов численного решения задачи Коши. Численный метод решения задачи Коши называется одношаговым, если для вычисления решения в точке x0  h используется информация о решении только в точке x 0 . В противном случае метод называется многошаговым. Лекция 6. Стр. 4 Если при этом yi   yi  k , yi  k 1 , ..., yi 1 , h , т.е. решение в точке xi явно выражается через значения в "предыдущих" точках, xi  k , xi  k 1 , ..., xi 1 то метод относится к явным методам. Если же yi   yi  k , yi  k 1 , ..., yi 1 , yi , h  , то метод неявный. Для оценки погрешности метода на одном шаге сетки обычно точное решение раскладывают по формуле Тейлора в окрестности узла xi : 1 y ( xi )h 2  ...  O (h p )  y i  f ( xi , y i )h  ...  O (h p ) . 2! Если значение y в точке xi известно, то в силу равенства y   f ( x, y ) значение производной y (x ) также можно считать известным. Для того чтобы вычислить производные y ( x), y ( x),... более высокого порядка, входящие в формулу, продифференцируем равенство y   f ( x, y ) по y ( xi 1 )  y ( xi  h)  y ( xi )  y ( xi )h  t , используя правило дифференцирования сложной функции. Получим y   f t  f y y   f t  f y f и т.д. Если расчетные формулы численного метода согласуются с разложением по формуле Тейлора до членов порядка h p , то число p называется порядком метода. Построение численных методов решения задачи Коши. Основной этап численного решения задачи Коши состоит в построении дискретной задачи Коши. Прежде чем описывать общий подход, построим простейшую дискретную задачу. y ( x  h)  y ( x ) . Обозначим yi  y ( xi ) и запишем уравнение y '  f ( x, y ) во h всех узлах сетки xi , заменяя точное значение производной y ' ( xi ) значением разностной Пусть y '  yi 1  yi . Подставив разностную производную в уравнение, получим h y  yi систему N линейных уравнений i 1  f ( xi , yi ) , i  1, 2, ..., N относительно N h неизвестных (значение y 0 известно из начальных условий). Эта система уравнений и представляет собой один из вариантов дискретной задачи Коши. Положив yi  y ( xi ) , i  1, 2, ..., N , получим формулу Эйлера численного решения задачи Коши приближенное решение yi 1  yi  hf ( xi , yi ) , вычисленное методом Эйлера. Итак, формулы yi 1  yi  hf ( xi , yi ) — расчетные формулы метода Эйлера приближенного решения задачи Коши y'  f ( x, y), y( x0 )  y0 . производной y ' ( xi )  Метод Эйлера допускает простую геометрическую интерпретацию. Пусть известна точка xi , yi  интегральной кривой уравнения y '  f ( x, y ) . Касательная к интегральной кривой уравнения, проходящая через эту точку определяется уравнением y  yi  f ( xi , yi )x  xi  . Следовательно, вычисленная методом Эйлера точка xi 1 , yi 1  , xi 1  xi  h , yi 1  yi  f ( xi , yi )h , лежит на этой касательной. Таким образом, после выполнения N шагов неизвестная интегральная кривая заменяется ломаной линией( ломаной Эйлера), для которой угловой коэффициент k n очередного n  го звена равен значению f ( xn , y n ) . Метод Эйлера представляет явный одношаговый метод. Для него погрешность аппроксимации имеет вид 1 y ( n ) h 2 . 2! Лекция 6. Стр. 5 Оценка погрешности. Для метода Эйлера справедлива такая оценка глобальной погрешности: max y( xn )  y n  C ( X ) 0 n  N M2 h, M 2  max y ( x) . [ x0 , X ] 2 Влияние вычислительной погрешности. Оценивая метод Эйлера необходимо учитывать, что при его реализации на ЭВМ неизбежно возникнут ошибки округления. В результате фактически вычисляемые значения y i * будут удовлетворять соотношению y *i 1  y *i  f ( xi , y *i )h   i . Величины  i учитывают вклад погрешностей округления. Оценка погрешности фактически вычисляемых значений y i * : M2  , M 2  [max h y ( x) . Оказывается, что полная 0 n  N x0 , X ] 2 h погрешность убывает только лишь при уменьшении шага h до некоторого значения hопт . max y( xn )  y *n  C ( X ) Вернемся к построению дискретной задачи Коши. Построение дискретной задачи заключается в следующем. Сначала строим на x0 , X  сетку x0 , x1 , ..., x N . В каждом узле сетки запишем дискретный аналог уравнения, заменяя производную y ' какой-либо разностной производной 1 k   j yi1 j , а правую часть заменяем какой-либо ее аппроксимацией — h j 0 f xi , yi   xi , yi 1k , yi 2k , ..., yi , yi 1 , h , построенной на тех же узлах. — y ' ( xi )  Полученная система N уравнений относительно N неизвестных значений 1 k  j yi1 j  xi , yi1k , yi2k , ..., yi , yi1 , h и есть дискретная задача Коши. h j 0 В частности, дискретная задача Коши для явного одношагового метода имеет вид  1 1 h  j yi 1 j  xi , yi , h  , или, что то же самое yi 1   1  yi  xi , yi , h .  h j 0 0 0 Приближенное решение задачи Коши в узлах сетки полагаем равным вычисленным значениям y i : yxi   yi . Основной характеристикой дискретной задачи Коши является степень аппроксимации. Например, об одношаговой явной дискретной задаче Коши говорят, что дискретное уравнение 1 1   j yi1 j  xi , yi , h аппроксимирует дифференциальное уравнение h j 0 y '  f ( x, y ) с порядком p, если max i зависит от h. 1 1   j y( xi 1 j )  xi , y( xi ), h   Ch p , константа C не h j 0 Лекция 6. Стр. 6 Величина max i 1 1   j y( xi 1 j )  xi , y( xi ), h  называется погрешностью аппроксимации. h j 0 Говорят, что численный метод решения задачи Коши сходится с порядком p, если max y ( xi )  y i  Ch p константа C не зависит от h. Величина y ( xi )  y i называется i локальной погрешностью метода. Численный метод решения задачи Коши называется устойчивым, если для всех достаточно малых h справедливо неравенство   1 i max yi  y *i  C  X  h max   j yi 1 j  xi , y *i , h   y 0  y *0  . Здесь константа  i h j 0  i   C ( X ) не зависит от h, а звездочкой помечены приближенные решения возмущенной задачи дискретной задачи Коши. Доказано, что для устойчивой дискретной задачи Коши порядок аппроксимации и порядок сходимости совпадают. Этот порядок называется порядком метода. Легко видеть, что описанный выше метод Эйлера является устойчивым методом первого порядка, а локальная погрешность метода Эйлера — величина  h 2 . Численные методы решения задачи Коши различаются используемыми формулами численного дифференцирования и способами аппроксимации правой части уравнения.   Решение задачи Коши методом Рунге-Кутты Методом Рунге—Кутты обычно называют одношаговый метод четвертого порядка, относящийся к широкому классу методов Рунге—Кутты. В этом методе величины yi 1 вычисляются по следующим формулам: h k1  2k 2  2k 3  k 4 , 6 h h  h h    k1  f ( xi , y i ), k 2  f  xi  , y i  k1 , k 3  f  xi  , y i  k 2 , k 4  f  xi  h, y i  hk3 . 2 2  2 2    y i 1  y i  Построить расчетные формулы метода Рунге-Кутты можно различными способами. Например, можно использовать следующий прием. Положив yxi 1   yxi   xi 1  f x, y( x)dx , имеем xi yxi 1   yxi   xi 1  f x, y( x)dx . Затем, заменив интеграл в правой части специально xi построенной квадратурной формулой (вычислить интеграл иначе невозможно, поскольку он зависит от неизвестного искомого решения y(x)). xi 1 Если заменить интеграл формулой левых прямоугольников  f x, y( x)dx hf ( x , y( x )) , i i xi получим метод Эйлера. xi 1 Если заменить интеграл формулой трапеций  f x, y( x)dx  2 ( f ( x , y( x )  f ( x h i i i 1 , y( xi 1 )) , xi получим y  xi 1   y  xi   h ( f ( xi , y ( xi )  f ( xi 1 , y ( xi 1 )) . Этот метод имеет второй порядок 2 точности, но является неявным. Заменим значение yi 1 на «предсказываемое» методом Эйлера. Получаем метод, который называют методом Эйлера-Коши: Лекция 6. Стр. 7 y  xi 1   y  xi   h ( f ( xi , y ( xi )  hf ( xi , y ( xi )) . 2 Метод Эйлера-Коши относятся к классу методов прогноза и коррекции. Эта идея лежит в основе многих численных методов решения задачи Коши. Отличаются методы друг от друга используемыми квадратурными формулами. В качестве упражнения полезно попытаться получить таким способом расчетные формулы метода Эйлера. Используя более или менее точные квадратурные формулы можно получить расчетные формулы разных порядков точности. Однако в практических научных расчетах формулы Рунге-Кутты порядка выше четвертого используются крайне редко из-за громоздкости и высоких требований к гладкости правой части. Метод Рунге-Кутты устойчив. Локальная погрешность метода равна Mh 4 . Практически оценить величину M достаточно сложно. При оценке погрешности обычно используют правило Рунге. Для этого сначала проводят вычисления с шагом h , а затем — с h h (h / 2) (h ) . Если y i — приближение, вычисленное с шагом h , а y 2 i — с шагом , то 2 2 16 ( h / 2) ( h / 2) (h)  y ( x 2i )  y 2i  yi . справедлива оценка y 2i 15 ( h / 2) (h) y 2i  yi h За оценку погрешности с шагом принимают величину max . i 2 15 шагом Используя оценку по правилу Рунге, можно строить адаптивные алгоритмы Рунге-Кутты. В адаптивных алгоритмах на каждом очередном шаге производят половинное деление шага до тех пор, пока max i y 2i ( h / 2)  yi 15 (h) не станет меньше заданной погрешности решения. Чтобы получить расчетные формулы метода Рунге—Кутты для систем дифференциальных уравнений, достаточно в расчетных формулах для уравнений первого порядка заменить y, f x, y , k1 , k 2 , k3 , k 4 на Y , F x, Y , k 1 , k 2 , k 3 , k 4 . Замечательная особенность методов Рунге-Кутты состоит в том, что при высокой точности он не требует вычисления производных правой части уравнения. Ведь известно, что численное дифференцирование следует применять с большой осторожностью.  