численные методы. лекции чм. 1. Постановка задачи Коши
 Скачать 4.45 Mb.
|
|
5. Формулы Рунге-Ромберга и экстраполяция по Ричардсону Пусть - значение численного решения задачи Коши (1.1) - (1.2), полученное в точке методом Рунге-Кутта р-го порядка точности. Если функция f непрерывно дифференцируема на на K K K K y 1 2 3 4 1 0 1 1 0 05 1 0 05 11 0 05 1 0 055 1105 0 1 1 0 1105 1 2105 1 0 1 6 1 2 11 2 1105 1 2105 111034 = + = = + + = = + + = = + + = = + + + + = , . , , , ( * . * . ) , K f x y K f x h y h K K f x h y h K K f x h y hK 1 1 1 2 1 1 1 3 1 1 2 4 1 1 3 0 1 11103 1 2103 2 2 0 1 05 11103 0 0605 1 3208 2 2 1 3258 1 4432 = = + = = + + = + + + = = = + + = = + + = ( , ) , ( , ) ( , ) , ( , ) y 2 11103 01 6 12103 2 13208 2 13258 14432 124281 = + + + + = (. * . * . ) y h k ( ) x k 13 отрезке [a, b] p+1 раз, то в любой точке вида (2.4) имеет место формула: , ( 5.1 ) где v(x) - некоторая функция от х, не зависящая от h. Уменьшив шаг в два раза, получим . ( 5.2 ) Вычитая (5.2) из (5.1), находим ( 5.3 ) или Из (5.2) и (5.3) следует, что с точностью до имеет место приближенная формула Рунге-Ромберга : . ( 5.4 ) Обозначив через левую часть в (5.4), для р=1 и р=4 соответственно имеем: . ( 5.5 ) Из (5.2) и (5.3) видно, что если , то правая часть в (5.4) имеет в точности р-й порядок относительно h и отличается от главной части погрешности, т.е. от ,на величину более высокого порядка относительно h. Формулу Рунге-Ромберга (5.4)рекомендуется применять только когда выполнено условие: . ( 5.6 ) x k y x y h v x h O h k k k p p ( ) ( ) ( ) ( ) = + + +1 y x y h v x h O h k k k p p ( ) ( ) ( ) = æ èç ö ø÷ + æ èç ö ø÷ + + 2 2 1 y h y h v x h O h k k k p p p 2 2 2 1 1 æ èç ö ø÷ - = æ èç ö ø÷ - + + ( ) ( ) ( ) ( ) ( ) v x h y h y h O h k p k p k p p ( ) ( ) 2 2 2 1 1 æ èç ö ø÷ = æ èç ö ø÷ - - + + O h p ( ) +1 y x y h y h y h k k k k p ( ) ( ) - æ èç ö ø÷ » æ èç ö ø÷ - - 2 2 2 1 e h k 2 æ èç ö ø÷ e h y h y h и e h y h y h k k k k k k 2 2 2 2 15 æ èç ö ø÷ » æ èç ö ø÷ - æ èç ö ø÷ » æ èç ö ø÷ - ( ) ( ) v x k ( ) ¹ 0 v x h k p ( ) / ( / ) 2 2 2 2 1 01 p k k k k y h y h y h y h ( ) ( ) ( ) - æ èç ö ø÷ - - < 14 Неравенство ( 5.6 ) может не выполняться по следующим причинам: 1)h слишком велико, 2) h слишком мало, 3)v(x) близко к нулю ( или равно нулю ). Умножим теперь обе части (5.2) на и из полученного равенства вычтем (5.1). В результате получим: . (5.7) Величина (5.8) называется уточненным(или экстраполированным) по Ричардсону значением величины Из (5.7) и (5.8) следует, что ; вместе с тем, согласно (5.1), если то , т.е. решение оказывается на порядок более точным, чем решение . При экстраполировании по Ричардсону также целесообразно проверять условие (5.6). 6.Численное решение систем дифференциальных уравнений первого порядка. Пусть дана система двух дифференциальных уравнений первого порядка: ( 6.1 ) Решением системы ( 6.1 ) называется пара функций , при подстановке которых в систему ( 6.1 ) получаются тождества: 2 p ( ) 2 1 2 2 1 p k p k k p y x y h y h O h - = æ èç ö ø÷ - + + ( ) ( ) ( ) ! ( ) ( ) y h y h y h k p k k p = æ èç ö ø÷ - - 2 2 2 1 y h k ( ) y x y h O h k k p ( ) ! ( ) ( ) - = +1 v x k ( ) ¹ 0 y x y h O h k k p ( ) ( ) ( ) - = ! ( ) y h k y h k ( ) ¢ = ¢ = ì í î y f x y y y f x y y 1 1 1 2 2 2 1 2 ( , , ); ( , , ). y x y x 1 2 ( ), ( ) ¢ º ¢ = y x f x y x y x y x f x y x y x 1 1 1 2 2 2 1 2 ( ) ( , ( ), ( )); ( ) ( , ( ), ( )) Y1 Y2 Y20 Y10 A0 15 Рис. 5 Любому решению системы ( 6.1 ) соответствует кривая в пространстве (x, , которая называется интегральной кривой. Для выделения единственного решения системы ( 6.1 ) фиксируют точку , как показано на рис.5. Числа называют начальными данными. Задача Коши для системы ( 6.1 ) ставится так: требуется найти решение системы ( 6.1 ), удовлетворяющее начальному условию , ( 6.2 ) Теорема. Пусть в некоторой области D пространства (x, непрерывны функции , а также частные производные этих функций , . Тогда для любой точки D система ( 6.1 ) имеет единственное решение, удовлетворяющее начальному условию ( 6.2 ). Рассмотрим теперь систему n дифференциальных уравнений ( 6.3 ) при начальных условиях . ( 6.4 ) Постановка задачи Коши ( 6.3 ) - ( 6.4 ) и теорема существования и единственности решения формулируются аналогично предыдущему ( в случае n = 2 получаются приведенные выше формулировки ). Введем векторы y y 1 2 , ) A x y y 0 0 10 20 ( , , ) y y x y y x 10 1 0 20 2 0 = = ( ) , ( ) x y y 0 10 20 , , y y x x 1 10 0 = = y y x x 1 20 0 = = y y 1 2 , ) f x y y 1 1 2 ( , , ) f x y y 2 1 2 ( , , ) ¶ ¶ ¶ ¶ f y f y 1 1 1 2 , ¶ ¶ ¶ ¶ f y f y 2 1 2 2 , ( , , ) x y y 0 10 20 Î ¢ = ¢ = ¢ = ì í ï ï î ï ï y f x y y y y f x y y y y f x y y y n n n n 1 1 1 2 2 2 1 2 3 1 2 ( , , ,...., ) ( , , ,...., ) ( , , ,...., ) y y y y y y x x x x n x x n 1 10 2 20 0 0 0 0 = = = = = = , , ......, X0 X 16 и запишем задачу ( 6.3 )-( 6.4 ) в векторной форме: . ( 6.5 ) Для решения задачи ( 6.5 ) численными методами разобьем отрезок [a, b], на которой ищется решение, на N равных частей узловыми точками: . ( 6.6 ) Отметим, что Метод Эйлера для задачи ( 4.5 ) определяется формулами , ( 6.7 ) где вектор берется из начального условия, а - векторы приближений к точному решению в узловых точках вида ( 6.6 ). Соотношения ( 6.7 ) могут быть выписаны покоординатно; при n=2, т.е. для задачи ( 6.1 )-( 6.2 ), формулы ( 6.7 ) представимы в виде ( 6.8 ) Метод Эйлера-Коши для задачи ( 6.5 ) может быть записан в векторной форме следующим образом: ( 6.9) В случае одного уравнения ( n = 1 ) из формул ( 6.9 ) получаются (4.4) . Аналогично , классический метод Рунге-Кутта для задачи Коши ( 6.5 ) имеет вид: , ( 6.10 ) где y y y y y x y x y x y x y x y x y x y x f x y f x y f x y f x y n n n n 0 10 20 0 1 2 1 2 1 2 = æ è ç ç ç çç ö ø ÷ ÷ ÷ ÷÷ = æ è ç ç ç çç ö ø ÷ ÷ ÷ ÷÷ = ¢ ¢ ¢ æ è ç ç ç çç ö ø ÷ ÷ ÷ ÷÷ = æ è ç ç ç çç ö ø ÷ ÷ ÷ ÷÷ , ( ) ( ) ( ) ( ) , '( ) ( ) ( ) ( ) , ( , ) ( , ) ( , ) ( , ) ¢ = = = y f x y y y x x ( , ), 0 0 x a kh k N k = + = , ,,..., 0 1 x a x b h b a N N 0 = = = - , , ( ) / y y h f x y k N k k k k + = + = - 1 0 1 1 ( , ) , , , . . . y 0 y y 1 2 , ,.... y x ( ) x x 1 2 , ,... y y hf x y y y y hf x y y k N k k k k k k k k k k 1 1 1 1 1 2 2 1 2 2 1 2 0 1 1 , , ( , , ); ( , , ); ,,... + + = + = + = - ì í ï îï y y h f x y f x y h f x y k N k k k k k k k k + + = + + + = - 1 1 2 0 1 1 [ ( , ) ( , ( , ) ) ] , , , . . y y h K K K K k N k k + = + + + + = - 1 1 2 3 4 6 2 2 0 1 1 ( ) , , , . . . , 17 ( 6.11 ) В случае n=2 , т.е. для задачи ( 6.1 ) - ( 6.2 ), формулы ( 6.10 ) в координатном виде записываются так: ( 6.12 ) k = 0,1,……,N-1. Здесь коэффициенты -координаты векторов из формул ( 6.11 ) видно, что коэффициенты зависят от номера шага k. Погрешность решения оценивается по формуле, аналогичной ( 5.8 ). А именно, пусть - значения численного решения в точке , полученные для шагов h и h/2 соответственно. Тогда погрешность решения в точке , найденного методом Рунге-Кутта ( 6.10 ), оценивается по формуле Рунге - Ромберга: ( 6.13 ) Здесь , ,где - точное решение задачи Коши ( 6.5 ) в точке В случае, когда применяется метод Эйлера, множитель 1/15 в ( 6.13 ) отсутствует. Пример 3. Решить методом Рунге-Кутта задачу Коши ( 6.14 ) K f x y K f x h y h K K f x h y h K K f x h y hK k k k k k k k k 1 2 1 3 2 4 3 2 2 2 2 = = + + = + + = + + ( , ), ( , ); ( , ), ( , ) y y h K K K K y y h K K K K k k k k 1 1 1 11 12 13 14 2 1 2 21 22 23 24 6 2 6 2 2 , , ( ) ; ( ) , + + = + + + + = + + + + K ij K i j j ( , ; , , , ); = = 1 2 1 2 3 4 K ij y h y h y h y h y h y h y h y h k k k nk k k k nk ( ) ( ) ( ) ( ) , ( / ) ( / ) ( / ) = æ è ç ç ç çç ö ø ÷ ÷ ÷ ÷÷ æ èç ö ø÷ = æ è ç ç ç çç ö ø ÷ ÷ ÷ ÷÷ 1 2 1 2 2 2 2 2 x k x k { } e y h y h k i n ik ik » - £ £ 1 15 2 1 max ( ) ( / ) e y x y h k k k k = - ( ) ( ) { } e y h y x k i n ik i k » - £ £ 1 15 1 max ( ) ( ) y x k ( ) x k ¢ = ¢ = - ì í î y y y y 1 2 2 1 , , y x y x 1 2 0 1 2 0 1 5 = = = = - . , 18 на отрезке [0, 0.5] с шагом 0.1. Сравнить результаты численного решения с точным решением. Решение. По условию, и, значит, При k=0,h=0.1 по формулам (6.11)находим Теперь по формулам (4.12) получаем Продолжая этот процесс, находим и т.д. Результаты вычислений представлены в таблице 4 .Точные значения вычислены по формулам Таблица 4 f x y y y f x y y y x y y 1 1 2 2 2 1 2 1 0 10 20 0 1 2 0 5 ( , , ) , ( , , ) , , . , = = - = = = - f x y y y y x y ( , ) , = - æ è ç ö ø ÷ = = = - æ è ç ö ø ÷ 2 1 0 0 1 2 0 5 ( ) K f x y y y K K 1 0 0 20 10 11 21 0 5 1 2 0 5 1 2 = = - æ è ç ö ø ÷ = - - æ è ç ö ø ÷ = - = - , ; . , . , K f x h y h K y h K y h K K K K f x h y h K y h K y 2 0 0 1 20 21 10 11 12 22 3 0 0 2 20 22 2 2 2 2 0 5 0 05 1 2 1 2 0 05 0 5 0 56 1175 0 56 1175 2 2 2 = + + æ èç ö ø÷ = + - - æ è ç ç çç ö ø ÷ ÷ ÷÷ = - + - - - - æ è ç ö ø ÷ = = - - æ è ç ö ø ÷ = - = - = + + æ èç ö ø÷ = + - , . ( . ) . ( . ) ; . , , , ( ) 10 12 13 23 4 0 0 3 20 23 10 13 2 0 5 0 05 1175 1 2 0 05 0 56 0 56875 1172 0 55875 117200 0 5 0 1 1172 1 2 0 1 0 55875 - æ è ç ç çç ö ø ÷ ÷ ÷÷ = - + - - - - æ è ç ö ø ÷ = - - æ è ç ö ø ÷ = - = - = + + = + - - æ è ç ö ø ÷ = - + - - - - æ è ç ö h K K K K f x h y hK y hK y hK . ( . ) . ( . ) ; , , , . ( . ) . ( . ) ø ÷ = = - - æ è ç ö ø ÷ = - = - 0 6172 114413 0 61720 114413 14 24 ; , K K y y 11 21 1 2 0 1 0 5 2 0 56 2 0 55875 0 6174 6 114409 0 5 0 1 1 2 2 1175 2 1172 114413 6 0 617302 = + - - × - × - = = - + - - × - × - = - ì í î . ( . ) / ; . ( . ) / y y y y 12 22 13 23 , , , y x y x k k 1 2 ( ) , ( ) y x x x y x x x 1 2 1 2 0 5 1 2 0 5 ( ) . cos . si n ; ( ) . si n . cos . = - = - - ì í î 19 0.1 0.2 0.3 0.4 0.5 1.14409 -.617302 1.07675 -.728436 .998644 -.832292 .910564 -.927832 .813387 -1.01410 1.14409 -.617302 1.07675 -.728436 .998644 -.832293 .910564 -.927832 .813386 -1.01410 Из таблицы видно, что в данном примере 7.Численное решение дифференциальных уравнений высших порядков Задача Коши для дифференциального уравнения n-го порядка ставится так: требуется решить уравнение (7.1) при начальных условиях (7.2) Пусть функция f(x) является решением задачи (7.1)-(7.2). Обозначим (7.3) Тогда и, следовательно, . Таким образом, задача (7.1)-(7.2) приводится к задаче Коши для системы дифференциальных уравнений первого порядка: (7.4) с начальными условиями (7.5) Здесь числа такие же как и в (7.2), а функции выражаются через у(х) по формулам (7.3). Пример 4. Решить методом Эйлера задачу Коши . (7.6) X K Y K 1 Y K 2 Y X Y X K K 1 2 ( ) ( ) e k » - 10 6 y f x y y y n n ( ) ( ) ( , , ' , . . . . , ) = -1 y x x y y x x y y x x y n n = = = = ¢ = = - - 0 0 0 0 1 0 0 1 , ' , . . . , ( ) ( ) y x y x y x y x y x y x y x y x n n 1 2 1 3 2 1 ( ) ( ) , ( ) ( ) , ( ) ( ) , . . . , ( ) ( ) = = ¢ = ¢ = ¢ - y x y x y x y x y x y x n n 2 3 1 ( ) ( ) , ( ) ( ) , . . . , ( ) ( ) ( ) = ¢ = ¢¢ = - ¢ = y y x n n ( ) ( ) ¢ = ¢ = ¢ = ¢ = ì í ï ïï î ï ï ï - y y y y y y y f x y y y n n n n 1 2 2 3 1 1 2 ( , , , . . . . , ) ( ) y x x y y x x y y x x y n n 1 0 0 2 0 0 0 0 1 = = = = ¢ = = - , , . . . . , ( ) y y y n 0 0 0 1 , , . . . , ( ) ¢ - y x y x y x n 1 2 ( ) , ( ) , . . . . , ( ) xy x y y y x y x ¢¢ - + ¢ + = = = ¢ = = - ( ) , , 1 0 1 1 1 1 20 на отрезке [1, 1.5] с шагом h=0.1. Сравнить результаты численного решения с точным решением. Решение. Положим Задача (5.6) приводится к следующей задаче Коши для системы двух дифференциальных уравнений первого порядка: Решим эту задачу методом Эйлера. В формулах (6.8) положим По условию h=0.1. При k=0 для из (6.8) получим Аналогично для находим и т.д. Точное решение задачи (7.6) имеет вид y(x)=2(x+1)-3exp(x-1). Результаты вычислений представлены в таблице 5 . Таблица 5 1.1 1.2 1.3 1.4 1.5 0.9 0.77 0.607 0.4077 0.16847 0.88449 0.73579 0.55042 0.32526 0.05386 Здесь к = 1,2, ....,5. Пример 5. Решить задачу Коши (7.7) (7.8) на отрезке[1, 1.5] с шагом h=0.1. Оценить погрешность найденного решения. Решение. Положим Задача Коши (7.7)-(7.8) приводится к задаче Коши для системы дифференциальных уравнений: y x y x y x y x y x y x y x 1 2 1 2 ( ) ( ) , ( ) ( ) ( ) , ( ) ( ) = = ¢ = ¢ ¢ = ¢¢ ¢ = ¢ = + æ èç ö ø÷ - ì í ï îï y y y x y y x 1 2 2 2 1 1 1 ; , y x y x 1 2 1 1 1 1 = = = = - , ( ) f x y y y f x y y x y y x x y y 1 1 2 2 2 1 2 2 1 0 10 20 1 1 1 1 1 ( , , ) ; ( , , ) / / , , , = = + - = = = - x x h 1 0 11 = + = . y y hf x y y y y hf x y y 11 10 1 0 10 20 22 20 2 0 10 20 1 0 1 1 0 9 1 0 12 1 1 1 3 = + = + - = = + = - + × - - = - ( , , ) . ( ) . ; ( , , ) . ( ( ) ) x x h 2 1 1 2 = + = . y y hf x y y y y hf x y y 21 11 1 1 11 12 22 21 2 1 11 12 0 2 0 1 1 3 0 77 1 3 0 1 3 3 163 = + = + - = = + = - + - = - ( , , ) . ( . ) . ; ( , , ) . ( . ) X K Y K Y X K ( ) y y y x x x x k k k k k k k = = + - - = + 1 2 1 3 1 1 0 1 , ( ) ( ) exp( ) , x y xy y I V 2 4 2 0 + ¢¢¢ + ¢¢ = y x y x y x y x = = ¢ = = ¢¢ = = ¢¢¢ = = 1 1 0 1 1 2 , y x y x y x y x y x y x 1 2 1 3 ( ) ( ) , ( ) ( ) ( ) , ( ) = = ¢ = ¢ = = ¢ = ¢¢ = ¢ = ¢¢¢ ¢ = y x y x y x y x y x y x y x I V 2 4 3 4 ( ) ( ) , ( ) ( ) ( ) , ( ) ( ) 21 с начальными условиями Решим эту задачу методом Рунге-Кутта. Таблица 6 1.1 1.2 1.3 1.4 1.5 .0102833 .0419953 .0958918 .1722490 .2710400 .0102833 .0420011 .0958984 .1722560 .7104700 Результаты вычислений приведены в таблице 6 ; там же указаны значения точного решения в данных узловых точках. Точное решение задачи (7.7)-(7.8): y(x) = 10 - 10x + 4ln|x| + 6ln|x| можно найти, если заметить, что уравнение (7.7) записывается в виде . Из таблицы получаем (k=1,2,....,5) Отметим, что оценить погрешности можно было также с помощью формулы (5.5). |