численные методы. лекции чм. 1. Постановка задачи Коши
 Скачать 4.45 Mb.
|
1 2 1. Постановка задачи Коши. Дифференциальное уравнение первого порядка, разрешенное относительно производной , имеет вид: ( 1.1 ) Решением дифференциального уравнения ( 1.1 ) называется функция y(x), подстановка которой в уравнение ( 1.1 ) обращает его в тождество: График решения y = y(x) называется интегральной кривой. Например, если f(x ,y)= y, то интегральные кривые описываются уравнениями вида: Задача Коши для дифференциального уравнения ( 1.1 ) состоит в том, чтобы найти решение уравнения ( 1.1 ), удовлетворяющее начальному условию: ( 1.2 ) Рис.1 Числа - называются начальными данными. Условие (1.2) означает, что . Выбор начального условия служит для выделения одной из интегральных кривых данного уравнения ( рис.1 ). Например, решением задачи Коши y’ = y; y(0) = 1 является функция Условия существования и единственности решения задачи Коши содержатся в следующей теореме. Теорема. Пусть функция f(x,y) непрерывна вместе со своей частной производной в области D. Тогда для всякой точки уравнение ( 1.1 ) имеет единственное решение у(х), удовлетворяющее начальному условию ( 1.2 ). y f x y ' ( , ) = y x f x y x ' ( ) ( , ( ) ) . = y Ce X = y x x y = = 0 0 x y 0 0 , y x y ( ) 0 0 = y e X = ¶ ¶ f y / ( , ) x y 0 0 Î D Y X x y 0 0 3 Некоторые дифференциальные уравнения решаются стандартными аналитическими методами. Для решения линейного дифференциального уравнения 1-го порядка y’+ p(x)y = q(x) ( 1.3 ) сначала находят общее решение соответствующего однородного уравнения y’+ p(x)y = 0 ( 1.4 ) по формуле , ( 1.5 ) а затем применяют метод вариации произвольной постоянной. Для этого заменяют константу С в формуле ( 1.5 ) функцией С(х) так, чтобы получить решение данной задачи Коши. Пусть, например, требуется получить решение задачи Коши y’- y = x, y(0) = 1 ( 1.6 ) По условию ,p(x)=-1,q(x)=x, Соответствующее однородное уравнение имеет вид: y’-y=0, а его общее решение находится по формуле ( 1.5 ): Применим метод вариации произвольной постоянной: После подстановки в исходное уравнение получим Отсюда Следовательно, общее решение данного дифференциального уравнения имеет вид По условию у=1 при х=0. Следовательно, С = 2. Решением данной задачи Коши является функция ( 1.7 ) 2. Понятие о численном решении задачи Коши Для многих задач Коши, представляющих практический интерес, невозможно найти точное решение. Приближенные решения задачи Коши находятся численными методами. Предположим, что точное решение задачи Коши (1.1)-(1.2) существует, единственно и обладает необходимыми свойствами гладкости. Это решение будет обозначаться у(х). Задача численного решения задачи Коши (1.1)-(1.2) состоит в том, чтобы получить таблицу приближенных значений искомого решения у(х)для заданных дискретных значений аргумента х: y C p t d t X X = - ò exp( ( ) ( ) ) 0 x y 0 0 0 1 = = , y C dt Ce X X = - - = ò exp( ( ) ) 1 0 y Cx e y C x e Cx e X X X = = + ( ) , ' ' ( ) ( ) C x e x C x xe X X ' ( ) , ' ( ) = = - C x xe dx xe e C X X X ( ) = = - - + - - ò y Cx e Ce x X X = = - - ( ) 1 y e x X = - - 2 1 4 ( 2.1 ) Точки ( 2.1 ) называют узловыми точками, а образованное ими множество называют сеткой. Начальная точка берется из условия ( 1.2 ), так что значение известно. Полученное численным методом приближенное значение для будет обозначаться , так что , k = 1, 2, …., N. ( 2.2 ) Обычно рассматривают равномерную сетку с шагом h , т.е. сетку для которой выполнено условие .( 2.3 ) Полагая , при условии ( 2.3 ) будем иметь . ( 2.4 ) где h = (b-a)/N. Узловые точки (2.4) делят отрезок [a ,b] на N частичных отрезков длины h. Приближенные значения находятся последовательно: сначала по известному значению вычисляется (первый шаг), затем находится (второй шаг)и т.д. В численном решении задачи Коши различают ошибки дискретизации и ошибки округления. Ошибка дискретизации есть свойство используемого метода. Это значит, что если бы все арифметические операции выполнялись абсолютно точно, то других ошибок, кроме ошибки дискретизации, не было бы. Ошибка округления зависит от ЭВМ и программы. Важно различать между собой два вида ошибок дискретизации: локальную ошибку дискретизации и глобальную ошибку дискретизации. Пусть - функция удовлетворяющая условиям , ( 2.5 ) т.е. есть решение дифференциального уравнения ( 1.1 ), определяемое не исходным начальным условием ( 1.2 ), а значением вычисленного решения в точке . Тогда локальная ошибка дискретизации определяется равенством , где - вычисленное приближенное значение для Локальная ошибка дискретизации - это ошибка, сделанная на данном шаге, при условии, что предыдущие значения точны и нет ошибок округления. Глобальная ошибка дискретизации есть разность между вычисленным решением (ошибки округления игнорируются) и точным решением: , Отметим, что , x x x N 0 1 , ,.... x 0 y x y ( ) 0 0 = y x k ( ) y k y y x k k » ( ) x x x x x x h N N 1 0 2 1 1 0 - = - = = - = > - a x b x N = = 0 , x a kh k N k = + = , ,,...., 0 1 y y y N 1 2 , ,...., y y x 0 0 = ( ) y y x 1 1 » ( ) y y x 2 2 » ( ) u x k ( ) u x f x u x u x y k k k k k ' ( ) ( , ( )), ( ) º = u x k ' ( ) x k d k d y u x k k k k = - + + 1 1 ( ) y k +1 y x k ( ) +1 e k e y x y k k k = - ( ) u x y x e d 0 1 0 ( ) ( ) , = = 5 Пусть численное решение задачи Коши (1.1)-(1.2) находят в узловых точках вида ( 2.4 ). Тогда погрешностью метода называют вектор Величину погрешности оценивают нормой Если шаг h достаточно мал и длина отрезка [a,b] не очень велика, то Численный метод называют сходящимся, если , при (тогда необходимо ). Говорят, что метод имеет р-й порядок точности, если существует число С такое, что ( 2.6 ) Число С может зависеть от функции f(x ,y), от начальных данных и длины отрезка [a,b], на котором ищется решение; вместе с тем, число С не должно зависеть от величины шага h. Напомним, что если - некоторая функция от h , то равенство означает, что отношение ограничено при . Поэтому условие ( 2.6 ) кратко записывается в виде 3. Метод Эйлера Простейшим численным методом решения задачи Коши (1.1)- (1.2) является метод Эйлера. Этот метод определяется формулами , ( 3.1 ) где узловые точки имеют вид ( 2.4 ) Формулы ( 3.1 ) получаются следующим образом. Пусть у(х) - точное решение задачи ( 1.1 ) - ( 1.2 ). По формуле Тейлора , ( 3.2 ) где точка расположена между и . Согласно ( 1.1 )- ( 1.2 ) имеем . Ввиду ( 2.4 ) имеем . Значит, из формулы ( 3.2 ) при получаем: . ( 3.3 ) Отсюда , ( 3.4 ) т.е. получилась первая из формул ( 3.1 ) . Заметим, что где , x k e h e e e N ( ) ( , ,...., ) = 1 2 { } e h e k N K ( ) max : , , . . = = 1 2 e d d d N N » + + + - 0 1 1 e h ( ) ® 0 h ® 0 N ® ¥ e h Ch P ( ) £ x y 0 0 , j( ) h j( ) ( ) h O h P = j( ) / h h P h ® 0 e h O h P ( ) ( ) = y y hf x y k N K K K K + = + = - 1 0 1 1 ( , ), ,,..., x k y x y x y x x x y x x ( ) ( ) ' ( ) ( ) ( ) ! ( ) = + - + ¢¢ - 0 0 0 0 2 2 x x x x 0 ¢ = = y x f x y y x y ( ) ( , ) , ( ) 0 0 0 0 0 x x h 1 0 - = x x = 1 y x y hf x y h y ( ) ( , ) ( ) 1 0 0 0 2 1 2 = + + ¢¢ x y x y y hf x y ( ) ( , ) 1 1 0 0 0 » = + 1 2 2 2 h y Ch ¢¢ £ ( ) , x { } C y x a x b = ¢¢ £ £ 1 2 max ( ) : 6 т.е при переходе от ( 3.3 ) к ( 3.4 ) было отброшено слагаемое вида Если в формуле ( 3.2 ) вместо взять , то с учетом равенств вместо( 3.4 ) получим Продолжая этот процесс, придем к формулам ( 3.1 ). Метод Эйлера имеет наглядный геометрический смысл. Рис. 2 Пусть АВ - интегральная кривая у=у(х), . Угловой коэффициент касательной к кривой АВ в начальной точке А равен .Прямая пересекает эту касательную в точке , ордината которой вычисляется по формуле ( 3.4 ). Через точку проходит интегральная кривая , где функция удовлетворяет условиям ( 2.5 ) при к = 1. Угловой коэффициент касательной к кривой в точке равен . Прямая пересекает эту касательную в точке , ордината которой вычисляется по формуле ( 3.5 ). Через точку проходит интегральная кривая , где функция удовлетворяет условиям ( 2.5 ) при к = 2. Повторяя проведенные построения, получим точку , а затем точки . В результате интегральная кривая АВ аппроксимируется ломанной Эйлера ( рис. 2 ). Oh ( ) 2 x 0 x 1 ¢ » » - = y x f x y y x y x x h ( ) ( , ) , ( ) , 1 1 1 1 1 2 1 y x y y hf x y ( ) ( , ) 2 2 1 1 1 » = + a x b £ £ ¢ = y x f x y ( ) ( , ) 0 0 0 x x = 1 A 1 y 1 A 1 y u x = 1 ( ) u x 1 ( ) y u x = 1 ( ) A 1 ¢ = u x f x y 1 1 1 1 ( ) ( , ) x x = 2 A 2 y 2 A 2 y u x = 2 ( ) u x 2 ( ) A 3 A A A N 4 5 , ,...., AA A A A N 1 2 3 7 Отметим на кривой АВ точки c координатами ( и положим . Глобальная ошибка дискретизации равна длине отрезка . На рис.2 видны также локальные ошибки дискретизации . Например, ошибка совпадает с длиной отрезка . При выводе формулы ( 2.4 ) было показано, что т.е. на первом шаге метода Эйлера допускается ошибка порядка Аналогично , если - функция, удовлетворяющая условиям ( 2.5 ) , то , Отсюда Здесь N = (b-a)/h. Следовательно, Аналогичные равенства верны для Поэтому . ( 3.6 ) Это значит, что метод Эйлера имеет первый порядок точности. На практике это проявляется в том, что при уменьшении шага h приближенное решение будет все более точным и при стремлении h к нулю будет сходиться к точному решению с линейной скоростью по h. В частности, при уменьшении шага h вдвое ошибка уменьшится примерно в два раза. Проиллюстрируем последнее обстоятельство на примере задачи Коши Используя различные шаги h , методом Эйлера найдем приближенные значения при х=1 ( табл.1 ). Число шагов N = 1/h. Таблица 1 h = 1.0 2.00 0.718 0.65 1/2 2.25 0.468 0.59 1/4 2.441 0.277 0.55 1/8 2.566 0.252 0.53 1/16 2.638 0.080 Точное решение при равно у(1)=2.718… Приближенное значение в точке , полученное при вычислении с шагом h , обозначено . В правом столбце приведено отношение B K x y x k N K K , ( ))( , ,..., ) = - 1 2 1 B B N = e y x y K K k = - ( ) A B K K d y u x 1 2 1 2 = - ( ) A C 2 2 d y y x O h 0 1 1 2 = - = ( ) ( ) h 2 u x K ( ) d y u x O h K K K K = - = + + 1 1 2 ( ) ( ) e d d d N Oh N N » + + + = + 0 1 1 2 ( ) e O h N = ( ) e e e N 1 2 1 , ,...., - e h Oh ( ) ( ) = ¢ = = = y y y X , | 0 1 y h N ( ) e h N ( ) y y h N ( ) ( ) 1 - e h e h N N ( / ) / ( ) 2 x N = 1 x N = 1 y h N ( ) 8 ошибок при последовательном делении пополам шага h. Из таблицы видно, что это отношение близко к 0.5 . Пример 1. Решить методом Эйлера задачу Коши , приняв шаг h=0.1 и N=4. Оценить погрешность метода в данной задаче. Решение. Узловые точки: По условию, f(x ,y) = x + y и,следовательно, формулы ( 3.1 ) принимают вид: Учитывая, что h=0.1 , , последовательно находим: , Точное решение данной задачи Коши определяется выражением ( 1.7 ). Округляя до шести значащих цифр, получим y(0.1)=1.11034, y(0.2)=1.24280, y(0.3)=1.39972, y(0.4)=1.58365. Отсюда видно, что На практике точное решение задачи Коши известно редко. Поэтому для определения достигнутой точности поступают следующим образом. Так как то начиная с некоторого h, первые значащие цифры перестают меняться - они отвечают точному решению. Произведя несколько расчетов с различными h и сравнив результаты, получают представление о достигнутой точности. Таблица 2 h y1 y2 y3 y4 0.1 1.10000 1.22000 1.36200 1.52820 0.01 1.10924 1.24038 1.3657 1.57773 0.001 1.11023 1.24256 1.39931 1.58305 0.0001 1.11033 1.24278 1.39968 1.58359 В табл. 2 приведены результаты численного решения методом Эйлера задачи Коши ( 1.6 )при четырех различных ¢ = + = = y x y y x , | 0 1 a x x x x x b = = = = = = = 0 1 2 3 4 0 01 0 2 0 3 0 4 , . , . , . , y y h x y k K K K K + = + + = 1 0 1 2 3 ( ), ,, , x y 0 0 0 1 = = , y 1 1 010 1 11 = + + = .( ) y y y 2 3 4 11 0 1 0 1 11 1 22 1 22 0 1 0 2 1 22 1 362 1 362 0 1 0 3 1 362 1 5282 = + + = = + + = = + + = . ( . . ) . , . ( . . ) , . ( . ) e( .) 01 0 06 £ y y x h k K - ® ® ( ) ,( ), 0 0 y K 9 значениях h. Приближенные значения даны для узловых точек Обсудим теперь кратко влияние ошибок округления , возникающих при компьютерной реализации метода Эйлера. Здесь имеется два источника ошибок округления. Во-первых, это ошибка, возникающая при вычислении . Обозначим эту ошибку . Во-вторых, ошибка , возникающая при вычислениях по самой формуле Эйлера. Таким образом, фактически вычисленные значения удовлетворяют соотношениям . ( 3.7 ) Во время вычислений происходит систематическое накопление ошибок, так как при вычислении исходные данные не являются точными и содержат погрешности, зависящие от неточности предшествующих вычислений. Вместе с тем, как уже отмечалось, ошибки дискретизации метода Эйлера стремятся к нулю при . Следовательно, за счет уменьшения шага, можно сделать ошибку дискретизации сколь угодно малой. Однако, чем меньше h , тем больше потребуется шагов по методу Эйлера и, вообще говоря, тем больше скажутся на полученном решении ошибки округления. На практике, когда при выполнении арифметических операций используются слова фиксированной длины, всегда существует такая величина шага h , меньше которой вклад ошибок округления начинает доминировать в суммарной ошибке. Эта ситуация изображена схематически на рисунке 3. Здесь - определяет минимальный шаг, который можно использовать в практических расчетах, Е - суммарная ошибка. В задачах, где не требуется слишком высокая точность, необходимый шаг обычно будет значительно больше, чем этот минимум, и основной вклад в полную ошибку будет вносить ошибка дискретизации. Такое поведение характерно и для других методов, хотя минимальный размер шага будет меняться от метода к методу и от задачи к задаче. 4. Методы Рунге-Кутта В формуле Эйлера ( 3.1 ) заменим некоторым средним значением функции f: , ( 4.1 ) где узловые точки имеют вид ( 2.4 ). В зависимости от способа вычисления среднего формула ( 4.1 ) определяет x x x x 1 2 3 4 01 0 2 0 3 0 4 = = = = . , . , . , f x y k k ( , ) e K h K y K [ ] y y h f x y K K K K K K + = + + + 1 ( , ) e h y K + 1 ( , ), x y k K K ³ 1 h ® 0 h 0 f x y k k ( , ) f k y y hf k N k k k + = + = - 1 0 1 1 , , ,......, x k f k 10 различные методы Рунге-Кутта. Метод Эйлера можно рассматривать как один из методов Рунге-Кутта, он получается из формулы ( 4.1 ) при = Метод Эйлера-Кошиопределяется парой расчетных формул: ( 4.2 ) ( 4.3 ) Этот метод получается из ( 4.1 ) при , т.е. - является средним арифметическим значением функции f, вычисленным в двух точках : . На отрезке интегральная кривая заменяется отрезком прямой с угловым коэффициентом, равным (рис.4). Формула ( 4.3 ) уточняет значение , найденное по формуле Эйлера (3.2) Рис.4 После подстановки ( 4.2 ) в ( 4.3 ) получим ( 4.4 ) Классический метод Рунге-Кутта определяется формулами: f k f x y k k ( , ) y y hf x y k k k k + = + 1 * ( , ) [ ] y y h f x y f x y k N k k k k k k + + + = + + = - 1 1 1 2 0 1 1 ( , ) ( , ), , ,......., * [ ] ( , ) ( , ) * f f x y f x y k k k k k = + + + 1 2 1 1 f k ( , ) , ( , ) * x y x y k k k k + + 1 1 [ ] x x k k , +1 f k y k + 1 * y y h f x y f x y hf x y k k k k k k k k + + = + + + 1 1 2 (( , ) ( , ( , )) Y A k X x k x K+1 y k Тангенс угла наклона f(X ,Y* ) k+1 k+1 Тангенс угла наклона=среднему арифметическому Тангенс угла наклона f(X , Y ) k k 11 , (4.5) где ), ( 4.6 ) Метод ( 4.5 ) получается из ( 4.1 ) при , т.е. является взвешенным средним значением функции f , вычисленным в четырех различных точках. Ввиду равенства числа являются значениями производной y’в четырех различных точках. На каждом шаге производная подсчитывается четыре раза: один раз в начальной точке, дважды в промежуточных точках и один раз в конечной точке. Затем найденные значения умножаются на соответствующие коэффициенты и прибавляются к . В результате получаются значения ( см. формулы (4.5),( 4.6 ) ). Классический метод Рунге-Кутта имеет четвертый порядок точности, т.е. для него . При уменьшении шага вдвое, ошибка метода ( 4.5 ) уменьшится примерно в раз. Пример 2. Методами Эйлера-Коши и Рунге-Кутта решить задачу Коши из примера 1. Решение. Требуется решить задачу Коши ( 2.8 ) при h=0.1, N=4 двумя способами: по формулам ( 4.2 ) - ( 4.3 )и по формулам ( 4.5 ) - ( 4.6 ). По условию, Узловые точки: По формулам ( 4.2 )-( 4.3 ) последовательно получаем и т.д. По формулам ( 4.5 )-( 4.6 ) при к=0 получим y y h K K K K k N k k + = + + + + = - 1 1 2 3 4 2 2 2 0 1 1 ( ), , ,......, K f x y K f x h y h K k k k k 1 2 1 2 2 = = + + ( , ), ( , K f x h y h K K f x h y hK k k k k 3 2 4 3 2 2 = + + = + + ( , ), ( , ) ( ) f K K K K k = + + + 1 6 2 2 1 2 3 4 f k ¢ = y f x y ( , ), K i i ( , . . . ) = 1 4 K i y k y k + 1 e h O h ( ) ( ) = 4 2 16 4 = f x y x y x y ( , ) , , = + = = 0 0 0 1 a x x x x x b = = = = = = = 0 1 2 3 4 0 0 1 0 2 0 3 0 4 , . , . , . , y y y y 1 1 2 2 1 010 1 11 1 0 05 0 1 01 11 111 111 0101 111 1231 111 0 05 01 111 0 2 1231 119205 * * .( ) . (( ) ( . .)) .( . . ) . (( . . ) ( . )) = + + = = + + + + = = + + = = + + + + = 12 Аналогично , при к=1 имеем и т.д. Результаты вычислений приведены в таблице 3 Таблица 3 x 0 0.1 0.2 0.3 0.4 y (Эйлер) 1.00000 1.10000 1.22000 1.36200 1.52820 у(Эйлер- Коши) 1.00000 1.11000 1.24205 1.39847 1.58180 у(Рунге- Кутта) 1.00000 1.11034 1.24281 1.39972 1.58365 у(Точное решение) 1.00000 1.11034 1.24280 1.39973 1.58365 Методы Рунге-Кутта построены для произвольного порядка точности p. Однако, чем выше порядок, тем больше требуется вычислений значений f на каждом шаге. В классическом методе Рунге-Кутта ( 4.5 ) число таких вычислений совпадает с порядком p и равно 4. Если 5 < p < 7 , то требуется p+1 таких вычислений, а для p > 8 число вычислений f равно p+2 . Поэтому методы Рунге-Кутта высокого порядка почти не применяются , а требуемая точность обычно достигается специальными способами коррекции шага. |