Русские Блоги
Вычисление интеграла с помощью модуля Scipy в Python
Метод получения интеграла в модуле Scipy на Python:
SciPy реализует базовое использование функции, чтобы найти интеграл функции, интеграл, продвинутая математика имеет много разговоров, основной смысл — найти площадь под кривой.
Среди них rn можно рассматривать как отклонение, которым обычно можно пренебречь, а wi можно рассматривать как вес.
SciPy предоставляет множество функций для вычисления различных типов интегралов, которые можно разделить на две категории в зависимости от разницы входных параметров: одна — это вход известной функции и верхний и нижний пределы интеграла; другая — набор входных точек.Это применимо к некоторым данным, собранным после физической реализации, но функция не может быть определена, но есть много точек данных, тогда какая площадь под огибающей этих точек также является интегральной проблемой, поэтому существует точка, установленная в функции SciPy Integral, параметр функции — это массив или список в форме.
1. Интеграл с известным типом функции
В этом разделе показано, как вычислить интеграл в SciPy в виде нескольких вопросов.
- Вопрос 1: Здесь предположим, что функция f (x) = x + 1, верхний и нижний пределы интеграла равны [1,2] математическое выражение:
Вы можете использовать подмодуль, интегрированный в модуль Scipyчетырехъядерная функцияНайти расчетное значение этой математической задачи.
Результат выполнения программы:
Результат выполнения программы:
Когда программа запущена:
Результат — бесконечность (бесконечность, бесконечность) и есть ошибка деления на 0! немного отредактируйте:
Мы можем нарисовать визуальную кривую этой функции:
import numpy as np, matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig, ax = plt.subplots(figsize=( 8, 3))
x = np.linspace(- 1, 1, 10000)
ax.plot(x, f(x), lw= 2)
ax.fill_between(x, f(x), color= ‘green’, alpha= 0.5)
ax.set_xlabel( " x x x ", fontsize= 18)
ax.set_ylabel( " f ( x ) f(x) f ( x ) ", fontsize= 18)
ax.set_ylim( 0, 25)
plt.show()
Получите следующую таблицу результатов:
2 Дайте интеграл от множества точек
В случае, если интегральная функция не может быть подтверждена, некоторые последовательности также могут быть интегрированы.
- Вопрос 4： Интегрировать Но есть 10 примеров данных этой функции, функция quad не f (x) = x, но соответствующие (xi,yi)。
3 Несколько точек
Двойной интеграл в SciPy можно вычислить с помощью функции dblquad, а тройной интеграл можно вычислить с помощью функции tplquad. Функцию nquad можно использовать для множественных интегралов от f (x1, x2, . xn).
- Двойная интегральная функция dblquad для вычисления, предположим, что есть функция f (x, y) необходимо вычислить свой двойной интеграл.
Exploring Integrals in Python
Want to learn more about integration and math behind it? Let’s explore.
Entering the Data science world will eventually lead you to explore and learn about the math behind it. Sometimes it can be overwhelming and challenging to understand math, so on that note, in this article, we will cover the basics of integration while playing in Python.
Without further ado, let’s start.
For better understanding, we will use a simple quadratic function with null points <-3,3>— points where the curve intersects the x-axis.
Graph of the function :
As you can see, it is nothing complicated, so let’s investigate more. We can ask ourselves, what about the white space under the curve? Here comes integration in our equation. It can be used for finding volumes, areas, central points, and many other useful values, but it is mainly known as the area under the curve (function).
A little visualization in case you didn’t understand:
Basics of integration
With understanding what integral represents, we can do a little bit of calculus. When resolving integral, we are doing the following — finding the function, which derivation is function before integration. I know it sounds a little bit confusing, trust me it is not difficult.
For a start, let’s take a simple constant function:
Now, integration calculus is simple :
I think you get the idea behind it. We added 1 to the exponent of variable x and then divided it with the same exponent value.
Note: Please, don’t get confused by the solution. I couldn’t just write ax because little mathematician in me will start to scream out. Constant C is called constant of integration, and it represents all numbers that can be in the solution. We don’t know if our solution is (ax + 2) or (ax + 3), and so on — therefore, we write + C to cover all possible variations.
Calculating integral of x :
Again we implemented the same method as before — adding 1 to the exponent and dividing entire value with exponent.
Indefinite and definite integrals
Indefinite integral we already calculated above, it is a function that answers the question, “What function, after derivation, gives f(x)?”. In contrast, a definite integral is a number that represents the area under the curve from x=a to x=b.
The calculation is the same until the last step; here is an example:
- For indefinite integral — the solution is calculated function :
2. For definite integral — there is additional calculation for boundaries that are set for integral [a,b]=[-2,2]. You can see that boundaries are input in the solution calculated above.
First, we put the upper boundary b in the x and write it down. Next, we put the lower boundary a in the x. Finally, we subtract the first and last expression to get the final result.
Are you good so far? Hold on. Here comes the exciting part.
Implementation in Python
With a little bit of research, you can find many Python libraries for integral calculation. In this “journey” we are interested in how functions look after integration but also results if implementing boundaries for the same functions. Yes, you can do it from scratch, but it would be a lot of work. For me, the most elegant one is SymPy. Installation of it is straightforward, just write in your notebook following code, and installation is done.
Now, let’s have fun. Remember the first example from the beginning of the article? Let me remind you :
Integral Calculus in Python
This post shows how to perform integral calculus of continuous and limited real functions of real variables in Python through the use of common Python libraries frequently used in scientific applications. The integral calculation techniques here are both primarily numerical since this site deals with computation, however some analytical techniques are also shown.
The post is organized by examples: each paragraph contains an example of an integral to compute and the related Python code snippet that calculates it using an appropriate library.
All of the various code snippets described in this post require Python version 3 and the NumPy library, while individually they require an additional library (and its dependencies, if any) between SciPy and SymPy.
We thank Prof. Fausta D'Acunzo of Preparazione 2.0 for theoretical support provided on multi-variable integral calculus.
To get the code see the Full code download paragraph at the bottom of this post.
Integration via SciPy
Integral of function of one variable (with finite extremes)
In integral calculus, the definite integral is an operator that, given a real-valued function of a real-valued variable and an interval $[a,b]$ (subset of the domain), associates to the function the area subtended by its graph in the interval $[a,b]$.
The SciPy library provides several numerical methods for computing the integral of such functions; the purpose of this chapter is to present a series of demos to show "by examples" the use of such methods; for complete documentation of these methods, the reader is invited to consult the official documentation of SciPy.
Let the following integral of a function of one variable be given: $$ \int_<1>^ <5>2 x e^ <-x>\,dx $$ whose analytical solution is $ \approx 1.3907 $ verifiable online via Wolfram Alpha.
Calculating the integral with quad
Below is the example of Python code that calculates the integral using the quad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with fixed_quad
Below is the example of Python code that calculates the integral using the fixed_quad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with quadrature
Below is the example of Python code that calculates the integral using the quadrature function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with romberg
Below is the example of Python code that calculates the integral using the romberg function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with trapezoid
The trapezoid function is a fixed-sample function integration method, so the code first discretizes the integration interval evenly spaced and for all discrete values of $x$ it computes the corresponding values of $y$ and then passes the two sets of discrete values $x$ and $y$ to the integration method.
Below is the example of Python code that calculates the integral using the trapezoid function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with cumulative_trapezoid
The function cumulative_trapezoid is also a fixed-sample function integration method, and so what was said about trapezoid applies.
Below is the example of Python code that calculates the integral using of the cumulative_trapezoid function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with simpson
The function simpson is also a fixed-sample function integration method, and so what was said about trapezoid applies.
Below is the example of Python code that calculates the integral using the simpson function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Integral of function of one variable with quad (with extreme infinity)
Set the following integral of a function of one variable be given: $$ \int_<1>^ <+\infty>2 x e^ <-x>\,dx $$ whose analytical solution is $ \approx 1.4715 $ verifiable online via Wolfram Alpha.
Below is the example of Python code that calculates the integral using the quad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculation of the length of a planar curve arc
Calculation of the length of a planar curve arc expressed in explicit form with quad
Let the following planar curve be given in explicit form: $$ y=e^ <-x^2>$$ by applying the corresponding formula above to calculate the length of the arc between $ x=-1 $ and $ x = 1 $, the integral is obtained: $$ \int_<-1>^ <1>\sqrt<1 + \left(\frac<\,d(e^<-x^2>)><\,dx>\right)^2> \,dx $$ whose analytical solution is $ \approx 2.4089 $ verifiable online via Wolfram Alpha.
Here is the example of Python code that calculates length of a planar curve arc expressed in explicit form using the quad function of SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculation of the length of a planar curve arc expressed in parametric form with quad
Let the following planar curve be given in parametric form: $$ x(t)=cos^3 t $$ $$ y(t)=sin^3 t $$ by applying the corresponding formula above to calculate the length of the arc between $ t=0 $ e $ t = 2\pi $, the integral is obtained: $$ \int_<0>^ <2 \pi>\sqrt<\left(\frac<\,d (cos^3 t)><\,dt>\right)^2 + \left(\frac<\,d (sin^3 t)><\,dt>\right)^2> \,dt $$ whose analytical solution is $ 6 $ verifiable online via Wolfram Alpha.
Here is the example of Python code that calculates length of a planar curve arc expressed in parametric form using the quad function of SciPy library: whose output is:
Here is the link to the code on GitHub.
Double integral of a function of two variables
In integral calculus, the definite double integral is an operator that, given a real-valued function of two real-valued variables and a set included in the domain, associates to the function the volume of the solid (called cylindroid) between the surface described by the function and the plane containing the given set.
Let the following double integral of a function of two variables be given: $$ \int_<1>^ <5>\int_
Calculating the integral with dblquad
Below is the example of Python code that calculates the integral using the dblquad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with nquad
Below is the example of Python code that calculates the integral using the nquad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Double integral of a function of two variables with nquad (other example)
Let the following double integral of a function of two variables be given: $$ \int_<1>^ <+\infty>\int_<1>^ <+\infty>2 x y e^ <-x y>\,dx dy $$ whose analytical solution is $ \approx 1.17453 $ verifiable online via Wolfram Alpha.
Below is the example of Python code that calculates the integral using the nquad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Triple Integral of a function of three variables
Let the following triple integral of a function of three variables be given: $$ \int_<1>^ <2>\int_
Calculating the integral with tplquad
Below is the example of Python code that calculates the integral using the tplquad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Calculating the integral with nquad
Below is the example of Python code that calculates the integral using the nquad function of the SciPy library: whose output is:
Here is the link to the code on GitHub.
Integration via SymPy
Integral of function of one variable (with finite extremes)
In integral calculus, the definite integral is an operator that, given a real-valued function of a real-valued variable and an interval $[a,b]$ (subset of the domain), associates to the function the area subtended by its graph in the interval $[a,b]$.
The SciPy library provides several numerical methods for computing the integral of such functions; the purpose of this chapter is to present a series of demos to show "by examples" the use of such methods; for complete documentation of these methods, the reader is invited to consult the official documentation of SciPy.
Let the following integral of a function of one variable be given: $$ \int_<1>^ <5>2 x e^ <-x>\,dx $$ whose analytical solution is $ \approx 1.3907 $ verifiable online via Wolfram Alpha.
Integral of a function of one variable with integrate(f, (x, a, b))
Below is the example of Python code that calculates the integral using the integrate(f, (x, a, b)) function of the SymPy library: whose output is:
Here is the link to the code on GitHub.
Integral of a function of one variable with integrate(f, (x, a, b)) via indefinite integral
Below is the example of Python code that calculates the integral using the integrate(f, x) function of the SymPy library: whose output is:
The program first calculates the indefinite integral and then applying the fundamental theorem of integral calculus, calculates the value of the integral.
Here is the link to the code on GitHub.
Double integral of a function of two variables
In integral calculus, the definite double integral is an operator that, given a real-valued function of two real-valued variables and a set included in the domain, associates to the function the volume of the solid (called cylindroid) between the surface described by the function and the plane containing the given set.
Let the following double integral of a function of two variables be given: $$ \int_<1>^ <4>\int_
Double Integral of a function of two variables with integrate(f, (x, xa, xb), (y, ya, yb))
Below is the example of Python code that calculates the integral using the integrate(f, (x, xa, xb), (y, ya, yb)) of the SymPy library: whose output is:
Here is the link to the code on GitHub.
Download of the complete code
The complete code is available at GitHub.
These materials are distributed under MIT license; feel free to use, share, fork and adapt these materials as you see fit.
Also please feel free to submit pull-requests and bug-reports to this GitHub repository or contact me on my social media channels available on the top right corner of this page.
Вычисление определённых интегралов: базовые алгоритмы
В этой публикации описаны простейшие методы вычисления интегралов функций от одной переменной на отрезке, также называемые квадратурными формулами. Обычно эти методы реализованы в стандартных математических библиотеках, таких как GNU Scientific Library для C, SciPy для Python и других. Публикация имеет целью продемонстрировать, как эти методы работают «под капотом», и обратить внимание на некоторые вопросы точности и производительности алгоритмов. Также хотелось бы отметить связь квадратурных формул и методов численного интегрирования обыкновенных дифференциальных уравнений, о которых хочу написать ещё одну публикацию.
Определение интеграла
Интегралом (по Риману) от функции на отрезке называется следующий предел:
где — мелкость разбиения, , , — произвольное число на отрезке .
Если интеграл от функции существует, то значение предела одно и то же вне зависимости от разбиения, лишь бы оно было достаточно мелким. 
Более наглядно геометрическое определение — интеграл равен площади криволинейной трапеции, ограниченной осью 0x, графиком функции и прямыми x = a и x = b (закрашенная область на рисунке).
Квадратурные формулы
Определение интеграла (1) можно переписать в виде
где — весовые коэффициенты, сумма которых должна быть равна 1, а сами коэффициенты — стремиться к нулю при увеличении числа точек, в которых вычисляется функция.
Выражение (2) — основа всех квадратурных формул (т.е. формул для приближенного вычисления интеграла). Задача состоит в том, чтобы выбрать точки и веса таким образом, чтобы сумма в правой части приближала требуемый интеграл как можно точнее.
Вычислительная задача
Задана функция , для которой есть алгоритм вычисления значений в любой точке отрезка (имеются в виду точки, представимые числом с плавающей точкой — никаких там функций Дирихле!).
Требуется найти приближённое значение интеграла .
Решения будут реализованы на языке Python 3.6.
Для проверки методов используется интеграл .
Кусочно-постоянная аппроксимация
Идейно простейшие квадратурные формулы возникают из применения выражения (1) «в лоб»:
Т.к. от метода разбиения отрезка точками и выбора точек значение предела не зависит, то выберем их так, чтобы они удобно вычислялись — например, разбиение возьмём равномерным, а для точек вычисления функции рассмотрим варианты: 1) ; 2) ; 3) .
Получаем методы левых прямоугольников, правых прямоугольников и прямоугольников со средней точкой, соответственно.
Для анализа производительности квадратурных формул построим график погрешности в координатах «число точек — отличие численного результата от точного».
Что можно заметить:
- Формула со средней точкой гораздо точнее, чем с правой или левой точками
- Погрешность формулы со средней точкой падает быстрее, чем у двух остальных
- При очень мелком разбиении погрешность формулы со средней точкой начинает возрастать
Первые два пункта связаны с тем, что формула прямоугольников со средней точкой имеет второй порядок аппроксимации, т.е. , а формулы правых и левых прямоугольников — первый порядок, т.е. .
Возрастание погрешности при измельчении шага интегрирования связано с нарастанием погрешности округления при суммировании большого числа слагаемых. Эта ошибка растёт как , что не даёт при интегрировании достигнуть машинной точности.
Вывод: методы прямоугольников с правой и левой точками имеют низкую точность, которая к тому же медленно растёт с измельчением разбиения. Поэтому они имеют смысл разве что в демонстрационных целях. Метод прямоугольников со средней точкой имеет более высокий порядок аппроксимации, что даёт ему шансы на использование в реальных приложениях (об этом чуть ниже).
Кусочно-линейная аппроксимация
Следующий логический шаг — аппроксимировать интегрируемую функцию на каждом из подотрезков линейной функцией, что даёт квадратурную формулу трапеций:
Иллюстрация метода трапеций для n=1 и n=2.
В случае равномерной сетки длины всех отрезков разбиения равны, и формула имеет вид
Построив график ошибки от числа точек разбиения, убеждаемся, что метод трапеций тоже имеет второй порядок аппроксимации и вообще даёт результаты, слабо отличающиеся от метода прямоугольников со средней точкой (в дальнейшем — просто метод прямоугольников). 
Контроль точности вычисления
Задание в качестве входного параметра числа точек разбиения не слишком практично, поскольку обычно требуется вычислить интеграл не с заданной плотностью разбиения, а с заданной погрешностью. Если подынтегральная функция известна наперёд, то можно оценить погрешность заранее и выбрать такой шаг интегрирования, чтобы заданная точность заведомо достигалась. Но так редко бывает на практике (и вообще, не проще ли при известной наперёд функции и сам интеграл протабулировать наперёд?), поэтому необходима процедура автоматической подстройки шага под заданную погрешность.
Как это реализовать? Один из простых методов оценки погрешности — правило Рунге — разность значений интегралов, рассчитанных по n и 2n точкам, даёт оценку погрешности: . Метод трапеций удобнее для удвоения мелкости разбиения, чем метод прямоугольников с центральной точкой. При расчёте методом трапеций для удвоения числа точек нужны новые значения функции только в серединах отрезков предыдущего разбиения, т.е. предыдущее приближение интеграла можно использовать для вычисления следующего.
Метод прямоугольников не требует вычислять значения функции на концах отрезка. Это означает, что его можно использовать для функций, имеющих на краях отрезка интегрируемые особенности (например, sinx/x или x -1/2 от 0 до 1). Поэтому показанный далее метод экстраполяции будет работать точно так же и для метода прямоугольников. Отличие от метода трапеций лишь в том, что при уменьшении шага вдвое отбрасывается результат предыдущих вычислений, однако можно утроить число точек, и тогда предыдущее значение интеграла также можно использовать для вычисления нового. Формулы для экстраполяции в этом случае необходимо скорректировать на другое соотношение шагов интегрирования.
Отсюда получаем следующий код для метода трапеций с контролем точности:
С таким подходом подынтегральная функция не будет вычисляться по нескольку раз в одной точке, и все вычисленные значения используются для окончательного результата.
Но нельзя ли при том же количестве вычислений функции добиться более высокой точности? Оказывается, что можно, есть формулы, работающие точнее метода трапеций на той же самой сетке.
Кусочно-параболическая аппроксимация
Следующим шагом аппроксимируем функцию элементами парабол. Для этого требуется, чтобы число отрезков разбиения было чётным, тогда параболы могут быть проведены через тройки точек с абсциссами <(x0=a, x1, x2), (x2, x3, x4), . (xn-2, xn-1, xn=b)>.
Иллюстрация кусочно-параболического приближения на 3 и 5 точках (n=2 и n=3).
Приближая интеграл от функции на каждом из отрезков [xk;xk+2] интегралом от параболической аппроксимации на этом отрезке и считая точки равномерно распределенными (xk+1=xk+h), получаем формулу Симпсона:
Из формулы (4) напрямую получается «наивная» реализация метода Симпсона:
Для оценки погрешности можно использовать точно так же вычисление интеграла с шагами h и h/2 — но вот незадача, при вычислении интеграла с более мелким шагом результат предыдущего вычисления придётся отбросить, хотя половина новых вычислений функции будет в тех же точках, что и раньше.
Бесполезной траты машинного времени, к счастью, можно избежать, если реализовать метод Симпсона более хитроумным образом. Присмотревшись повнимательнее, заметим, что интеграл по формуле Симпсона может быть представлен через два интеграла по формуле трапеций с разными шагами. Яснее всего это видно на базовом случае аппроксимации интеграла по трём точкам :
Таким образом, если реализовать процедуру уменьшения шага вдвое и хранить два последних вычисления методом трапеций, метод Симпсона с контролем точности реализуется более эффективно.
Сравним эффективность метода трапеций и парабол:
Как видим, обоими методами ответ можно получть с достаточно высокой точностью, но количество вызовов подынтегральной функции разительно отличается — метод более высокого порядка эффективнее в 32 раза!
Построив график погрешности интегрирования от числа шагов, можно убедиться, что порядок аппроксимации формулы Симпсона равен четырём, т.е. ошибка численного интегрирования (а интегралы от кубических многочленов с помощью этой формулы вычисляются с точностью до ошибок округления при любом чётном n>0!). 
Отсюда и возникает такой рост эффективности по сравнению с простой формулой трапеций.
Что дальше?
Дальнейшая логика повышения точности квадратурных формул, в целом, понятна — если функцию продолжать приближать многочленами всё более высокой степени, то и интеграл от этих многочленов будет всё точнее приближать интеграл от исходной функции. Этот подход называется построением квадратурных формул Ньютона-Котеса. Известны формулы вплоть до 8 порядка аппроксимации, но выше среди весовых коэффициентов wi в (2) появляются знакопеременные члены, и формулы при вычислениях теряют устойчивость.
Попробуем пойти другим путём. Ошибка квадратурной формулы представляется в виде ряда по степеням шага интегрирования h. Замечательное свойство метода трапеций (и прямоугольников со средней точкой!) в том, что для неё этот ряд состоит только из чётных степеней:
На нахождении последовательных приближений к этому разложению основана экстраполяция Ричардсона: вместо того, чтобы приближать подынтегральную функцию многочленом, по рассчитанным приближениям интеграла строится полиномиальная аппроксимация, которая при h=0 должна давать наилучшее приближение к истинному значению интеграла.
Разложение ошибки интегрирования по чётным степеням шага разбиения резко ускоряет сходимость экстраполяции, т.к. для аппроксимации порядка 2n нужно всего n значений интеграла методом трапеций.
Если считать, что каждое последующее слагаемое меньше предыдущего, то можно последовательно исключать степени h, имея приближения интеграла, рассчитанные с разными шагами. Поскольку приведённая реализация легко позволяет дробить разбиение вдвое, удобно рассматривать формулы для шагов h и h/2.
Легко показать, что исключение старшего члена погрешности формулы трапеций в точности даст формулу Симпсона:
Повторяя аналогичную процедуру для формулы Симпсона, получаем:
Если продолжить, вырисовывается такая таблица:
| 2 порядок | 4 порядок | 6 порядок | . |
|---|---|---|---|
| I0,0 | |||
| I1,0 | I1,1 | ||
| I2,0 | I2,1 | I2,2 | |
| . | . | . |
В первом столбце стоят интегралы, вычисленные методом трапеций. При переходе от верхней строки вниз разбиение отрезка становится вдвое мельче, а при переходе от левого столбца вправо повышается порядок аппроксимации интеграла (т.е. во втором столбце находятся интегралы по методу Симпсона и т.д.).
Элементы таблицы, как можно вывести из разложения (5), связаны рекуррентным соотношением:
Погрешность приближения интеграла можно оценить по разности формул разных порядков в одной строке, т.е.
Применение экстраполяции Ричардсона вместе с интегрированием методом трапеций называется методом Ромберга. Если метод Симпсона учитывает два предыдущих значения по методу трапеций, то метод Ромберга использует все ранее вычисленные методом трапеций значения для получения более точной оценки интеграла.
Дополнительный метод добавляется в класс Quadrature
Проверим, как работает аппроксимация высокого порядка:
Убеждаемся, что, по сравнению с методом парабол, число вызовов подынтегральной функции снизилось ещё в 8 раз. При дальнейшем увеличении требуемой точности преимущества метода Ромберга проявляются ещё заметнее:
alt=»image» />
Некоторые замечания
Замечание 1. Количество вызовов функции в этих задачах характеризует число суммирований при вычислении интеграла. Уменьшение числа вычислений подынтегрального выражения не только экономит вычислительные ресурсы (хотя при более оптимизированной реализации и это тоже), но и уменьшает влияние погрешностей округления на результат. Так, при попытке вычислить интеграл тестовой функции метод трапеций зависает при попытке достигнуть относительной точности 5×10 -15 , метод парабол — при желаемой точности 2×10 -16 (что является пределом для чисел в двойной точности), а метод Ромберга справляется с вычислением тестового интеграла вплоть до машинной точности (с ошибкой в младшем бите). То есть, повышается не только точность интегрирования при заданном числе вызовов функции, но и предельно достижимая точность вычисления интеграла.
Замечание 2. Если метод сходится при задании некоторой точности, это не означает, что вычисленное значение интеграла имеет ту же самую точность. В первую очередь, это относится к случаям, когда задаваемая погрешность близка к машинной точности.
Замечание 3. Хотя метод Ромберга для ряда функций работает почти магическим образом, он предполагает наличие у подынтегральной функции ограниченных производных высоких порядков. Это значит, что для функций с изломами или разрывами он может оказаться хуже простых методов. Например, проинтегрируем f(x)=|x|:
Замечание 4. Может показаться, что чем выше порядок аппроксимации, тем лучше. На самом деле, лучше ограничить число столбцов таблицы Ромберга на уровне 4-6. Чтобы понять это, посмотрим на формулу (6). Второе слагаемое представляет собой разность двух последовательных элементов j-1-го столбца, поделенную на примерно 4 j . Т.к. в j-1-м столбце находятся аппроксимации интеграла порядка 2j, то сама разность имеет порядок (1/ni) 2j
4 —ij . C учётом деления получается
7 второе слагаемое в (6) теряет точность после приведения порядков при сложении чисел с плавающей точкой, и повышение порядка аппроксимации может вести к накоплению ошибки округления.
Замечание 5. Желающие могут ради интереса применить описанные методы для нахождения интеграла и эквивалентного ему . Как говорится, почувствуйте разницу.
Заключение
Представлено описание и реализация базовых методов численного интегрирования функций на равномерной сетке. Продемонстрировано, как с помощью несложной модификации получить на базе метода трапеций класс квадратурных формул по методу Ромберга, что значительно ускоряет сходимость численного интегрирования. Метод хорошо работает для интегрирования «обычных» функций, т.е. слабо меняющихся на отрезке интегрирования, не имеющих особенностей на краях отрезка (см. Замечание 5), быстрых осцилляций и т.д.
Продвинутые методы численного интегрирования для более сложных случаев можно найти в книгах из списка литературы (в [3] — с примерами реализации на C++).