Вычисление интеграла по методу симпсона. Оценка точности вычисления «неберущихся» интегралов Вывод Формулы Симпсона

Метод трапеций

Разобьем отрезок на равных частей при помощи точек:

Метод трапеций заключается в замене интеграла суммой:


Абсолютная погрешность приближения, полученного по формуле трапеций, оценивается с помощью формулы, где.

Метод парабол (метод Симпсона)

а) Через любые три точки с координатами проходит только одна парабола.

б) Выразим площадь под параболой на отрезке через:

Учитывая значения и из пункта а) следует:

в) Разобьем отрезок на равных частей при помощи точек:

Метод парабол заключается в замене интеграла суммой:

Для приближенных практических расчетов применяется формула:


Абсолютная погрешность вычисления по формуле (4) оценивается соотношением, где.

Оценка точности вычисления «неберущихся» интегралов

В данной работе вычисление абсолютной и относительной погрешности проводится при условии, что известно точное значение определенного интеграла. Однако не всякая первообразная, даже тогда, когда она существует, выражается в конечном виде через элементарные функции. Таковы первообразные, выраженные интегралами, и т.д. Во всех подобных случаях первообразная представляет собой некоторую новую функцию, которая не сводится к комбинации конечного числа элементарных функций.

Определенные интегралы от таких функций можно вычислить только приближенно. Для оценки точности вычисления в таких случаях используют, например, правило Рунге. В данном случае интеграл вычисляется по выбранной формуле (прямоугольников, трапеций, парабол Симпсона) при числе шагов, равном n, а затем при числе шагов, равном. Погрешность вычисления значения интеграла при числе шагов, равном, вычисляется по формуле Рунге:, для формул прямоугольников и трапеций, а для формулы Сипсона. Таким образом, интеграл вычисляется для последовательных значений числа шагов, ..., где - начальное число шагов. Процесс вычислений заканчивается, когда для очередного значения будет выполнено условие, где - заданная точность.

Для того чтобы не вычислять один и тот же интеграл по нескольку раз для разных разбиений отрезка интегрирования, можно вычислить шаг интегрирования заранее.

Пример. Выбрать шаг интегрирования для вычисления интеграла с точностью 0,01 пользуясь квадратурными формулами прямоугольников, трапеций, Симпсона.

Квадратурная формула прямоугольников.

Вычислим, при каком шаге погрешность будет составлять 0,01:

подынтегральный трапеция парабола неберущийся

Поскольку, то.

При шаге отрезок разбивается на равностоящих узлов.

Квадратурная формула трапеций.

Поскольку, .

При шаге,отрезок разбивается на равностоящих узлов.

Квадратурная формула Симпсона.

Вычислим, при каком шаге погрешность составит 0,01:

При шаге, отрезок разбивается на равностоящих узлов.

Как и следовало ожидать, наименьшее количество равностоящих узлов получается при вычислении интеграла по квадратурной формуле Симпсона.

Студенту предлагается работа, состоящая из четырех этапов:

  • 1 этап - точное вычисление определенного интеграла.
  • 2 этап - приближенное вычисление определенного интеграла одним из методов: прямоугольников или трапеций.
  • 3 этап - приближенное вычисление определенного интеграла методом парабол.

4 этап - расчет и сравнение абсолютной и относительной ошибок приближенных методов: , где - точное решение интеграла, - значение интеграла, полученное с помощью приближенных методов.

Построение графика подынтегральной функции.

Варианты и образец выполнения РГР приведены ниже.

Варианты

№ варианта

Образец выполнения РГР

Задание. Вычислить интеграл

1. Точное вычисление:


2. Приближенное вычисление с помощью формул прямоугольников:

Составим таблицу:

По первой формуле прямоугольников получаем:

0,1 = 0,1·3,062514 = 0,306251.

По второй формуле прямоугольников получаем:

0,1 = 0,1· 4,802669 = 0,480267.

В данном случае первая формула дает значение интеграла с недостатком, вторая - с избытком.

3. Приближенное вычисление по формуле трапеций:

В нашем случае получаем:

0,1 = =0,1 = 0,1·4,095562 = =0,409556.


Вычислим относительную и абсолютную погрешности.

4. Приближенное вычисление по формуле Симпсона:

В нашем случае получаем:


Вычислим относительную и абсолютную погрешности.

В действительности, = 0,40631714.

Таким образом, при разбиении отрезка на 10 частей по формуле Симпсона мы получили 5 верных знаков; по формуле трапеций - три верных знака; по формуле прямоугольников мы можем ручаться только за первый знак.

Кафедра «Высшей математики»

Выполнил: Матвеев Ф.И.

Проверила: Бурлова Л.В.

Улан-Удэ.2002

1.Численные методы интегрирования

2.Вывод формулы Симпсона

3.Геометрическая иллюстрация

4.Выбор шага интегрирования

5.Примеры

1. Численные методы интегрирования

Задача численного интегрирования заключается в вычислении интеграла

посредством ряда значений подынтегральной функции .

Задачи численного интегрирования приходится решать для функций, заданных таблично, функцией, интегралы от которых не берутся в элементарных функциях, и т.д. Рассмотрим только функции одной переменной.

Вместо функции, которую требуется проинтегрировать, проинтегрируем интерполяционный многочлен. Методы, основанные на замене подынтегральной функции интерполяционным многочленом, позволяют по параметрам многочлена оценить точность результата или же по заданной точности подобрать эти параметры.

Численные методы условно можно сгруппировать по способу аппроксимации подынтегральной функции.

Методы Ньютона-Котеса основаны на аппроксимации функции

полиномом степени . Алгоритм этого класса отличается только степенью полинома. Как правило, узлы аппроксимирующего полинома – равноотносящие.

Методы сплайн-интегрирования базируются на аппроксимации функции

сплайном-кусочным полиномом.

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

Методы Монте-Карло используются чаще всего при вычислении кратных интегралов, узлы выбираются случайным образом, ответ носит вероятностный характер.


суммарная погрешность погрешность усечения

погрешность округления

Независимо от выбранного метода в процессе численного интегрирования необходимо вычислить приближенное значение интеграла и оценить погрешность. Погрешность уменьшается при увеличении n-количества

разбиений отрезка

. Однако при этом возрастает погрешность округления

за счет суммирования значений интегралов, вычисленных на частичных отрезках.

Погрешность усечения зависит от свойств подынтегральной функции и длины

частичного отрезка.

2. Вывод формулы Симпсона

Если для каждой пары отрезков

построить многочлен второй степени, затем проинтегрировать его и воспользоваться свойством аддитивности интеграла, то получим формулу Симпсона. Рассмотрим подынтегральную функцию на отрезке . Заменим эту подынтегральную функцию интерполяционным многочленом Лагранжа второй степени, совпадающим с в точках :

Проинтегрируем

:

и называется формулой Симпсона.

Полученное для интеграла

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

Оценим теперь погрешность интегрирования по формуле Симпсона. Будем считать, что у

на отрезке существуют непрерывные производные . Составим разность

К каждому из этих двух интегралов уже можно применить теорему о среднем, поскольку

непрерывна на и функция неотрицательна на первом интервале интегрирования и неположительна на втором (то есть не меняет знака на каждом из этих интервалов). Поэтому:

(мы воспользовались теоремой о среднем, поскольку

- непрерывная функция; ).

Дифференцируя

дважды и применяя затем теорему о среднем, получим для другое выражение: , где

Из обеих оценок для

следует, что формула Симпсона является точной для многочленов степени не выше третьей. Запишем формулу Симпсона, напрмер, в виде: , .

Если отрезок

интегрирования слишком велик, то его разбивают на равных частей (полагая ), после чего к каждой паре соседних отрезков , ,..., применяют формулу Симпсона, именно:

Запишем формулу Симпсона в общем виде.

Для нахождения определенного интеграла методом трапеций площадь криволинейной трапеции также разбивается на n прямоугольных трапеций с высотами h и основаниями у 1 , у 2 , у 3 ,..у n , где n - номер прямоугольной трапеции. Интеграл будет численно равен сумме площадей прямоугольных трапеций (рисунок 4).

Рис. 4

n - количество разбиений

Погрешность формулы трапеций оценивается числом

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

Формула Симпсона

Если для каждой пары отрезков построить многочлен второй степени, затем проинтегрировать его на отрезке и воспользоваться свойством аддитивности интеграла, то получим формулу Симпсона.

В методе Симпсона для вычисления определенного интеграла весь интервал интегрирования разбивается на подинтервалы равной длины h=(b-a)/n. Число отрезков разбиения является четным числом. Затем на каждой паре соседних подинтервалов подинтегральная функция f(x) заменяется многочленом Лагранжа второй степени (рисунок 5).

Рис. 5 Функция y=f(x) на отрезке заменяется многочленом 2-го порядка

Рассмотрим подынтегральную функцию на отрезке. Заменим эту подынтегральную функцию интерполяционным многочленом Лагранжа второй степени, совпадающим с y= в точках:

Проинтегрируем на отрезке.:

Введем замену переменных:

Учитывая формулы замены,


Выполнив интегрирование, получим формулу Симпсона:

Полученное для интеграла значение совпадает с площадью криволинейной трапеции, ограниченной осью, прямыми, и параболой, проходящей через точки На отрезке формула Симпсона будет иметь вид:

В формуле параболы значение функции f(x) в нечетных точках разбиения х 1 , х 3 , ..., х 2n-1 имеет коэффициент 4, в четных точках х 2 , х 4 , ..., х 2n-2 - коэффициент 2 и в двух граничных точках х 0 =а, х n =b - коэффициент 1.

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

Если функция f(x) имеет на непрерывную производную четвертого порядка, то абсолютная величина погрешности формулы Симпсона не больше чем

где М - наибольшее значение на отрезке . Так как n 4 растет быстрее, чем n 2 , то погрешность формулы Симпсона с ростом n уменьшается значительно быстрее, чем погрешность формулы трапеций.

Вычислим интеграл

Этот интеграл легко вычисляется:

Возьмем n равным 10, h=0.1, рассчитаем значения подынтегральной функции в точках разбиения, а также полуцелых точках.

По формуле средних прямоугольников получим I прям =0.785606 (погрешность равна 0.027%), по формуле трапеций I трап =0.784981 (погрешность около 0,054. При использовании метода правых и левых прямоугольников погрешность составляет более 3%.

Для сравнения точности приближенных формул вычислим еще раз интеграл

но теперь по формуле Симпсона при n=4. Разобьем отрезок на четыре равные части точками х 0 =0, х 1 =1/4, х 2 =1/2, х 3 =3/4, х 4 =1 и вычислим приближенно значения функции f(x)=1/(1+x) в этих точках: у 0 =1,0000, у 1 =0,8000, у 2 =0,6667, у 3 =0,5714, у 4 =0,5000.

По формуле Симпсона получаем

Оценим погрешность полученного результата. Для подынтегральной функции f(x)=1/(1+x) имеем: f (4) (x)=24/(1+x) 5 , откуда следует, что на отрезке . Следовательно, можно взять М=24, и погрешность результата не превосходит величины 24/(2880 4 4)=0.0004. Сравнивая приближенное значение с точным, заключаем, что абсолютная ошибка результата, полученного по формуле Симпсона, меньше 0,00011. Это находится в соответствии с данной выше оценкой погрешности и, кроме того, свидетельствует, что формула Симпсона значительно точнее формулы трапеций. Поэтому формулу Симпсона для приближенного вычисления определенных интегралов используют чаще, чем формулу трапеций.

Использование трех точек для интерполирования подынтегрального выражения позволяет использовать параболическую функцию (полином второй степени). Это приводит к формуле Симпсона приближенного вычисления интеграла.

Рассмотрим произвольный интеграл

Воспользуемся заменой переменной таким образом, чтобы границы отрезка интегрирования вместо стали [-1,1], для этого введем переменную z:

Тогда и

Рассмотрим задачу интерполирования полиномом второй степени (параболой) подынтегральной функции, используя в качестве узлов три равноудаленные узловые точки – z = -1, z = 0, z = +1 (шаг равен 1, длина отрезка интегрирования равна 2). Обозначим соответствующие значения подынтегральной функции в узлах интерполяции

Система уравнений для нахождения коэффициентов полинома

Проходящего через три точки , и

примет вид

или

Коэффициенты легко могут быть получены

Вычислим теперь значение интеграла от интерполяционного многочлена

Путем обратной замены переменной вернемся к исходному интегралу. Учтем, что

Получим формулу Симпсона для произвольного интервала интегрирования:

При необходимости, исходный отрезок интегрирования может быть разбит на N сдвоенных отрезков, к каждому из которых применяется формула Симпсона. Шаг интерполирования при этом составит

Для первого отрезка интегрирования узлами интерполирования будут являться точки a, a+h, a+2h, для второго – a+2h, a+3h, a+4h, третьего a+4h, a+5h, a+6h и т.д. Приближенное значение интеграла получается суммированием N площадей:

В данную сумму входят одинаковые слагаемые (для внутренних узлов с четным значением индекса - 2i). Поэтому можно перегруппировать слагаемые в этой сумме таким образом

Что эквивалентно

Так как

Погрешность этого приближенного метода уменьшается пропорционально длине шага интегрирования в четвертой степени, т.е. при увеличении числа интервалов вдвое ошибка уменьшается в 16 раз

Увеличение точности

Здесь мы рассмотрим так называемый процесс Эйткена. Он дает возможность оценить погрешность метода и указывает алгоритм уточнения результатов. Расчет проводится последовательно три раза при различных шагах разбиения h 1 , h 2 , h 3 , причем их отношения постоянны: h 2 / h 1 = h 3 / h 2 = q (например, при делении шага пополам q=0.5). Пусть в результате численного интегрирования получены значения интеграла I 1 , I 2 , I 3 . Тогда уточненное значение интеграла вычисляется по формуле

а порядок точности используемого метода численного интегрирования определяется соотношением

.

Уточнение значения интеграла можно также проводить методом Рунге-Ромберга.

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

На практике нередко встречаются случаи, когда подынтегральная функция меняется по-разному на отдельных участках отрезка интегрирования. Это обстоятельство требует такой организации экономичных численных алгоритмов, при которой они автоматически приспосабливались бы к характеру изменения функции. Такие алгоритмы называются адаптивными (приспосабливающимися). Они позволяют вводить разные значения шага интегрирования на отдельных участках отрезка интегрирования. Это дает возможность уменьшить машинное время без потери точности результатов расчета. Подчеркнем, что этот подход используется обычно при задании подынтегральной функции y=f(x) в виде формулы, а не в табличном виде.

Рассмотрим принцип работы адаптивного алгоритма. Первоначально отрезок разбиваем на n частей. В дальнейшем каждый такой элементарный отрезок делим последовательно пополам. Окончательное число шагов, их расположение и размеры зависят от подынтегральной функции и допустимой погрешности e .

К каждому элементарному отрезку применяем формулы численного интегрирования при двух различных его разбиениях. Получаем приближения для интеграла по этому отрезку:

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

Процесс деления отрезка пополам и вычисления уточненных значений продолжается до тех пор, пока их разность станет не больше некоторой заданной величины d i, зависящей от e и h:

.

Аналогичная процедура проводится для всех n элементарных отрезков. Величина принимается в качестве искомого значения интеграла. Условия и соответствующий выбор величин d i обеспечивают выполнение условия

x i -1/2 =(x i +x i -1)/2 – середина i -го отрезка

Представим на отрезке [x i -1 , x i ] подынтегральную функцию f (x ) в виде полинома третьей степени P i (x ). Этот полином должен быть равен значениям подынтегральной функции в точках сетки и в середине отрезка: P i (x i - 1)=f (x i -1)– равенство полинома значению функции на левой границе i -го отрезка,

P i (x i- 1/2) = f (x i -1/2), P i (x i ) = f (x i ).

Такой полином можно записать, например, следующим образом:

P i (x )=a +b(x-x i -1)+c(x-x i -1)(x-x i -1/2),

здесь a , b, c – неизвестные коэффициенты, подлежащие определению.

Введем обозначение для ширины i -го отрезка: h i =x i -x i -1 ,

тогда (x-x i -1/2)= h i /2, а (x i -1/2 -x i -1)= h i /2.

Запишем значения полинома на левой, правой границах и в середине i -го отрезка

P i (x i ) = a +b*h i + c*h i *h i /2 = f (x i )= f i (1)

P i (x i- 1) = a = f (x i -1)= f i -1 (2)

P i (x i- 1/2)= f (x i -1/2)= a +b*h i /2 = f i -1/2 (3)

Из соотношения (2) следует a = f i -1 ,

из выражения (3) легко увидеть, что b= h i (f i -1/2 - f i )/2,

из выражения (1) получаем c=2 (f i -a -b h i )/h i 2 , подставим в выражение для коэффициента c выражения для коэффициентов a и b, в результате получим:

c=2(f i - f i -1) /h i 2 (2/h i )(2/h i )(f i -1/2 - f i -1 ) ,

c=2 [f i - f i -1 -2 f i -1/2 +2 f i -1 ] /h i 2 ,

c=2 [f i - 2 f i -1/2 +f i -1 ] /h i 2 .

Подставим найденные коэффициенты a , b, c в выражение для полинома:

P i (x )= f i -1 + 2(f i -1/2 - f i -1)( x -x i -1) /h i + 2 [f i - 2 f i -1/2 +f i -1 ] ( x -x i -1) ( x -x i -1/2)/h i 2

Перейдем от переменной x к переменной t= x -x i -1

Тогда dt = dx , а при x = x i -1 ; t=0, при x = x i ; t=h i при

x = x i -1/2 =x -( x i -x i -1)/2= x -x i /2 -x i -1 /2= x - x i -1 +x i -1 /2-x i /2=t-h i /2

Тогда на i -ом интервале значение интеграла с учетом введенных обозначений, можно записать:

Подставим в выражение для значения коэффициентов a,b и c

Таким образом,

S i – представляет собой значение интеграла на i -ом отрезке. Для получения интеграла на отрезке от a до b, необходимо сложить все S i

Если h i =h для любого i =1,…, N, тогда и формулу Симпсона можно упростить

(4)

Формулу (4) можно упростить, для этого раскроем скобки в выражении под знаком суммирования

Выделим из первой суммы значение функции в точке x =a

,

а из последней суммы – значение функции в точке x =b

В результате получаем рабочую формулу Симпсона для равномерной сетки.

Учтем, что , , получим окончательное выражение для формулы Симпсона

В первой сумме формулы (5) вычисляют сумму значений функции во всех внутренних узлах отрезка , вторая сумма вычисляет сумму значений функции в средних точках i -ых отрезков.



Если середины отрезков включить в сетку наряду с узлами, тогда новый шаг h 0 = h/2 = (b-a)/(2*n), а формула (5) может быть записана в виде:

Рассмотрим . Значение данного интеграла легко найти аналитически и оно равно -0,75. Метод Симпсона для подынтегральной функции в виде полинома степени 3 и ниже дает точное значение.

Алгоритм вычисления этого интеграла методом Симпсона (формула (5)).

цикл по i от 1 до n-1

конец цикла

цикл по I от 1 до n

конец цикла

s=h*(f0+2*s1+4*s2+fn)/6

функция f1

параметры x

возврат x^3+3*x^2 + x*4 - 4

Пример программы вычисления интеграла методом Симпсона на языке VFP (по формуле (6)):

SET DECIMALS TO 10

? "I=",simpson(0,2,20)

PROCEDURE simpson

PARAMETERS a,b,n

S_четные=0

S_нечетные=0

for x=a+h TO b-h STEP 2*h

S_нечетны = S_нечетные + 4*f(x)

for x=a+2*h TO b-h STEP 2*h

S_четные = S_четные + 2*f(x)

S=f(a)*h/3+(S_четные+S_нечетные)*h/3+f(b)*h/3

Пример решения на языке VBA :

"процедура проверки правильности вычисления значения интеграла по его первообразной

s_четные = 0

s_нечетные = 0

For x = a + h To b - h Step 2 * h

s_нечетные = s_нечетные + 4 * f(x)

Debug.Print "s_нечетные = " & s_нечетные

For x = a + 2 * h To b - h Step 2 * h

s_четные = s_четные + 2 * f(x)

Debug.Print "s_четные=" & s_четные

s = h / 3 * (f(a) + (s_четные + s_нечетные) + f(b))

Debug.Print "Метод Симпсона: s= " & s

Debug.Print "Значение первообразной: s_test= ” & s_test(b-a)

Результат работы программы на VBA:

s_нечетные = 79,9111111111111

s_четные=36,0888888888889

Метод Симпсона: s= 2,66666666666667

Значение первообразной: s_test= 2,66666666666667

Контрольные вопросы



1. Что такое определенный интеграл?

2. Привести алгоритм метода прямоугольников.

3. На интервале функция f(x) монотонно возрастает. I 1 – значение интеграла функции f(x) на отрезке , вычисленное по методу левых прямоугольников, I 0 – значение интеграла функции f(x) на отрезке , вычисленное по методу средних прямоугольников. Будут ли отличаться значения интеграла, вычисленные этими методами? Если значения различны, то какое из них больше? Чем определяется разница?

4. Оценить погрешность для вычисления интеграла методом правых прямоугольников для монотонно убывающей функции.

5. Привести алгоритм метода трапеций.

6. Привести алгоритм метода Симпсона.

7. Как определить погрешность вычисления интеграла итерационными методами?

8. Какой из методов имеет наименьшую погрешность вычисления определенного интеграла?

9. Получить формулу метода Симпсона.

Задания

Вычислить следующие интегралы методами: прямоугольников, трапеций, Симпсона с точностью 0,001 и оценить погрешность результатов вычислений этими методами.

2.

3.

4.

5.

10.

11.

14.

18.