1626435587-51311eae4652e8ad616b5bdef025cbb3 (844239), страница 13
Текст из файла (страница 13)
В случае, если порядок численного метода (итерационного процесса)не известен, можно провести аналогичноерассмотрение, считая одной из неизвестных в системе (35, 36), дополненной уравнением вида = (2 ℎ) + · (2 ℎ) . Данный алгоритмизвестен каки рассмотрен в [2].Заметим также, что если проводить вычисления на последовательности сеток ℎ, ℎ, . .
., ℎ, можно добиться повышения порядка точности численного метода на .a prioriпроцесс Эйткена4.6. Интегралы по бесконечной областиЧасто в физике возникают несобственные интегралы по бесконечной области: по всему пространству, по бесконечному времени и т. п.62В данном разделе кратко рассмотрены основные приёмы, используемых для нахождения численных значений таких интегралов на ЭВМ.В случае, если подынтегральная функция убывает достаточно быстро (экспоненциально) при больших значениях аргумента, можно использовать самый простой способ — ограничение области∫︀ ∞ интегрирования. Так, пусть необходимо вычислить интеграл () и пусть∃ > 0, и 0 : | ()| ≤ − ∀ > 0 .
Заменяя верхний пределинтегрированияконечным значением ≥ 0 , мы допустим ошибку∫︀ ∞ = () . Оценим погрешность, мажорируя интеграл: || ≤∫︀ ∞ −= − . Требуя, чтобы погрешность метода не превосходила некоторого , можно оценить верхний предел интегрирования: = 1 ln(/). Видно, что уменьшение приводит к медленному (логарифмическому) росту верхнего предела интегрирования . Например,при = 1, = 1, = 10−15 получаем = 34,5.Подчеркнём, что указанный способ неприменим для функций, убывающих на бесконечности медленно (по степенному закону).
Рассмотрим | ()| ≤ − ∀ > 0 , > 1. В этом случае мажорантная оценкадля погрешности ∫︀замены верхнего предела интегрирования на будет∞иметь вид || = | () | ≤ | −1|1− . Накладывая условие || ≤ ,получаем1⃒⃒ −1⃒⃒⃒⃒=⃒.( − 1) ⃒Для = 1, = 2, = 10−8 имеем = 108 . Большая величина полученного значения обуславливает необходимость использования более эффективных методов интегрирования медленно убывающих функций.Для вычисления несобственных интегралов функций, имеющих набесконечности степенну́ю асимптотику, обычно используют замену переменной.
Например, делая замену = /(1 − ), получаем∫︁∞=∫︁1 () = (()) .(1 − )200В результате имеем интеграл по конечной области, однако как правило,он также является несобственным ввиду наличия особенности в подынтегральном выражении при → 1−0. Способам вычисления интеграловтакого типа посвящён следующий пункт.634.7. Интегралы с особенностью∫︀ Пусть требуется вычислить значение интеграла () в конечных пределах, при этом функция () может иметь интегрируемыеособенности в точках и .В случае, если () имеет лишь устранимые особенности, для вычисления интеграла достаточно переопределить функцию в особыхточках либо в∫︀ их малых окрестностях.
Например, пусть требуется вы1числить = 0 sin . Подынтегральная функция имеет устранимуюособенность в точке = 0, поэтому в программном коде следует доопределить подынтегральную функцию в нуле:double f(double x) { return (x == 0) ? 1 : sin(x)/x; }Если проверка условия x==0 даст логическую истину, функция f вернёт1, в противном случае будет вычислено отношение sin(x)/x.В основу указанного решения положено сравнение двух чисел с плавающей точкой на точное равенство — неиссякаемый источник ошибокчисленных алгоритмов, см. п. 4 на с.
39. Как следствие, описанный способ может плохо работать или даже приводить к сбоям в программе вслучае, когда погрешность вычисления числителя дроби существенна,а знаменатель близок к нулю. Более надёжным решением является замена подынтегральной функции суммой первых членов ряда Тейлора вокрестности устранимой особой точки. Для рассмотренного выше примера имеем:sin 246=1−+−+ (8 ).6120 5040√Откуда делаем вывод, что при || < 6 значение подынтегральнойфункции будет неотличимо от 1 (здесь — машинное эпсилон), что позволяет безопасно вычислять подынтегральную функцию следующимобразом:double f(double x){const double x0 = 3.6e-8; //3.6e-8 ~ sqrt(6*DBL_EPSILON)return (fabs(x) < x0) ? 1 : sin(x)/x;}Рассмотрим теперь другой случай: пусть функция () имеетнеустранимую интегрируемую особенность, так что | | неограниченновозрастает в окрестности конечного множества точек на промежуткеинтегрирования.
Без ограничения общности можно считать, что ()64имеет особенность только на концах промежутка — в противном случаеследует разбить отрезок на несколько частей, представив интеграл ввиде суммы интегралов по отрезкам с концами в особых точках подынтегральной функции.Для вычисления интегралов в указанных условиях существует рядприёмов. Например, можно использовать аддитивное выделение особенности, представляя () = () + (), где |()| < ∞, а () интегрируется аналитически. Другой способ состоит в мультипликативномвыделении особенности: () = ()·(); полученный интеграл с весовой функцией может быть вычислен с помощью квадратурных формулГаусса-Кристоффеля [2].
Недостатком этого и ряда других приёмов,описанных в [2], является необходимость поиска нестандартного решения для каждой задачи. В главе 4.5 книги [4] рассмотрен простой алгоритм, свободный от указанного недостатка. Он позволяет вычислятьв общем виде несобственные интегралы и обладает при этом высокойскоростью сходимости. Суть его состоит в следующем. Пусть требует∫︀ ся вычислить интеграл = () .
Сделаем под интегралом заменупеременной:11 = ( + ) + ( − ) · th ,(39)22Тогда+∞∫︁∫︁ − = (()) = 1 + (()) ·· 2 ,2ch −∞(40)−где th и ch — гиперболический тангенс и косинус , соответственно.Интегрирование в (40) ведется в бесконечных пределах, однако подынтегральная функция быстро убывает ввиду наличия множителя ch−2 ,что позволяет заменить пределы интегрирования на конечные, от −до , и вычислить полученный интеграл с помощью любой квадратурной формулы (например, формулы Симпсона).
Заметим, что погрешность при вычислении интеграла будет складываться из двух составляющих, = 1 + 2 . Первая связана с отбрасыванием экспоненциальноубывающих бесконечных «хвостов», вторая — с заменой интеграла в конечных пределах на значение квадратурной формулы. Увеличивая ,можно сделать ошибку 1 сколь угодно малой. Однако при фиксированном числе точек увеличение будет приводить к росту величинышага сетки ℎ = 2 /( − 1) и, как следствие, к увеличению остаточного члена 2 в квадратурной формуле. Напротив, уменьшая шаг сеткиℎ, можно уменьшить остаточный член 2 , однако при = const это65приведёт к уменьшению и росту 1 .
Следовательно, должен существовать некоторый компромисс между указанными двумя источниками погрешности, позволяющий минимизировать величину суммарнойошибки = 1 + 2 . Условие 1 ≈ 2 приводит к оценке оптимальнойвеличины шага [4].(41)ℎ≈ √2Заметим, что для различных подынтегральных функций, обладающихразличной расходимостью в точках = , , коэффициент в (41) может варьироваться, поэтому выражение (41) даёт, как правило, субоптимальную величину шага. Однако ввиду быстрой сходимости рассмотренного метода и удобства использования универсальных квадратурных формул, оптимизацию в каждом отдельном случае обычно невыполняют.Также заметим, что хотя замена (39) даёт < () < ∀, из-заограниченной точности машинных вычислений при достаточно больших || возможно выполнение условий () = или () = .
Посколькуподынтегральная функция может иметь особенность в точках = , ,вычисление в этих точках будет приводить к сбоям в программе илиполучению в качестве результата не-числовых значений inf или NaN.Следовательно, для корректной работы алгоритма необходимо отслеживать в программе выполнение условий < () < .4.8. Сравнение методовВ данном разделе мы рассмотрели лишь наиболее простые квадратурные формулы, ограничившись лишь случаем равномерных сеток(формулы). На практике широко используется формула Симпсона, реже применяются формулы,имеющий первый алгеб(︀ )︀раический порядок точности (Σ = ℎ2 ), тогда как формулы нулевого порядка практически не используются ввиду крайне низкой скорости сходимости. Исключением является вычисление нормы функций наравномерных сопряжённых сетках, связанных преобразованием Фурье.В этом случае формулы правых,левых и центральных прямоуголь∑︀ников и трапеций дают 1| |2 .
Указанная сумма сохраняется придискретном преобразовании Фурье (равенство Парсеваля), что делает использование формул прямоугольников предпочтительнее методаСимпсона для данной задачи.При интегрировании функций, имеющих гладкие производные высоких порядков, могут применяться квадратурные формулы более высоких порядков точности. Наиболее употребительными являются квад-Ньютона-Котеса66ратурные формулы Гаусса:∫︁ () ≈∑︁ ( ),=1где весовые коэффициенты и расположение узлов , = 1, .
. . , выбираются с целью достижения максимально возможного22 порядкаточности [A10] при заданном числе узлов сетки . Кроме того, формулы Гаусса являются устойчивыми независимо от , тогда как коэффициенты в формулах Ньютона-Котеса высокого порядка точности( ≥ 10) становятся знакопеременными [A10], что приводит к резкомуросту погрешности квадратурной формулы (численной неустойчивости). Анализ свойств формул Гаусса с доказательством соответствующих теорем можно найти в монографии [A10].
В книге [4] приведёнготовый к использованию код программы на языке Си для вычисленияинтегралов по десяти точкам с использованием квадратурной формулыГаусса.Подчеркнём, что использование формул высокогоещё неозначает получение высокойчисленного решения. Применение формул высокого порядка даёт выигрыш в точности только длядостаточно гладких функций, которые притом хорошо аппроксимируются полиномами.Помимо значений подынтегральной функции, квадратурные формулы могут включать также значения её производных [2], что позволяет повысить их точность, однако усложняет вычисления.Заметим также, что хотя приведённые в данном разделе квадратурные формулы и позволяют решить подавляющее большинство практических задач, в некоторых случаях эффективным решением являетсяиспользование нестандартных формул.
Это связано,∫︀ прежде всего, свычислением интегралов с весовой функцией вида ()() . Весовая функция () может иметь особенность или, например, быстро осциллировать на фоне медленно меняющейся амплитуды ().В этом случае применение специальных квадратурных формул позволяет существенно повысить эффективность расчётов за счёт использования более крупного шага численного интегрирования.
Дляряда часто встречающихся весовых функций готовые решения можно∫︀ найти в литературе. Например, для вычисления интегралов вида() exp() используются[2].точностиформулы Филонапорядка22 Алгебраический порядок точности формул Гаусса равен 2 − 1, так что остаточный член составных формул Σ = (ℎ2 ).67Упражнения1. Полагая функцию () представимой в виде суммы ряда Тейлорана промежутке интегрирования [, ], показать, что если некото∫︀ ∑︀рая квадратурная формула () = ( ) + имеет алгебраический порядок точности , то остаточный член этой формулы = (ℎ+2 ), а остаточный член построенной на её основесоставной формулы Σ = (ℎ+1 ).a priori2. Оценитьчисло точек, в которых необходимо вычислитьподынтегральную функцию для нахождения значения функциистандартного нормального распределения Φ0,1 (2) с точностью нехуже 0,1 % с использованием формулы Симпсона:1Φ0,1 () = √2∫︁(︂ 2 )︂exp −.2(42)−∞∫︀ /43. Вычислить значение интеграла 0 sin различными методами.
Построить графики апостериорной оценки погрешности интегрирования в зависимости от числа точек , в которых вычислялась подынтегральная функция. При построении графикаиспользовать двойной логарифмический масштаб (set logscalexy в gnuplot). Сравнить с графиками степенных функций, построенных в тех же осях, объяснить результат.4.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
