WWW.DISSERS.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА

   Добро пожаловать!


Pages:     || 2 | 3 | 4 |
ВОРОНЕЖСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ Факультет прикладной математики, информатики и механики ГУДОВИЧ Н.Н.

ИЗБРАННЫЕ ВОПРОСЫ КУРСА ЧИСЛЕННЫХ МЕТОДОВ Выпуск 4. Численное интегрирование.

ВОРОНЕЖ 2002 2 10. Интерполяционные квадратурные формулы.

Пусть f – непрерывная на отрезке [a,b] функция: f C[a,b].

Рассмотрим определённый интеграл b (I f ) = f ( x )dx.

( 1.1) a Если f 0, то величина (1.1) есть площадь криволинейной трапеции, ограниченной сверху графиком функции f ( рис. 1.1).

Напомним, что под квадратурой плоской фигуры в узком смысле понимают задачу построения квадрата, площадь которого равна площади этой фигуры, а под квадратурой фигуры в широком смысле – просто задачу нахождения площади фигуры. Поэтому формулы для приближённого вычисления интеграла (1.1) называют квадратурными формулами.

Зачем нужны формулы для приближенного вычисления интеграла (1.1) Дело в том, что воспользоваться формулой Ньютона-Лейбница b (f x )dx F( b ) -= F a( ) ( 1.2 ) a невсегда возможно.

Во-первых, первообразная F(x) для заданной подинтегральной функции f(x) может оказаться неизвестной для вычислителя.

Во-вторых, эта первообразная может оказаться неэлементарной функцией, т.е. не принадлежать классам степенных функций, многочленов, рациональных функций, тригонометрических и обратных тригонометрических функций, показательных и логарифмических функций, гиперболических и обратных гиперболических функций, а также функций, получаемых из перечисленных с 3 помощью четырёх арифметических действий и суперпозиций, взятых в конечном числе. В этом случае нахождение значений первообразной F(a), F(b), фигурирующих в формуле Ньютона-Лейбница (1.2), может оказаться сложной задачей.

Отметим попутно, что первообразная F может оказаться неэлементарной и в случае, когда сама функция f – элементарна. Классический пример такого рода даёт функция (f x ) = e- x2, элементарная как суперпозиция показательной и степенной функций. Её первообразная – встречающийся в теории вероятностей “ интеграл ошибок “ erf(x) – неэлементарная функция.

В-третьих, если функция f задана не аналитически, а таблицей своих значений в некоторых точках отрезка [a,b], то найти её первообразную и воспользоваться формулой (1.2) непредставляется возможным.

Во всех этих случаях интеграл (1.1) приходится вычислять приближённо, заменяя подинтегральную функцию f близкой функцией, интеграл от которой уже может быть вычислен по формуле Ньютона-Лейбница; если в качестве такой приближающей функции берётся интерполяционный многочлен, то формулу для приближённого значения интеграла называют интерполяционной квадратурной формулой.

Итак, пусть на отрезке [a,b] выбраны попарно различные точки x0,x1,...,xn, (1.3) и пусть в этих точках заданы значения функции f :

f(x0 ), f(x1 ),..., f(xn ). (1.4) Имеем:

f(x) = pn (x) + rn (x), (1.5) где pn (x) = pn (x; {xi }; f ) - интерполяционный многочлен, построенный по значениям (1.4) функции f в узлах (1.3), а rn (x) = rn (x; {xi }; f ) - погрешность интерполяции в точке x. Подстановка (1.5) в (1.1) даёт b b b (f x )dx pn( x )dx += n r ( x )dx.

a a a Первое слагаемое в правой части полученного равенства b (I f ) = pnn ( x )dx (1 6. ) a есть точное значение определённого интеграла от интерполяционного многочлена, вычисляемое по формуле Ньютона-Лейбница и принимаемое в качестве приближённого значения определённого интеграла от исходной функции f ; величина же b (R f ) = ( x )dx (1 7. ) nn r a есть погрешность этого приближённого значения.

Установим структуру выражения (1.6). Записывая интерполяционный многочлен pn (x) в форме Лагранжа и подставляя это представление в (1.6), получим:

n x( - x ) j =0j b b n n kj (I f ) pnn ( x )dx == dx) f ( x )= Akk, ( 0k f ( xk ) n == 0k a a x( - x ) jk =0j kj гдеиспользованы обозначения n x( - x ) j =0j b kj Ak = dx. (1 8. ) n a x( - x ) jk =0j kj Следовательно, приближенное значение In ( f ) искомого интеграла I ( f ) есть линейная комбинация n (I f ) f ( xkn )= Ak ( 1.9 ) k =значений функции f в точках (1.1) с коэффициентами (1.8).

Определение 1.1. Сумму в правой части равенства (1.9) называют квадратурной суммой, а фигурирующие в ней точки xk и коэффициенты Ak - соответственно квадратурными узлами и квадратурными коэффициентами.

Замечание 1.2. Квадратурные коэффициенты (1.8) представляют собой определённые интегралы от многочленов. Так как нахождение первообразной от многочлена не представляет труда, эти коэффициенты легко могут быть вычислены по правилу Ньютона-Лейбница. Отметим, что квадратурные коэффициенты, как это видно из формулы (1.8), независят от значений функции f в узлах интерполяции, а определяются исключительно расположением этих узлов на отрезке [a,b]. Следовательно, их достаточно вычислить лишь один раз, а затем использовать для вычисления определённых интегралов от всех функций f, значения которых в этом набореузлов заданы.

Замечание 1.3. Квадратурные коэффициенты полезно нормировать, воспользовавшись линейной заменой переменного x=a + (b – a)t, (1.10) порождающей взаимно однозначное отображение отрезка [0,1] оси t на отрезок [a,b] оси x. Обозначая через tm прообразы узлов xm, подставляя в (1.8) вытекающие из (1.10) равенства x – xj = ( b – a )( t – tj ), xk – xj = ( b – a )( tk – tj ), dx = ( b – a )dt и пересчитывая пределы интегрирования, получим n t( - t ) j =0j kj (A b -= a ) dt, k = 0,1,..., n.

k n t( - t ) jk =0j kj Следовательно, квадратурную сумму (1.9) можно представить в виде n (I f ) ( b -= a ) f ( xkn ) Bk, (1.11) k =где Bk - нормированные квадратурные коэффициенты:



n t( - t ) j =0j kj Bk = dt, k = 0,1,..., n. (1.12 ) n t( - t ) jk =0j kj Укажем наиболее употребительные интерполяционные квадратурные формулы.

При n=0 подинтегральная функция f заменяется интерполяционным многочленом нулевой степени, т.е. константой p0 (x) = f(x0 ).

Геометрически это соответствует тому (см. рис.1.2 ), что площадь криволинейной трапеции, ограниченной сверху графиком функции f, заменяется площадью прямоугольника с длинами сторон b – a, f(x0 ). Поэтому полученная приближённая формула b (f x )dx ( b - a ) f ( x0 ) ( 1.13 ) a называется формулой прямоугольника.

Замечание 1.4. В качестве квадратурного узла x0 обычно берут левый конец отрезка [a,b], правый конец этого отрезка или его середину.

Соответственно ( см. рис. 1.3 ) получают формулу левого, правого и центрального прямоугольника:

- a ) f ( a ) b( b (f x )dx b( - a ) f ( b ) (. 1.14 ) a + ba - a ) f ( ) b( При n=1 подинтегральная функция f заменяется интерполяционным многочленом первой степени, т.е. линейной функцией. Если при этом в качестве узлов интерполяции взяты концы отрезка ( x0 = a, x1 =b ), то геометрически это соответствует тому ( рис. 1.4 ), что в качестве приближения к площади криволинейной трапеции принимается площадь прямолинейной трапеции с основаниями f(a), f(b) и высотой b – a. А так как площадь прямолинейной трапеции равна произведению полусуммы оснований на высоту, приходим к формуле:

b (f a ) + f ( b ) 1 (f x )dx ( b - a ) b( -= a )( f ( a ) + f ( b )), ( 1.15 ) 2 a которую естественно назвать формулой трапеции; квадратурные коэффициенты в этом случае имеют значения B0 = B1 = 1/2.

Наконец, если n=2, то функция f заменяется многочленом второй степени, графиком которой является парабола; если при этом в качестве квадратурных узлов приняты концы отрезка и его середина ( x0 = a, x1 = (a+b)/2, x2 = b ), то получается квадратурная формула, называемая формулой параболы, или формулой Симпсона ( рис. 1.5 ).

Вывести формулу параболы из геометрических соображений, как это делалось для формул прямоугольника и трапеции, затруднительно, поэтому нам придётся воспользоваться формулами (1.11), (1.12). Вычисления по формулам (1.12) с учётом равенств t0 =0, t1 =1/2, t2 =1 дают:

1 1 t( )( t -- 1) t( t )( t -- t21 ) 3 1 B0 = dt = dt 2 t2 -= ( 2 t + 2 dt) = 6, t( t10 )( t0 -- t2 ) 0 0 ( )( -- 1) 1 1 t( t )( t -- t20 ) (t t - 1) B1 = dt = dt -= 4 ( t2 - t )dt =, 1 t( t01 )( t1 -- t2 ) 0 0 ( - ) 2 1 1 - ) (t t t( t )( t -- t10 ) 1 B2 = dt = dt 2 ( t2 -= )t dt = ;

t( t02 )( t2 -- t1 ) 2 0 0 (1 ) следовательно, формула Симпсона имеет вид:

b - ab (f x )dx f( ( a ) 4 f (( a ++ b ) / 2 ) + f ( b )). (1.16 ) a Обратимся теперь к вопросу о погрешности квадратурной формулы.

Вспоминая формулу для погрешности интерполяционного многочлена n( +1 ) n (f ( x )) (r x ) = x( - x ) n j n( + 1)! =0j и подставляя это выражение в (1.7), получим b n n( +1 ) (R f ) = (f ( x )) x - x )dx, (1.17 ) n ( j n( + 1)! =0j a где (x) – некоторая точка отрезка [a,b]. Так как модуль интеграла не превосходит интеграла от модуля, а модуль произведения равен произведению модулей сомножителей, имеем:

b n n( +1 ) (R f ) (f ( x )) ( x - x ) dx.

n j n( + 1)! =0j a Заменяя здесь модуль (n+1) – вой производной максимальным по отрезку [a,b] значением и замечая, что расстояние ( x – xj ) между точками x, xj отрезка [a,b] непревышает длины b – a отрезка, приходим к выводу о том, что для любой функции f класса C n+1[a,b] справедлива оценка n( 1 ) (R f ) max f ( x ) ( b - a )n++. (1.18 ) n n( + 1)! xa b В заключение коснёмся вопроса о сходимости при n (т.е. при неограниченном увеличении числа квадратурных узлов ) приближённого значения интеграла In ( f ) к его точному значению I( f ), т.е. вопроса о сходимости глобально-интерполяционного квадратурного процесса на классе непрерывных на отрезке [a,b] функций.

При исследовании глобальной интерполяции ставился вопрос о том, существует ли алгоритм задания узлов интерполяции n {n xi )n( }, =0i при котором последовательность интерполяционных многочленов pn ({xi(n)};f ), (n=0,1,2,... ), отвечающих произвольной непрерывной на отрезке [a,b] функции f, сходится при n к исходной функции f равномерно на этом отрезке [a,b], т.е. сходится к f в метрике пространства C[a,b] :

pf ({ xi )n( }, f ) - 0 при n для Cf [ a,b].

n [C a,b ] Ответ на этот вопрос отрицательный: при любом способе задания узлов найдётся ( для каждого способа, вообще говоря, своя ) функция f, для которой такая сходимость места неимеет. А так как рассматриваемый в данном пункте способ приближённого вычисления определённых интегралов основан на глобальной интерполяции, отмеченный только что факт расходимости глобальных интерполянтов наводит на мысль о возможных осложнениях и в задаче вычисления интегралов. Такие сложности действительно имеют место, хотя и в более слабой форме: в отличие от процесса глобальной интерполяции, где все способы задания узлов, условно говоря, «плохие», в случае глобальноинтерполяционного квадратурного процесса имеются и «хорошие» способы задания узлов, при использовании которых приближённые значения интеграла In (f) сходятся к точному значению I(f) для любой непрерывной функции f :

In ( f ) I( f ) при n для f C[a,b]. (1.19) Для того, чтобы соотношение (1.19) имело место, квадратурные коэффициенты, а они, как отмечено выше, определяются исключительно способом задания квадратурных узлов, должны удовлетворять некоторому условию; это условие мы сейчас и сформулируем.





Заметим,что выражение (1.1) можно рассматривать как значение на элементе f C[a,b] оператора I, сопоставляющего функции f значение определённого интеграла от этой функции по отрезку [a,b]. В силу известных свойств определённого интеграла ( интеграл от суммы двух функций равен сумме интегралов от слагаемых, постоянный множитель можно выносить за знак интеграла ) этот оператор – линейный. Аналогично, является линейным и оператор In, сопоставляющий функции f значение квадратурной суммы (1.9).

Напомним, что операторы, областью значений которых является пространство вещественных чисел R, называются функционалами, так что I,In - линейные функционалы на пространстве непрерывных функций C[a,b].

Поскольку нормой вещественного числа как элемента пространства R является его абсолютная величина, для функционалов I, In справедливы неравенства b b (f x )dx (f x ) dx (I f ) a a supI == sup sup f max f ( x ) max f ( x ) Cf [ a,b ] f f [C a,b ] xa b xa b b b max f ( x ) dx max f ( x ) dx xa b xa b a a sup = sup sup ( b -= a ) = b - a, max f ( x ) max f ( x ) f f f xa b xa b n n (f xk )n( ) Ak )n( (f xk )n( ) Ak )n( (I f ) n =0k =0k supI == sup sup n f max f ( x ) max f ( x ) f f f xa b xa b n max f ( x ) Ak )n( n n xa b =0k sup sup ( Akn( ) ) == Ak )n(.

max f ( x ) f f 0k == 0k xa b Из этих неравенств следует, что I, In – ограниченные функционалы, причём n bI - a, In Ak )n(. ( 1.20 ) =0k Заметим, что на самом деле здесь имеют место равенства.

Действительно, полагая f0 (x) 1, получим неравенство (I f ) (I f0 ) - ab supI = = -= ab, f f0 f противоположное левому из неравенств (1.20).

Аналогично, задавая функцию f1 в квадратурных узлах формулой (f xk )n( ) sign Ak )n(, k == 0,1,...,n и доопределяя её между соседними узлами по линейности, а между крайними узлами и концами отрезка константами, равными значениям f1 в крайних узлах, получим непрерывную на [a,b] функцию с нормой, равной единице. С использованием этой функции получаем неравенство n sign( Ak )n( ) Ak )n( ) n (I f ) (I f1 ) nn =0k )n( supI = = =,A n k f f1 f =0k противоположное правому из неравенств (1.20).

Следовательно, n )n( = AI. ( 1.21) n k =0k Отметим, что мы добавили в обозначения квадратурных узлов и квадратурных коэффициентов верхний индекс n, чтобы подчеркнуть тот факт, что при каждом n они представляют собой свои наборы точек и чисел.

Теорема 1.5. Для сходимости глобально-интерполяционного квадратурного процесса на классе C[a,b], т.е. для справедливости соотношения (1.19), необходимо и достаточно, чтобы величины (1.21) были равномерно по n ограничены:

n )n( CA <, C не зависит от.n (1.22 ) k k =Доказательство. На языке функционального анализа соотношение (1.19) означает сильную сходимость функционалов In к функционалу I, т.е.

сходимость значений In (f) функционалов In к значению I(f) функционала I на произвольном элементе пространства C[a,b]. Согласно теореме БанахаШтейнгауза для такой сходимости необходимо и достаточно одновременное выполнение двух условий:

1. Сильная сходимость In к I на каком-либо подмножестве пространства C[a,b], плотном в этом пространстве.

2. Равномерная по n ограниченность норм функционалов In.

Рассмотрим подмножество пространства C[a,b] – совокупность всех многочленов. Согласно теореме Вейерштрасса из математического анализа любая непрерывная на отрезке [a,b] функция может быть с любой точностью ( в метрике пространства C[a,b] ) приближена многочленом, а это и означает плотность множества многочленов в пространстве C[a,b]. С другой стороны, поскольку для многочлена f степени m при n m интерполяционные многочлены pn (f) совпадают с самим многочленом f, при этих n интерполяционная квадратурная формула даёт точное значение интеграла. По этой причинедля любого многочлена f величины In ( f ) сходятся при n к величине I ( f ), и первое из перечисленных выше условий теоремы БанахаШтейнгауза оказывается выполненым.

Что же касается второго условия, то в силу (1.21) оно в точности совпадает с условием (1.22) ; последнее условие, таким образом, и есть необходимое и достаточное условие сходимости глобально-интерполяционного квадратурного процесса на классе C[a,b].

Замечание 1.6. На практике часто встречается случай, когда подинтегральная функция задаётся таблицей своих значений в равноотстоящих узлах )n( ax += ih, i = 0,1,...,n, h = ( b - a ) / n ;

i интерполяционные квадратурные формулы с таким способом задания узлов называют формулами Ньютона-Котеса. Примером такой формулы служит формула Симпсона (1.16), отвечающая случаю n=2 ; достаточно часто используется также формула «трёх восьмых», отвечающая n=3. Формулы Ньютона-Котеса с большими значениями n не применяются, поскольку этот квадратурный процесс не удовлетворяет условию сходимости (1.22) : в двумерном массиве квадратурных коэффициентов Ak(n) по мереувеличения n начинают встречаться коэффициенты со сколь угодно большими модулями.

Повышение точности вычисления интеграла в случае равноотстоящих узлов достигается не за счёт использования глобальных интерполянтов с высокой степенью n, а за счёт применения локальной интерполяции с малым n и большим числом частичных отрезков разбиения.

20. Локально-интерполяционные квадратурные формулы.

Pages:     || 2 | 3 | 4 |










© 2011 www.dissers.ru - «Бесплатная электронная библиотека»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.