Д.П. Костомаров, А.П. Фаворский - Вводные лекции по численным методам (pdf) (1113729), страница 21
Текст из файла (страница 21)
(73)∂v∂v ⎝2α2α⎠ ⎝ 2α ∂v⎠Подставим полученные выражения для отдельных слагаемых в формулу (70). Врезультате она примет вид рекуррентной формулыzi +1 = {1 + hF ′ ( ui + θ i zi )} zi − ψ i h , 0 ≤ i ≤ n − 1 ,(74)которую нужно дополнить нулевым начальным условием (38). Использовать этуформулу для последовательного вычисления значений сеточной функции z нельзя: вее правую часть входят неизвестные аргументы: ψ i , ui , θ i .
Однако эту системурекуррентных равенств можно заменить системой рекуррентных неравенств дляпоследующей оценки zi .- 101 -Предположим, как и при исследовании метода Эйлера, что частная производная∂f( x, u ) в интересующей нас области изменения ее аргументов ограничена (40).
Тогда∂uс учетом формулы (73) для производной F ′ ( v ) получим11 + hF ′ ( ui + θ i z ) ≤ 1 + Ch + C 2 h 2 < eCh = q , q > 1 .(75)2С учетом этого рекуррентные равенства (74) можно заменить рекуррентныминеравенствами(76)zi +1 ≤ q zi + ψ c h ,которые полностью совпадают с неравенствами (42) предыдущего раздела. Мы ужезнаем, что из них следует оценка нормы погрешности решения через нормупогрешности аппроксимации уравнения(77)z c ≤ leCl ψ c .Теперь нужно оценить норму погрешности аппроксимации уравнения (68).Предположим, что функция f ( x, u ) имеет в интересующей нас области изменениясвоих аргументов непрерывные вторые производные и, следовательно, решениедифференциального уравнения u ( x ) трижды непрерывно дифференцируемо. Этопозволяет написать следующие разложения Тейлора11ui +1 = u ( xi + h ) = ui + u′ ( xi ) h + u′′ ( xi ) h 2 + u′′′ ( xi ) h 3 ,(78)26hhh ⎧ ∂f∂f⎫⎛⎞f ⎜ xi +, ui +f ( xi , ui ) ⎟ = f ( xi , ui ) +⎨ ( xi , ui ) + ( xi , ui ) f ( xi , ui ) ⎬ +2α2α2α ⎩ ∂x∂u⎭⎝⎠(79)⎫h2 ⎧ ∂2 f∂2 f∂2 f2+ 2 ⎨ 2 ( x%i , u%i ) + 2( x%i , u%i ) f ( x%i , u%i ) + 2 ( x%i , u%i ) f ( x%i , u%i ) ⎬ ,8α ⎩ ∂x∂x∂u∂u⎭гдеhhxi = xi + θ i h , x%i = xi + θ%i, u%i = ui + θ%if ( xi , ui ) ,2α2α0 < θ i < 1 , 0 < θ%i < 1 .Здесь последние слагаемые в обоих разложениях представляют собой остаточныечлены в форме Лагранжа, которые берутся в неизвестных нам промежуточных точках.Подставим разложения (78), (79) в формулу (68) для погрешностиаппроксимации дифференциального уравнения (27) и примем во вниманиесоотношения, вытекающие из этого уравненияu ′ ( xi ) = f ( xi , ui ) ,(80)∂f∂fu ′′ ( xi ) = ( xi , ui ) + ( xi , ui ) f ( xi , ui ) .∂x∂uБлагодаря (80) члены нулевого и первого порядков относительно h сокращаются иостаются только члены второго порядка, обязанные своим происхождениемостаточным членам в разложениях (78), (79).
В результате получается следующеепредставление для погрешности аппроксимации уравнения- 102 -⎧11ψ i = h ⎨ u′′′ ( xi ) −8α⎩62⎡ ∂2 f⎢ ∂x 2 ( x%i , u%i ) +⎣(81)⎫⎤∂ f∂ f+2( x%i , u%i ) f ( x%i , u%i ) + 2 ( x%i , u%i ) f 2 ( x%i , u%i )⎥ ⎬ .∂x∂u∂u⎦⎭Функции, входящие в правую часть этого соотношения, по предположениюнепрерывны и ограничены в интересующей нас области изменения своих аргументов.Это позволяет заменить равенство (81) неравенством(82)ψ i ≤ ψ c ≤ Mh 2 ,где M - константа, мажорирующая выражение в фигурных скобках формулы (81).Подставляя оценку (82) в неравенство (77), получим(83)z c ≤ MleCl h 2 .Таким образом, при h → 0 погрешность аппроксимации уравнения и, какследствие, погрешность решения стремятся к нулю со скоростью h2 . Это означает, чторазностное уравнение (61), полученное по схеме Рунге-Кутта, имеет второй порядокточности относительно h .Второй порядок точности лучше, чем первый, однако практика показывает, чтоэтой точности также недостаточно.
Наиболее часто при проведении реальныхрасчетов используется схема Рунге-Кутта четвертого порядка точности следующеговидаyi +1 − yi 1= ( k1 + 2k2 + 2k3 + k4 ) ,(84)h6гдеhh ⎞⎛k1 = f ( xi , yi ) , k2 = f ⎜ xi + , yi + k1 ⎟ ,22 ⎠⎝(85)hh ⎞⎛k3 = f ⎜ xi + , yi + k2 ⎟ , k4 = f ( xi + h, yi + hk3 ) .22 ⎠⎝Если в схеме второго порядка точности на каждом шаге функцию f ( x, y )приходилось вычислять два раза, то здесь – четыре раза. Однако это усложнениесхемы расчета окупается высокой точностью.
На более подробном обсуждении схемы(84), (85) останавливаться не будем и ограничимся конкретным примером.Задача 4.Построить решение задачи Коши (51), (52) на отрезке [0, 2] с шагом h = 0.25 посхеме Рунге-Кутта второго порядка типа «предиктор-корректор» (65) и по схемеРунге-Кутта четвертого порядка (84), (85). Сравнить результаты расчетов междусобой и с аналитическим решением задачи (53).Результаты расчетов по этой задаче приведены в таблице 2. Здесь в первомстолбце даны значения переменной xi , во втором и третьем столбцах – результатырасчетов по схемам Рунге-Кутта второго и четвертого порядков, в последнем шестомстолбце – значения аналитического решения (53) в узлах сетки. В четвертом и пятом22- 103 -столбцах приведены результаты расчетов по методу Адамса.
Они будут обсуждаться вследующем разделе.Таблица 2.xiu ( xi )Р.К. - II Р.К. - IVАд. - IIАд. - IV0,00 0,000000 0,000000 0,000000 0,000000 0,0000000,25 0,031250 0,032593 0,031250 0,032593 0,0325940,50 0,133057 0,136099 0,130859 0,136099 0,1361020,75 0,314791 0,319962 0,309692 0,319962 0,3199661,00 0,587068 0,594879 0,578331 0,594826 0,5948851,25 0,961913 0,972975 0,948662 0,972847 0,9729841,50 1,452948 1,467988 1,434141 1,467772 1,4680001,75 2,075605 2,095486 2,050001 2,095159 2,0955012,00 2,847365 2,873107 2,813492 2,872644 2,873127Сравнение результатов второго столбца таблицы 1, рассчитанных по методуЭйлера с шагом h = 0.25 , с результатами второго и третьего столбца таблицы 2показывает как уменьшается погрешность при фиксированном шаге h по мереперехода к более точным методам.
Так метод Рунге-Кутта четвертого порядка,несмотря на достаточно крупный шаг, дает погрешность решения z c = 0.00002 . Этона много лучше, чем при расчете по схеме Эйлера с шагом h = 0.01 (см. четвертыйстолбец в таблице 1). В то же время при расчете по схеме Эйлера было сделано двестишагов с однократным вычислением функции f ( x, y ) на каждом шаге, а при расчете посхеме Рунге-Кутта – восемь шагов с четырехкратным вычислением функции f ( x, y )на каждом шаге. Таким образом, более сложный, но и более совершенный методпозволяет при меньшем объеме вычислений получить более точный результат.В заключение сделаем следующее замечание.
Априорные оценки погрешностипо схеме Эйлера (50) или Рунге-Кутта (83) представляют теоретический интерес. Ониопределяют скорость, с которой погрешность стремится к нулю при h → 0 . Однако напрактике оценки подобного типа неэффективны, поскольку содержат производныеискомого решения u ( x ) . Обычно точность численного решения задачи устанавливаютс помощью апостериорных оценок, основанных на сравнении результатов расчетов сшагом h и h / 2 . Процедура их вывода и применения была описана в предыдущейглаве в связи с задачей численного интегрирования.2.4. Метод Адамса.Адамс – английский астроном и математик XIX века, который много занималсянебесной механикой. При изучении траекторий планет ему постоянно приходилосьчисленно интегрировать уравнения их движения. Желая минимизировать объемвычислений, Адамс разработал один из наиболее экономичных методов численногорешения дифференциальных уравнений, к обсуждению которого мы теперьпереходим.Пусть u ( x ) - решение дифференциального уравнения (27).
Для производнойэтой функции имеет место равенство- 104 -u ′ ( x ) = f ( x, u ( x ) ) = F ( x ) .(86)Интегрируя его между двумя точками сетки, получим соотношениеui +1 = ui +xi +1∫ F ( x ) dx .(87)xiМы не можем использовать это соотношение непосредственно для перехода впроцессе решения задачи от i -ой точки сетки к ( i + 1) -ой, поскольку функция F ( x )нам не известна. Чтобы сделать следующий шаг, нужно приближенно заменить этуфункцию на такую функцию, которую можно вычислить.
Опишем, как эта проблемарешается в методе Адамса.Пусть в процессе численного решения задачи мы довели расчет до точки xi . Врезультате проведенных расчетов нам оказались известными величины y j и f ( x j , y j ) ,0 ≤ j ≤ i . Возьмем некоторое фиксированное целое число m ≤ i и построиминтерполяционный многочлен m -ой степени, принимающий в точках x j , i − m ≤ j ≤ iзначения f ( x j , u j )Pm ( x j ) = f ( x j , u j ) , i − m ≤ j ≤ i .(88)Его можно записать по формуле ЛагранжаPm ( x ) =i∑ f ( x , y )Q ( x) ,j =i − mjj(89)m, jгде Qm , j ( x ) специальные многочлены видаQm , j ( x ) =( x − xi−m )L ( x − x j−1 )( x − x j+1 )L ( x − xi )( x j − xi−m )L( x j − x j−1 )( x j − x j+1 )L( x j − xi ),(90)которые мы уже рассматривали в третьей главе.Главная идея метода Адамса заключается в том, чтобы для расчета yi +1использовать формулу типа (87), приближенно заменяя в ней функцию F ( x ) наинтерполяционный многочлен Pm ( x ) , составленный согласно (89) по результатампредыдущих вычислений.
Это приводит к рекуррентной формулеyi +1 = yi +xi +1∫Pm ( x ) dx = yi +xii∑ a f (x , y ),j =i − mjjj(91)гдеaj =xi +1∫ Q ( x ) dx .m, j(92)xiРассмотрим более подробно данную схему численного решения задачи Коши впростейших случаях m = 0 и m = 1 , когда технические трудности не закрываютпрозрачную идею метода.
При m = 0 для аппроксимации функции F ( x ) используетсяполином нулевой степени, т. е. постояннаяF ( x ) ≈ P0 = f ( xi , yi ) .- 105 -В этом случае формула (91) переходит в рекуррентную формулу метода Эйлераyi +1 = yi + hf ( xi , yi ) ,обеспечивающую первый порядок точности. Такой результат сам по себе тривиален.Мы привели его только для того, чтобы показать, что для метода Адамса, как и дляметода Рунге-Кутта, исходной точкой является схема Эйлера.Перейдем к исследованию варианта m = 1 . В этом случае для аппроксимациифункции F ( x ) используется полином первой степени, построенный по значениямфункции f в двух точках ( xi −1 , yi −1 ) и ( xi , yi ) :x − xi −1x − xiP1 ( x ) = f ( xi , yi )− f ( xi −1 , yi −1 ).hhПодставляя его в формулу (91) и проводя интегрирование, получим1⎧3⎫yi +1 = yi + ⎨ f ( xi , yi ) − f ( xi −1 , yi −1 ) ⎬ h .(93)22⎩⎭Отметим следующую особенность рекуррентной формулы (93).
Для расчетаочередного значения сеточной функции yi +1 нужно знать ее значения в двухпредыдущих точках yi и yi −1 . Таким образом, формула (93) начинает работать толькосо второй точки. Вычислить по ней y1 нельзя. Это значение решения разностнойзадачи приходится вычислять каким-нибудь другим методом, например, методомРунге-Кутта.Рекуррентную формулу (93) можно записать в виде разностного уравненияyi +1 − yi 31= f ( xi , yi ) − f ( xi −1 , yi −1 ) .(94)22hПодсчитаем для него погрешность аппроксимации дифференциального уравненияu − u ⎧311⎫ u − u ⎧3⎫ψ i = i +1 i − ⎨ f ( xi , ui ) − f ( xi −1 , ui−1 ) ⎬ = i +1 i − ⎨ u′ ( xi ) − u′ ( xi −1 )⎬ .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
