Учебник - Практические занятия по вычислительной математике - Аристова (1238839), страница 28
Текст из файла (страница 28)
(1.28)), это означает, что происходит увеличение порядка аппроксимации не на единицу, как следовало ожидать,а на две. Это свойство является характерной чертой формул Ньютона–Котеса четных порядков аппроксимации. Применяя данную процедуру кформулам Симпсона, можно получить вычисление интеграла с шестымпорядком аппроксимации:(2.9)I I 4 (h / 2) I 4 (h / 2) I 4 (h) 15 O(h6 ) .Повторяя эту процедуру дальше, получим алгоритм Ромбергауточнения вычисления интеграла. Этот алгоритм может рассматриватьсякак экстраполяция вычисленных значений интегралов на последовательности вдвое сгущающихся сеток в нулевой шаг интегрирования. Важно,что экстраполяция должна производиться не по величине шага h, а повеличине h2.VII.3. Квадратурные формулы ГауссаЕсли свобода выбора узлов интегральной формулы в руках вычислителя, то можно ставить задачу следующим образом: выбрать узлы ивеса квадратурной формулы для обеспечения алгебраической точности187VII.
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕданной квадратуры как можно более высокого порядка. Иными словами,мы будем подбирать узлы квадратурной формулы xk и веса квадратурнойформулы ck, чтобы квадратурная формула с 2n параметрами11nf ( x)dx ci f ( xi )(3.1)i 1была точна для многочленов как можно более высокой степени.
Для указанного интервала интегрирования эта задача решается наиболее просто,в других случаях нужно делать линейную замену переменных.Теорема. Если в качестве узлов квадратурной формулы (38)взять нули полиномов Лежандра qn(x): qn(xi) = 0, i = 1, …, n, а в качествевесов интегралы от базисных функций многочленов Лагранжа li(x) степени n – 1, а именно1( x x1 )( x x2 )...( x xi 1 )( x xi 1 )...( x xn )ci dx ,(3.2)( xi x1 )( xi x2 )...( xi xi 1 )( xi xi 1 )...( xi xn )1то квадратурная формула (38) будет точна при подстановке в нее вместо f (x) любого многочлена степени не выше 2n – 1.Приведем узлы и веса квадратур для нескольких первых квадратурГаусса:n 2, x1 3 / 3, x2 3 / 3, c1 c2 1;n 3, x1 x3 3 / 5, x2 0, c1 c3 5 / 9, c2 8 / 9,n 4, x1 x4 0.861136, x2 x3 0.339981,c1 c4 0.347855, c2 c3 0.652145.Для вычисления интеграла по произвольному отрезку [a,b] системаортогональных многочленов на заданном отрезке получается из многочленов Лежандра после линейной замены:ba ba(3.3)t [1,1] x t , x [ a, b] .22Для оценки остаточного члена интегрирования по формулам Гауссаможно воспользоваться оценкой(b a) 2 n 1 (n !) 4 (2 n )rnГ f ( ), [a, b] .(3.4)(2n 1)[(2n)!]3Приведем вычисленные значения коэффициента перед производнойв оценке остаточного члена для нескольких первых квадратурных формул Гаусса при интегрировании по отрезку [–1,1]:188VII.3.1.
Квадратурные формулы Гаусса–Кристоффеля1 (4)1f ( ), r3Г f (6) ( ),13515750(3.5)213r4Г f (8) ( ), r5Г f (10) ( ),...34728751237732650Видим быстрое затухание ошибки при достаточной гладкости интегрируемой функции.r2Г VII.3.1. Квадратурные формулы Гаусса–КристоффеляВ некоторых приложениях возникает необходимость вычисленияbинтегралов с заданной весовой функцией p(x): f ( x) p( x)dx . При некоaтором наборе весовых функций могут оказаться полезными квадратурные формулы Гаусса–Кристоффеля видаbanf ( x) p( x)dx ci f ( xi ) .(3.6)i 1Теорема. Квадратурная формула (3.6) точна для многочленовстепени не выше 2n – 1, если ее узлами xi являются корни многочленаQn(x) из семейства многочленов, ортогональных на промежутке интегрирования (a, b) с весом p(x), а весовыми коэффициентами являютсячисла1( x x1 )( x x2 )...( x xi 1 )( x xi 1 )...( x xn )ci p( x)dx .(3.7)(xi x1 )( xi x2 )...( xi xi 1 )( xi xi 1 )...( xi xn )1VII.4.
Приемы вычисления несобственных интеграловРассмотрим сходящиеся интегралы следующих двух типов:b f ( x)dx, причемf ( x) при x a (первый тип);(4.1)a f ( x)dx, (второй тип).(4.2)aЗамечание. Второй интеграл, вообще говоря, может быть сведенк первому заменой переменной интегрирования: t = 1/x. Поэтому покабудем говорить об интегралах первого типа.189VII. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕОчевидно, непосредственное использование квадратурных формултрапеций и Симпсона для вычисления таких интегралов невозможно (таккак точка x = a является для этих формул узлом интегрирования).
Пометоду прямоугольников вычисления формально провести можно, норезультат будет сомнительным, так как оценка погрешности теряетсмысл (производные подынтегральной функции не ограничены).Продемонстрируем приемы, которые позволяют получать в подобных случаях надежные результаты, на примере интеграла1I cos x0xdx.а) Иногда подходящая замена переменной интегрирования позволяет вообще избавиться от особенности.В рассматриваемом примере после замены x = t2 получаем1I 2 cos t 2 dt ,0и интеграл вычисляется с требуемой точностью по любой из квадратурных формул.б) Та же цель (избавление от особенности) достигается иногда предварительным интегрированием по частям:1I 0cos xxdx 2 x cos x1012 x sin x dx.0Последний интеграл формально может быть вычислен стандартнымобразом, но оценка погрешности для любой квадратурной формулы будет иметь лишь первый порядок, так как при x = 0 не существует втораяпроизводная от подынтегральной функции.
Проводя еще раз интегрирование по частям, придем к интегралу от дважды непрерывно дифференцируемой функции, который с гарантированной точностью может бытьвычислен по формулам трапеций или прямоугольников.в) Если упомянутыми простыми средствами избавиться от особенности не удается, то прибегают к универсальному методу выделенияособенности. В рассматриваемом случае представим интеграл в видесуммы двух интегралов:I I1 I 2 , I1 0cos xxqdx, I 2 190cos xxdx.VII.4.
Приемы вычисления несобственных интеграловВторой интеграл особенности не содержит и вычисляется по любойквадратурной формуле. Вопрос о выборе величины обсуждается ниже.Первый интеграл с требуемой точностью вычисляем аналитически,используя представление подынтегральной функции в окрестности особой точки (x = 0) в виде отрезка ряда по степеням x, который получимпосле замены cos x соответствующим рядом Тейлора: 1 x 2 x 4 (1)m x 2 mI1 02!4!( 2m)!x 1 2 9 / 2 (1) m4! 9dx 2 1 2 5 / 2 2! 511 2m1/ 2 .(2m)! 2m1 2Важно, что подобное аналитическое представление в малой окрестности особой точки можно получить практически во всех конкретныхслучаях. Как это сделать – зависит от квалификации вычислителя.Допустим, что мы решили ограничиться в полученном представлении первыми m слагаемыми.
При этом для данного примера мы допускаем погрешность, которая не превосходит последнего приведенного взаписи для I1 слагаемого в силу того, что ряд для – знакопеременный.Следовательно, для выбора двух параметров ( и m) имеем следующийкритерий:11 2m1/ 2 (2m)! 2m1/ 22(4.3)(ε/2 отводится в качестве допустимого уровня погрешности при вычислении I2).Таким образом, один из параметров (m или ) можно задавать посвоему усмотрению, второй – определяется из неравенства (4.3). Приэтом нужно принять в расчет следующее соображение.Если δ << 1, то существенно ухудшается оценка погрешности длялюбой квадратурной формулы, которую мы предполагаем использоватьдля вычисления I2, так как в качестве коэффициента при hp (где p – порядок точности выбранной формулы) фигурирует максимальное на отрезке[δ, 1] значение p-й производной от подынтегральной функции, котороепри x = δ в рассматриваемом случае имеет порядок δ – (p + 0.5).Кроме того, при вычислении интеграла I2 придется вычислятьподынтегральную функцию f (x) от аргумента либо равного (для формул трапеций и Симпсона), либо очень близкого к (для формулы прямоугольников с центральной точкой).
Но значение f () при малом мо-191VII. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕжет быть настолько большим (в рассматриваемом случае f () ~ 1 /||), что абсолютная погрешность функции f () не позволит вычислитьинтеграл с требуемой точностью при заданной длине мантиссы и выбранном шаге интегрирования.Следовательно, целесообразно задать «не слишком малое» (например, = 0.1), а затем m найти из условия (4.3).Замечание. Разумеется, если поиск последовательных членовразложения подынтегральной функции затруднителен, то приходитсяограничиваться доступными членами.
В этом случае из условия типа(4.3) находится параметр .Аккуратное вычисление интеграла с особенностью может быть выполнено гораздо более экономичными средствами. Это достигается спомощью приема регуляризации (метод Канторовича), или выделенияособенности.
Поясним его в более общей ситуации. Пусть требуетсявычислить интеграл:1f (t )0tdt ,где f (t) – гладкая функция. Регуляризация состоит в том, что проделывается тождественное преобразование1011f (t )dt f (t ) (t ) t 1/2dt (t )t 1/2dt .t00Функция φ(t) выбирается такой, чтобы первый интеграл правой части не содержал особенности и при небольшом объеме вычислений достаточно точно определялся хотя бы по формуле Симпсона.
Второй интеграл особенность содержит, но вычисляется аналитически. В данномслучае цель будет достигнута, если в качестве φ(t) взять отрезок рядаТейлора f (t) в точке t = 0. Это приводит к вычислению1f (t ) f (0) tf '(0)0t110tdt f (0) 1dt f '(0) tdt.0В примере, с которого мы начинали f (t) = cos t, приходим к вычислению11001/21/2 (cos t 1)t dt t dt.192VII.5. Вычисление интегралов от быстроосциллирующих функцийВторое слагаемое есть 2, первое вычислим по формуле Симпсона: сначала с шагом 0.5, что даст значение 1.807967, а затем с шагом 0.25, чтодаст значение 1.808850.