Д.П. Костомаров, А.П. Фаворский - Вводные лекции по численным методам (pdf) (1113729), страница 22
Текст из файла (страница 22)
(95)h2h2⎩2⎭⎩2⎭Предположим, что функция f ( x, u ) имеет в интересующей нас области измененияаргументов непрерывные вторые производные, так что решение задачи u ( x ) триждынепрерывно дифференцируемо. Запишем разложения Тейлора11ui +1 = ui + u′ ( xi ) h + u′′ ( xi ) h 2 + u′′′ xi + θ% i h h3 ,26(96)1%u′ ( xi −1 ) = u′ ( xi ) − u′′ ( xi ) h + u′′′ xi − θ% i h h 2 .2Подставляя их в формулу (95), получим1⎧1⎫ψ i = ⎨ u′′′ xi + θ%i h + u′′′ xi − θ%%i h ⎬ h2 .(97)64⎩⎭Отсюда можно написать оценку5ψ i ≤ ψ c ≤ M 3h 2 ,(98)12где M 3 - постоянная, мажорирующая третью производную функции u ( x ) :())(()()- 106 -u′′′ ( x ) ≤ M 3 , x0 ≤ x ≤ x0 + l .(99)Мы видим, что разностное уравнение метода Адамса, соответствующее случаю m = 1 ,аппроксимирует дифференциальное уравнение (27) со вторым порядком точностиотносительно h .
Как и в случае метода Рунге-Кутта, это обеспечивает второй порядокточности для погрешности решения z c при предположении, что значение y1 , котороерассчитывается нестандартно, вычислено со вторым порядком точности.Процесс построения более точных схем можно продолжить за счет увеличенияm . При m = 2 получается схема третьего порядка точности, при m = 3 - четвертого ит.д.
Схема четвертого порядка, как и в методе Рунге-Кутта, является наиболееупотребительной, поэтому мы коротко остановимся на ее выводе и обсуждении.Если написать интерполяционный полином третьей степени P3 ( x ) (89) на сеткеиз четырех точек xi , xi −1 , xi −2 , xi −3 и провести интегрирование (92), то рекуррентнаяформула (91) примет вид:59379⎧ 55⎫yi +1 = yi + h ⎨ f ( xi , yi ) −f ( xi −1 , yi −1 ) +f ( xi −2 , yi −2 ) −f ( xi −3 , yi −3 ) ⎬ . (100)242424⎩ 24⎭Приведем еще одно форму записи этой формулы через так называемые конечныеразности153yi +1 = yi + hfi + h 2 ∆1 fi + h3∆ 2 fi + h 4 ∆ 3 fi ,(101)2128гдеf i = f ( xi , yi ) ,1{ f ( xi , yi ) − f ( xi−1, yi−1 )} ,h(102)1∆ 2 f i = 2 { f ( xi , yi ) − 2 f ( xi −1 , yi −1 ) + f ( xi −2 , yi −2 )} ,h1∆ 3 fi = 3 { f ( xi , yi ) − 3 f ( xi −1 , yi −1 ) + 3 f ( xi −2 , yi −2 ) − f ( xi −3 , yi −3 )}.hПервая, вторая и третья разности (102) приближенно соответствуют первой, второй итретьей производной функции F ( x ) = f ( x, u ( x ) ) .
Эквивалентность формул (100) и(101) легко проверить непосредственно. Формула (101) иногда более удобна дляорганизации вычислительного процесса и контроля точности.Особенность метода Адамса проявляется в формуле (100) еще сильнее, чем вформуле (93). Здесь для расчета очередного значения yi +1 нужно знать значения y вчетырех предыдущих точках - yi , yi −1 , yi−2 , yi−3 . Таким образом, формула (100)начинает работать только с четвертой точки.
Вычислить по ней y1 , y2 , y3 нельзя. Этизначения решения разностной задачи приходится рассчитывать другим методом,например, методом Рунге-Кутта.Перейдем к обсуждению точности схемы (100). Если функция f ( x, u ) имеетнепрерывные четвертые производные по своим аргументам в интересующей насобласти их изменения, так что решение задачи u ( x ) пять раз непрерывно∆1 f i =- 107 -дифференцируемо, то разностное уравнение (100) аппроксимирует дифференциальноеуравнение (27) с четвертым порядком точности относительно h .
Доказательство этогоутверждения проводится также, как и для схемы второго порядка (93), только теперь вразложениях типа (96) нужно удерживать больше членов. Четвертый порядокточности при аппроксимации уравнения обеспечивает четвертый порядок точностидля погрешности решения z c при предположении, что начальные значения дляметода Адамса y1 , y2 , y3 вычислены с такой же точностью. Они рассчитываютсянезависимо и при этом важно, чтобы начальный этап вычислительного процесса невнес такую погрешность, которая исказит все последующие результаты.Задача 5.Построить решение задачи Коши (51), (52) на отрезке [0,2] с шагом h = 0.25 посхеме Адамса второго (93) и четвертого (100) порядка.
Сравнить результатырасчетов между собой, с результатами расчетов по схеме Рунге-Кутта и саналитическим решением задачи.Результаты расчетов приведены в четвертом и пятом столбцах таблицы 2. Всоответствии с заданием, нужно сравнивать четвертый столбец со вторым и шестым, апятый – с третьим и шестым. Напомним, что в шестом столбце приведеноаналитическое решение (53) рассматриваемой задачи, так что сравнение с нимпозволяет судить о точности приближенного решения по схеме Рунге-Кутта и схемеАдамса.Расчет по схеме Адамса второго порядка точности начинается с y2 , четвертого с y4 . Значение y1 в четвертом столбце, y1 , y2 , y3 в пятом столбце рассчитывались посхеме Рунге-Кутта соответствующего порядка, поэтому в таблице они оказываютсяодинаковыми с соответствующими данными второго и третьего столбцов. Сравнениерезультатов проведенных расчетов двумя методами с аналитическим решением задачипоказывает, что их точность примерно одинакова.Сравним схемы четвертого порядка точности в методе Рунге-Кутта (84) иАдамса (100) с точки зрения организации вычислительного процесса.
Чтобы сделатьодин шаг по методу Рунге-Кутта, необходимо вычислить функцию f ( x, y ) четырераза (85), а в методе Адамса только один раз. В трех предшествующих точках функцияf ( x, y ) была уже вычислена на предыдущих шагах и вычислять ее снова нетнеобходимости. В этом заключается главное достоинство метода Адамса, котороеособенно высоко ценилось в докомпьютерную эру.Главный недостаток метода Адамса мы уже отмечали: при его применениипервые шаги приходится делать с помощью другого метода, например, с помощьюметода Рунге-Кутта и только после этого можно перейти на расчет по схеме Адамса.Таким образом, программа решения задачи Коши по методу Адамса должна включатьв себя как элемент программу метода Рунге-Кутта для расчета начальной стадиивычислительного процесса.С этой особенностью метода Адамса связана еще одна проблема.
При численноминтегрировании дифференциального уравнения часто приходится менять шаг h . Вметоде Рунге-Кутта это не составляет труда, поскольку каждый шаг делается- 108 -независимо от предыдущего. В методе Адамса ситуация иная. Здесь нужно либоизначально программировать весьма сложные формулы расчета с переменным шагом,либо после каждой смены шага заново проводить расчет первых трех точек по методуРунге-Кутта. Только после этого можно переходить на стандартный счет по методуАдамса. Эти недостатки приводят к тому, что сегодня при компьютернах расчетахпредпочтение часто отдается более удобному методу Рунге-Кутта.§3. Численное решение краевой задачи для линейногодифференциального уравнения второго порядка.Рассмотрим следующую задачу для линейного дифференциального уравнениявторого порядка:(103)u′′ − q ( x ) u = − f ( x ) , a < x < b,(104)u ( a ) = u1 , u ( b ) = u2 .Здесь два дополнительных условия заданы в граничных точках отрезка [ a, b ] , поэтомузадачу (103), (104) называют краевой.Пусть функции f ( x ) и q ( x ) непрерывны на отрезке [ a, b ] , причем(105)q ( x ) ≥ q0 > 0 .При сделанных предположениях, как известно из курса дифференциальныхуравнений, решение задачи (103), (104) существует и является единственным.Перейдем к обсуждению вопросов, связанных с его расчетом с помощьючисленного метода.
Возьмем некоторое целое число n , введем шаг h = ( b − a ) / n ипостроим сеткуxi = a + ih, 0 ≤ i ≤ n .(106)Заменим дифференциальное уравнение (103) его разностным аналогом. В результатеполучим следующую задачу:yi −1 − 2 yi + yi +1− qi yi = − f i , 1 ≤ i ≤ n − 1 ,(107)h2y0 = u0 , yn = u2 .(108)Здесь qi = q ( xi ) , f i = f ( xi ) , граничные условия (108) для сеточной функции { yi } взятытакими же, что и в дифференциальной задаче.Разностные уравнения (107) можно переписать в видеyi −1 − ( 2 + qi h 2 ) yi + yi +1 = − fi h 2 , 1 ≤ i ≤ n − 1.(109)Мы получили линейную систему из ( n − 1) -го уравнения с ( n − 1) -им неизвестным yi ,1 ≤ i ≤ n −1 . Значения y0 и yn неизвестными не являются: они задаются граничнымиусловиями (108).Между разностными схемами для задачи Коши и для краевой задачи естьсущественное различие.
В первом случае для определения сеточной функции { yi } мыимели рекуррентные соотношения, которые позволяли последовательно рассчитатьвсе ее значения. Такие разностные схемы называются явными. В краевой задаче (107),- 109 -(108) сеточная функция { yi } определяется из решения системы линейныхалгебраических уравнений. Такая разностная схема называется неявной.Из записи разностных уравнений в форме (109) видно, что мы получили системууравнений с трехдиагональной матрицей с диагональным преобладанием:диагональный элемент ( 2 + qi h 2 ) больше суммы двух других элементов той же строки,равной 2. Системы такого типа мы уже встречали в третьей главе в связи с задачейинтерполяции кубическим сплайном.
Диагональное преобладание гарантируетсуществование и единственность решения системы, которое может быть построенометодом прогонки.Перейдем к обсуждению основного вопроса: с какой точностью сеточнаяфункция { yi } , полученная в результате решения задачи (107), (108), приближаетрешение краевой задачи (103), (104).
Пусть u ( x ) решение исходной краевой задачи.Обозначим через ui = u ( xi ) его значения в узлах сетки и введем две сеточныефункции: погрешность решения и погрешность аппроксимации уравненияzi = yi − ui , 0 ≤ i ≤ n ,(110)u − 2ui + ui +1− qiui + f i , 1 ≤ i ≤ n − 1 .ψ i = i −1(111)h2Выразим из соотношения (110) yi через ui и zi и подставим в разностноеуравнение (107). Оставим члены, содержащие zi , слева, а остальные члены перенесемнаправо. В результате получимzi −1 − 2 zi + zi +1⎧ ui −1 − 2ui + ui +1⎫−qz=−−qu+f(112)⎨⎬ = −ψ i , 1 ≤ i ≤ n − 1 .iiiii2h2h⎩⎭Граничные условия в дифференциальной и разностной задачах совпадают, так чтозначения сеточной функции zi в граничных точках будут нулевымиz0 = z n = 0 .(113)Мы не можем рассчитать погрешность {zi } , решая задачу (112), (113), поскольку вправые части уравнений входят неизвестные величины ui и ψ i .
Однако задача (112),(113) позволяет оценить погрешность.Пусть максимальное по модулю число zi соответствует индексу i = j :z c = z j ≥ zi , 0 ≤ i ≤ n .(114)В граничных точках zi обращается в ноль (113), так что индекс j не равен ни нулю,ни n . Рассмотрим уравнение (112) для этого значения индекса и запишем его в виде:(115)( 2 + q j h2 ) z j = z j−1 + z j+1 + ψ j h2 .Возьмем модуль от обеих частей равенства и оценим правую часть сверху( 2 + q j h 2 ) z j = ( 2 + q j h 2 ) z c ≤ z j−1 + z j+1 + ψ j h 2 ≤ 2 z c + 2 ψ c h 2илиzc≤1ψ c.q0(116)- 110 -Здесь мы сократили одинаковые члены слева и справа, разделили обе частинеравенства на множитель q j h 2 и заменили q j в знаменателе на минимальновозможное значение функции q ( x ) на отрезке [ a, b ] , равное q0 (105). Таким образомнам удалось оценить погрешность решения z c через погрешность аппроксимацииуравнения ψ c .Для оценки погрешности аппроксимации уравнения предположим, что функцииf ( x ) и q ( x ) дважды непрерывно дифференцируемы на отрезке [ a, b ] .