Амосов А.А., Дубинский В.А., Копченова Н.В. Вычислительные методы для инженеров (1994) (1095853), страница 73
Текст из файла (страница 73)
с Ь,у с .... В этом случае формула (13.37) приводит к следующему методу уточнения, кото.„рый называют также ристодори экстраполяции Ричардсона1. Пусть шаг ' Л измельчается по правилу Л/ — — Ь1 1/2, 1' = 1, 2, ..., Х. Сначала полагают 1Ь = 1". Для вычисления всех последующих приближений исполь- о 2 — ! Пример 13.8. Так как для формулы трапециЯ равенство (13.38) имеет место, причем Ц = 2, Щ = 4, ..., йу = 2У, то к ней можно применить экстраполяцию Ричардсона. В результате получается так называемый метод Ро.аоер~а.
Существуют стандартные программы вычисления интегралов методом Ромберга. Правда, следует отметить, что эффективность этого метода не велика. Как правило, лучший результат дает применение квадратурных формул Гаусса или рассматриваемых ниже адаптивных процедур численного интегрирования. 3 а м е ч а н и е. Первый же шаг метода Ромберга приводит к уточнению квадратурной формулы трапеций, совпадающему с формулой Симпсона. Действительно, ~й Х, + ~„т~и-1.
~ ~ + ~„ 4. Адаптивные процедуры численною интегрирования. До сих пор нам было удобнее рассматривать методы численного интегрирования с постоянным шагом. Однако они обладают значительно меньшей эффективностью и используются в вычислительной практике существенно реже, чем методы с переменным шагом. Объясняется это тем, что, равномерно распределяя узлы по отрезку интегрирования, мы полностью игнорируем особенности поведения подынтегральной функции. В то же время интуитивно ясно, что на участках плавного изменения функции достаточно поместить сравнительно небольшое число узлов, разместив значительно большее их число на участках резкого изменения функции. Распределение узлов интегрирования в соответствии с характером поведения подынтегральной функции часто позволяет при том же общем числе узлов получить значительно более высокую точность.
Тем не менее выбор соответствующего неравномерного распределения узлов интегрирования является очень сложной задачей и вряд ли мог быть широко использован на практике, если бы для решения этой задачи не удалось привлечь ЭВМ. Современные процедуры численного интегрирования (адаптпивпые хоадратурпме про1раммы) используют некоторый алгоритм автомати- 398 ь ~1" — Щ,т)йх] < а. а (13.39) тг Здесь 4' — приближенное значение интеграла 1; = [ 1(т)с1х,вычислях; ~ емое по некоторой формуле. Типичная программа разбивает исходный отрезок [а, 6] на элементарные отрезки [х; ~, х;] так, чтобы для погрешностей В; = 1; — ф выполнялось неравенство Е !Я;~ <е. (13.40) При этом каждый из элементарных отрезков получается, как правило, делением пополам одного из отрезков, найденных на более раннем шаге алгоритма. Заметим, что неравенство (13.40) выполняется, если каждая из погрешностей В; удовлетворяет условию ~Л,~ < Ь;а/($ — а).
а п е а действительно, тогда ~1 — 1"~ = ~ Е (1, — .Ьг)] 4 Е ~я,) < — Е Ь, = 1= ' ' $=1 ' 6-д 1=1 ' 6 Рассмотрим, например, одну из простейших адаптивных процедур, основанную на формуле трапеций (13.41) и использующую для контроля точности составную формулу трапеций с шагом Ь;/2: 4 ~ 1~~' = ~Д-~ + 2Я - ~ >г + Д). Ь; (13.42) Заметим, что если значение 1"' . уже найдено, то для нахождения тр,( значения 1Ь требуется лишь одно дополнительное вычисление функции 1 в точке т, . ~,2.
399 ческого распределения узлов интегрирования. Пользователь такой программы задает отрезок [а, 5], правило вычисления функции 1 и требуемую точность а > О. Программа стремится, используя по возможности минимальное число узлов интегрирования, распределить их а так, чтобы найденное значение 1Ь = Е 1~' удовлетворяло неравенству й1 Можно показать, что для оценки погрешностей квад ратурных формул интерполяционного типа на элементарных отрезках правило Рунге сохраняет силу. Поэтому погрешность приближенного значения ~ можно оценить по формуле (13.43) В рассматриваемой адаптивной процедуре последовательно выбирают точки х; и вычисляют значения последнее из которых Я„совпадает с Ф и принимается за приближенное значение интеграла 1.
Перед началом работы полагают Яс —— О, хс = а и задают некоторое начальное значение шага Ь~. Опишем в-й шаг процедуры в предположении, что значение Я; ~ уже найдено и очередное значение шага Ь; определено. 1с. По формулам (13.41) — (13.43) вычисляют значения ф,, ф еь 2с. Если ~ я;~ > Ь,яДЬ вЂ” а), то шаг Ь, уменьшают в 2 раза и повторяются вычисления и, 1э. Зс. После того как очередное дробление шага приводит к выполнению условия ~с;~ ~ Ь,яДЬ вЂ” а), вычисляют значения г, = ю, ~ + Ьь Я, = =Ь;,+ф 4с.
Если х; + Ь; > 6, то полагают Ь;,~ — 6 — г,. В противном случае полагают Ь;,~ —— Ь,. На этом ~-й шаг завершается. 3 а м е ч а н и е 1. В некоторых адаптивных процедурах при выполнении условия типа ~с;~ < Ь,г/($ — а) очередное значение шага удваивается: Ь;,~ — 2 Ь,. Таким образом, шаг интегрирования 'в зависимости от характера поведения подынтегральной функции может не только измельчаться, но и укрупняться. 3 а и е ч а н и е 2.
Известно Я, что при оптимальном распределении узлов интегрирования модули погрешностей, приходящихся на элементарные отрезки интегрирования, должны быть примерно одинаковыми. Распределение узлов, полученное с помощью адаптивных процедур, как правило, не является таковым. Однако для большинства функций оно является вполне удовлетворительным.
400 $13.5. Вычисление интегралов в нерегулярных случаях Нередко приходится вычислять интегралы /Г (2)йт а (13.44) 2 Пример 13.9. Функция — е * имеет особенность в точке с = 0. Попробуем применить для вычисления интеграла г 1 = / — е ~ Йс в 1.689677 О1- формулу прямоугольников с постоянным шагом: (13А5) а 1 2 1 в фр — Х е ~1-1 12 11 1=1 Й -1/2 (13.46) Результаты вычислений для нескольких значений шага й приведены в табл. 13.3. Заметим, что значения фр сходятся к 1 очень медленно. Теоретический анализ погрешности показывает, что она убывает пропорционально лишь а1 ~2 а не а2 как в регулярном случае. Т а б л и ц а 13.3 .байр 27 101 1.9 101 1,4 101 9.6 ° 10 2 0.200 0.100 0.050 0.025 1.420 1.499 1,555 1.594 401 от функций, имеющих те или иные особенности.
Например, сама функция Р (или ее производная некоторого порядка) имеет участки резкого изменения, точки разрыва или является неограниченной. Такие функции плохо аппроксимируются многочленами и поэтому для вычисления соответствующих интегралов может оказаться неэффективным непосредственное применение стандартных квадратурных формул, рассчитанных на возможность кусочно-многочленной аппроксимации функции Г. Случай, когда в интеграле (13.44) промежуток интегрирования бесконечен, также требует специального рассмотрения, Укажем на некоторые подходы к вычислению интегралов в нерегулярных случаях, позволяющие учесть особенности поведения функции Г и благодаря этому значительно сократить затраты машинного времени или достичь большей точности. 1.
Разбиение промеясутка интегрирования на части. Пусть подынтегральная функция Г является кусочнс~гладкой и сс < ср < ... < ср— известные точки разрыва функции Г либо ее производных. В этом случае имеет смысл представить интеграл (13.44) в виде суммы: с /Г (х) с) х = / Г (х)с) х + ( Г (х) с1 х + ... + / Г (х)А . (13.47) с 1 00 ных интегралов вида ] Г (х)дх. Если требуется вычислить такой интега рал с точностью е > О, то его представляют в виде суммы: ОЭ ь Оо ]Г (х)дх = ]Г (х)с$х + ~Г (х)с1х. Затем благодаря выбору достаточно сс я большого Ь добиваются выполнения неравенства ~ ~Г (х)с1х] < е/2 и 6 вычисляют интеграл ] Г (х)с1х с точностью е/2.
а 2. Выделение веса. В некоторых случаях подынтегральная функция допускает разложение на два сомножителя: Г (х) = р (х) /(х), где функция р (х) является достаточно простой и имеет те же особенности, что и Г (х), а /(х) — гладкая функция. Тогда имеет смысл рассматривать интеграл (13.44) в виде Ь Х = ~р (х)/(х)с1х. а (13.48) Здесь функция р (х) называется весовой функцией (или веео.а).
При построении численных методов вычисления интеграла (13.48) весовая 402 Вычисление каждого из входящих в сумму (13.47) интегралов представляет собой стандартную задачу, так как на каждом из частичных отрезков [а, с1], (с1, сх], ..., (с „, Ь] подынтегральная функция является гладкой. Разбиение промежутка интегрирования на части может оказаться полезным приемом и в других случаях, например тогда, когда имеет смысл на разных частях применять различные квадратурные формулы.
Другой пример дает стандартный прием вычисления несобствен- г1 = Е ( р(г)Рм, (х)йх а=ь га-.,\ приводит к следующей квадратурной формуле интерполяционного типа: 1 м 16 = Е Й; Е аф(х; ~ д + 1~й;/2). ,у=о (13.49) Здесь коэффициенты а," вычисляются по формуле 1 Я3 ~ .ф а11 = — Г р (х; и + щ2) и — Ы 2 6=о $ 6~1 Если функция 1" имеет на отрезке 1а, Ь1 непрерывную производную порядка т + 1, то для погрешности формулы (13.49) верна оценка !1 16! ~ СаМ .1 l! р (х) ! 1хЬтаХ о Пример 13.10. Выведем аналог формулы прямоугольников с постоянным шагом для вычисления интеграла ~ Эдмунд Никола Лагерр (1834 — 1886) — французский математик.