Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
8.3.1.3. Составная формула трапеций. . i f x x i-1 x 1 1 0 0 1 1 1 ( ) 2 2 2 2 i i i n n n n i i h I f f f f f f I h f f h f 3 3 ' ' '' ' ' 2 2 1 1 ( ) 1 ( ) ( ) ( ) ( ) ( ) 12 12 12 n n сост тр i i h b a b a R f f f f h n n , ξ [a,b], ξ i [x i-1 ,x i ]. Порядок точности формулы равен 2. 54 2.4. Составная формула центральных прямоугольников. X f x X X i-1 i-1/2 i 1/ 2 1/ 2 3 / 2 1/ 2 1/ 2 1 ( ) i i n n i i I hf I h f f f h f Погрешность для простой формулы центральных прямоугольников получена ранее. Проведя преобразования, аналогичные сделанным в предыдущем параграфе, найдём ' ' 2 ( ) ( ) ( ) 24 сост ц пр b a R f f h , ξ [a,b]. Если число элементарных отрезков n = 2m чётное, можно избежать использования половинных узлов, для чего частичный интеграл следует вычислять на отрезках [х 2i-2 , x 2i ] i=1,…m. Средняя точка этого отрезка – x 2i-1 , а его длина равна 2h. В этом случае составная формула центральных прямоугольников приобретает вид 2 1 1 3 2 1 1 2 1 1 2 ; 1, 2 ( ) 2 i i m i n i i I hf i m I h f f f f h f Так как шаг h нов. =2h, 2 2 2 ( ) ( ) ( ) ( ) ''( ) ( ) ''( )(2 ) ( ) ''( ) 24 24 6 сост ц пр нов сост ц пр сост ц пр b a b a b a R f f h R f f h R f f h Порядок точности формулы равен 2. ξ [a,b]. 55 8.3.1.4. Составная формула Симпсона. 1 2 ( ) i i x i x I P x dx Здесь P 2 (x) – интерполяционный полином 2-ой степени, построенный на узлах x i-1 , x i-1/2 , x i 1 2 1 1/ 2 1 2 1/ 2 1/ 2 1/ 2 2 2 1 1/ 2 2 ( ) ( ) ( ) / 2 ( ) ( 4 ) 6 i i i i i i i i i i x i i i x f f f f f P x f x x x x h h h P x dx f f f h=x i -x i-1 Составная квадратурная формула Симпсона имеет вид 1 0 1/ 2 1 3/ 2 2 1 1/ 2 0 1/ 2 1 1 1 ( 4 2 4 2 ... 2 4 ) ( 4 2 ) 6 6 n n n i n n n n i i i i i h h I I f f f f f f f f f f f f . Погрешность формулы 5 (4) (4) 4 1 ( ) ( ) ( ) ( ) 2880 2880 n сост симпс i i h b a R f f f h , ξ [a,b]. Если число элементарных отрезков n = 2m чётное, можно избежать использования половинных узлов аналогично тому, как это описано в разделе «формула центральных прямоугольников». Составная квадратурная формула Симпсона приобретает вид: 1 0 2 1 2 1 1 ( 4 2 ) 3 m m n i i i i h I f f f f с погрешностью (4) 4 ( ) ( ) ( ) 180 сост симпс b a R f f h , ξ [a,b]. Порядок точности квадратурной формулы Симпсона равен 4-ём. 8.4. Квадратурные формулы наивысшей степени точности (формулы Гаусса). Степень точности квадратурной формулы может быть повышена путем специального способа выбора узлов. В квадратурных формулах Гаусса весовые коэффициенты и узлы 56 подбираются так, чтобы формула была точна для полинома наивысшей возможной степени N. Пусть весовая функция удовлетворяет условиям: | ( ) | , i 0, ( ) 0 b b i a a p x x dx p x dx .(*) Тогда параметры A k и x k квадратурной формулы с n узлами можно выбрать так, чтобы она была точна для полиномов N=2n-1 степени. Теорема 1. Пусть квадратурная формула имеет стандартный вид 1 ( ) ( ) ( ) (1) b n k k k a p x f x dx A f x Тогда, чтобы (1) имела степень точности (2n-1) необходимо и достаточно выполнение следующих условий: 1) Формула (1) должна быть интерполяционной, т. е. её весовые коэффициенты должны вычисляться по формулам ( ) ( ) ; ( ) '( ) b k k k a x A p x dx x x x 2) Узлы интерполирования x k должны быть такими, чтобы многочлен ( ) ( ) ( ) b a p x x Q x dx , ω(x)=(x-x 1 )…(x-x n ) был ортогонален с весом p(x) ко всякому многочлену Q(x) степени, меньшей n, т. е. чтобы ( ) ( ) ( ) 0 b a p x x Q x dx Доказательство. Необходимость. Предположим, что (1) точна для всех полиномов степени (2n-1). Тогда она должна быть точна и для полиномов степени (n-1), а, следовательно, согласно доказанной ранее теореме о степени точности ИКФ эта формула должна быть интерполяционной. 57 Предположим, что Q(x) = Q n-1 (x) любой полином степени (n-1). Тогда ω(x)Q(x) есть полином степени (2n-1). Пусть f(x) = ω(x)Q(x) и, следовательно, 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 ( ) ( ) ( ) 0, b b n k k k k a a b a p x f x dx p x x Q x dx A x Q x p x x Q x dx что соответствует условию (2) теоремы. Необходимость доказана. Достаточность. Рассмотрим любой многочлен f(x) степени (2n-1). Разделим его на ω(x). Тогда f(x) = ω(x)Q i (x)+r i (x), где Q i (x) и r i (x) многочлены степени i (n-1). Поэтому 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b b b n i i k k i k k a a a p x f x dx p x x Q x dx p x r x dx A x r x Причем, последнее равенство выполнено точно, т. к. первый интеграл равен нулю по условию теоремы, а для второго интеграла ИКФ является точной. Достаточность доказана. Свойства многочлена ω(x) уточняются дополнительными теоремами. Теорема 2. Пусть выполнены условия (*). Тогда существует единственный многочлен ω(x) указанного выше вида, ортогональный с весом p(x) ко всякому многочлену степени меньше n. Теорема 3. Пусть p(x) > 0 на отрезке [a,b], а многочлен ω(x) ортогонален с весом p(x) любому многочлену Q n-1 (x). Тогда корни ω(x) лежат внутри интервала [a,b] и различны между собой. Заметим, что требование нахождения корней ω(x) внутри интервала [a,b] становится необходимым, если вне этого интервала f(x) не определена. 58 Возникает вопрос об отыскании многочленов ω(x), удовлетворяющих условиям, сформулированным в теоремах. Рассмотрим случай p(x)=1, x [a,b]. Используя замену переменной x=(a+b)/2+t(b-a)/2, задачу интегрирования преобразуют к виду 1 1 ( ) ; [ 1,1]. 2 2 2 b a b a a b b a f x dx f t dt t Известно, что ортогональную на интервале [-1,1] систему многочленов с постоянным весом образуют многочлены Лежандра. Для их вычисления удобна рекурентная формула (n+1) ℒ n+1 (t)+(2n+1)t ℒ n +n ℒ n-1 (t) = 0 , ℒ 0 =1, ℒ 1 =t. Другой возможный способ - вычисление по формуле 2 n 1 ( 1) ( ) 2 ! n n n n d t t n dt Полиномы Лежандра 1) ортогональны на интервале [-1,1] любому многочлену Q i (t) степени, меньшей n: 1 1 ( ) ( ) 0, i n Q t t dt 2) корни t k действительные, простые и принадлежат интервалу [-1,1]. Полиномы Лежандра удовлетворяют условиям, сформулированным в теоремах. Таким образом, узлы КФ Гаусса могут совпадать с корнями полиномов Лежандра. Веса A k вычисляются по формуле 2 2 n 2 , k=1,...,n. (1 )[ '( )] k k k A t t Заметим, что ортогональность сохраняется при умножении функции или вектора на скалярный коэффициент. Поэтому в 59 качестве ω(t) можно использовать нормированный полином Лежандра ℒ n /a n , где a n коэффициент при старшей степени полинома Лежандра. Погрешность ИКФ Гаусса определяется по формуле 2 1 (2 ) 2 2 1 2 ( ) ( ) ( ), [a,b] 2 2 ( !) (n)= (2 1)(2 )! (2 !) n n n n b a R f n f n n n n Напомним, что под n подразумевается число узлов ИКФ. Величина α(n) быстро уменьшается сростом n, что видно из таблицы N 2 3 4 5 6 α(n) 1 135 1 15750 1 3472875 9 1 1.24 *10 9 1 649 *10 Значения весов и узлов КФ Гаусса на интервале [-1,1] не зависят от вида подынтегральной функции. В настоящее время построены таблицы узлов и весов квадратур Гаусса по крайней мере до n=4096 c 20 двадцатью десятичными знаками. Если p(x)= 2 1 ( 1 ) x , то за узлы x k следует принимать корни полинома Чебышева. В этом случае 1 2 1 1 (2 ) 2 1 1 2 1 ( ) ( ), x cos 2 1 ( ) ( ), [-1,1]. 2 (2 )! n k k k n n n k f x dx f x n n x R f f n В приведенных формулах x [-1,1]. Произвольный отрезок может быть преобразован к стандартному соответствующей заменой переменных. 60 8.5. Практическая оценка погрешности (правило Рунге). Пусть I h приближенное значение интеграла I, вычисленное с помощью какой либо составной квадратурной формулы. Предположим, что для погрешности этой формулы справедливо представление I - I h = Ch k + o(h k ), (1) где C ≠ 0 и k > 0 величины, не зависящие от h. Запись o(h k ) означает, что o(h k ) пренебрежимо мало по сравнению с h k при h ⟶0. Тогда величина Ch k называется главным членом погрешности квадратурной формулы. Заметим, что из (1) следует |I –I h | h k с некоторой постоянной > |C|. Поэтому число k является ни чем иным как порядком точности соответствующей квадратурной формулы. Если подынтегральная функция достаточно гладкая, то для каждой из составных квадратурных формул существует главный член погрешности. В ходе доказательства этого утверждения разработана методика выделения главного члена погрешности. Процедуру выделения главного члена погрешности продемонстрируем на примере составной квадратурной формулы трапеций. Ранее было получено 3 2 2 ' ' ' ' , 0 1 1 ( ) ( ) ( ) ''( ) 12 12 12 b n n сост тр i i n h a h h h R f f hf f x dx Таким образом, R сост.тр. = Ch 2 + o(h 2 ), где 1 ''( ) 12 b a C f x dx При половинном шаге из (1) следует I - I h/2 = (1/2 k )Ch k + o((h/2) k ) (2) 61 При достаточно малом h бесконечно малыми в(1) и (2) можно пренебречь. В этом случае, вычитая из (1) (2), получим I h/2 - I h ≈ (1/2 k )Ch k (2 k -1) ≈ (I-I h )(2 k -1)(1/2 k ) ≈( I - I h/2 ) (2 k -1). I-I h/2 ≅ /2 2 1 h h k I I (3) I - I h ≅ /2 2 2 1 h h k k I I (4) k =1 для формул левых и правых прямоугольников; k = 2 для формул центральных прямоугольников и трапеций; k = 4 для формулы Симпсона. 8.6. Адаптивные процедуры численного интегрирования. Пользователь такой процедуры задает отрезок [a,b], правило вычисления подынтегральной функции и требуемую точность ε > 0. Программа стремится, используя по возможности минимальное число узлов интегрирования, распределить их так, чтобы найденное значение удовлетворяло неравенству 1 | ( ) | i b n h i i a I I f x dx . Здесь i h i I приближенное значение интеграла I i на отрезке [x i-1 ,x i ], вычисляемое по выбранной квадратурной формуле. Программа подбирает длину очередного отрезка так, чтобы выполнялось неравенство 1 n i i I . Это условие будет выполнено, если погрешность на очередном отрезке i i h b a . Действительно, в этом случае 1 1 n n i i i i h b a 62 Алгоритм 1.Начальная установка ε, a1=a, b1=b, h=b1-a1, s=0. 2.Используя выбранную составную квадратурную формулу, рассчитывают значения интеграла I h и I h/2 , оценивают погрешность ε i по правилу Рунге. 3. Если i i h b a , полагают h=h/2,b1=a1+h. Новый отрезок есть [a1,b1]. Возвращаются к п.2. 4.Если условие п.3 не выполнено, то s=s+I h/2 , a1=b1, b1=b, h=b- a1 (h длина очередного отрезка). 5.Если h > 0, возвращаются к п.2. Иначе s есть искомый результат. 8.7. Обусловленность квадратурных формул интерполяционного типа. Значения подынтегральной функции всегда вычисляются с погрешностью. Предположим, что |f(x) – f*(x)| △(f*) для всех x [a,b]. Тогда | ( ) * ( ) | ( ) (f *) b b a a f x dx f x dx b a Для квадратурной формулы имеем 1 1 1 | ( ) * ( ) | | | (f *) n n n i i i i i i i i A f x A f x A Таким образом, абсолютное число обусловленности квадратурной формулы ν = 1 | | n i i A 63 Все ИКФ точны для многочленов нулевой степени и поэтому b- a= 1 1* b n i i a dx A Если все веса ИКФ положительны, то ν = 1 | | n i i A = 1 n i i A =b-a. Если среди весов имеются отрицательные, то 1 | | n i i A > 1 n i i A =b-a. Например, для формул Ньютона -Котеса при n = 20, ν = 8.3(b-a), a при n = 30, ν = 560(b-a). Весовые коэффициенты формулы Гаусса все положительны при любом n. |