В данной лекции будет рассматриваться задача численного интегрирования. Формулы
численного интегрирования функций одного переменного называют квадратурными формулами. Задача приближенного вычисления определенного интеграла (на отрезке или по многомерной области) фактически разбивается на две самостоятельные подзадачи. Первая — это интегрирование
Вторая задача — подсчет значения определенного интеграла от известной
функции. При этом самая ресурсоемкая операция с точки зрения вычислений — подсчет значения функции. Желательно построить численный метод, позволяющий получать как можно более высокую точность при наименьшем количестве
вычислений, при этом выбор узлов квадратурных формул целиком в руках вычислителя. В этом случае наиболее эффективными окажутся
Простейшую квадратурную формулу (формулу численного интегрирования) можно получить следующим образом. Пусть необходимо вычислить интеграл
$$I = \int\limits_{a}^{b}{f(t)}dt.$$Положим, что f(t) на рассматриваемом отрезке [a, b] не изменяется ( $$f(t) \approx const$$ ). Тогда $$I = f(\xi )(b - a), \xi \in \left[{a, b}\right].$$ Если $$$ \xi = \frac{a + b}{2} $,$$ то получим формулу прямоугольников с центральной точкой
Конечно, для константы приведенная выше формула точна — говорят, что
построенная квадратурная формула будет точна на полиномах степени 0. Легко можно доказать, что формула прямоугольников с центральной точкой будет давать точное значение и в случае линейной функции. Для всех других функций эту формулу будем рассматривать как приближенную.
Если предположить, что функция f(t) на отрезке интегрирования [a, b] достаточно близка к линейной, то можно заменить приближенное значение интеграла i площадью трапеции с высотой (b - a) и основаниями f(a) и f(b). Тогда получается формула трапеций
Введем на отрезке интегрирования сетку, определим значения функции в узлах сетки. Узлы в дальнейшем будем именовать узлами квадратурной формулы (или квадратуры). Пусть, как и в задаче интерполяции, имеется совокупность узлов $$\left\{{t_n}\right\}_{n = 0}^{N}, t_n = a + n{\tau}, {\tau}= (b - a)/{N}, t \in \left[{a, b}\right].$$ Пусть также задана таблица $$f_n = \left\{{f(t_n)}\right\}_{n = 0}^{N}.$$ Отрезок [tk, tk + 1] далее иногда будем называть элементарным отрезком.
Заменим подынтегральную функцию ее интерполяционным полиномом в форме Лагранжа. Будем полагать, что
$$\int\limits_{a}^{b}{f(t)} dt \approx \int\limits_{a}^{b}{L_n} (t)dt.$$Рассмотрим некоторые частные случаи.
Формула трапеций. На отрезке [tk, tk + 1] проводим замену подынтегральной функции интерполяционным полиномом первой степени:
после чего, выполнив интегрирование по элементарному отрезку, получим приближенное значение интеграла на [tk, tk + 1]:
После суммирования интегралов по всем элементарным отрезкам [tk, tk + 1] получаем формулу трапеций для отрезка [a, b]:
На равномерной сетке (сетке с равноотстоящими узлами) при $$\tau _{k} = \tau = (b - a)/N$$ полученная формула принимает вид
$$$ I \approx \frac{{\tau}}{2}\sum\limits_{k = 0}^{N - 1}{\left[{f(t_k) + f(t_{k + 1})}\right]} = \frac{{\tau}}{2}\left[{f(t_0 ) + 2f(t_1 ) + \ldots + 2f(t_{N - 1}) + f(t_N)}\right]. $$$Формула Симпсона. Заменим подынтегральную функцию f(t) на отрезке [tk - 1, tk] интерполяционным полиномом (в форме Лагранжа) второй степени. Для простоты положим $$\tau = (b - a)/N = t_{k} - t_{k - 1} = const$$ для всех k — сетка на отрезке интегрирования равномерная. Тогда
После вычисления интеграла от полинома получим приближенное значение интеграла по элементарному отрезку
$$$ I_k \approx \frac{{\tau}}{6}\left[{f_{k - 1} + 4f_{k - 1/2} + f_k}\right]. $$$Суммируя по всем элементарным отрезкам [tk - 1, tk], получим
где $$$ f_k = f(t_k), f_{k + 1/2} = f(\frac{{t_k + t_{k + 1}}}{2}) $.$$ Формулу Симпсона можно также записать, не используя дробных индексов:
$$$ I \approx \frac{{\tau}}{6}(f_0 + 4f_1 + 2f_2 + 4f_3 + \ldots + 2f_{N - 2} + 4f_{N - 1} + f_N), $$$если локальную формулу получать путем интегрирования интерполяционного
полинома второй степени по отрезку [tk - 1, tk + 1]:
где F(t) — [tk - 1, tk + 1] по точкам tk - 1, tk, tk + 1. В этом случае N — число разбиений отрезка на элементарные отрезки — должно быть четным.
Еще одна используемая на практике
Конечно, существуют и формулы интерполяционного типа более высоких порядков. Они не применяются на практике по следующим обстоятельствам. Любая формула интерполяционного типа записывается в виде
$$I \approx \sum\limits_{k = 0}^{M}{\alpha_k f_k}.$$Во всех приведенных выше формулах коэффициенты $$\alpha_k$$ были
положительными. Такие квадратурные формулы называются
Для степени интерполяционного полинома более 7 среди коэффициентов встречаются отрицательные. Д.Пойа показал, что
где $$\alpha_{nk}$$ — веса квадратурной формулы, получающейся
при замене подынтегрального выражения интерполяционным полиномом степени n. Такое увеличение суммы абсолютных значений коэффициентов связано с быстрым ростом
Подробнее о свойствах
Погрешность квадратурных формул может быть оценена, например, с использованием остаточного члена интерполяционного полинома:
$$$ \varepsilon_k = \left|\int\limits_{t_k}^{t_{k + 1}}R_N(t) dt\right| \le \frac{\max\limits_{[t_k, t_{k + 1}]}\left|f^{(N + 1)}(\xi)\right| }{(N + 1)!} {\tau}\max\limits_{[t_k, t_{k + 1}]}\left|{\mathop \Pi\limits_{k = 0}^{N}(t - t_k)}\right|. $$$Из последней формулы следует, что квадратурная формула точна, если
подынтегральная функция является многочленом степени не выше N.
Получим, например, локальную оценку погрешности для формулы трапеций, используя формулу для остаточного члена интерполяционного полинома 1 - го порядка:
$$$ \varepsilon_k \le \int\limits_{t_k}^{t_{k + 1}}{\frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{2}} \max \left|{(t - t_k)(t - t_{k + 1})}\right|dt = \frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{{12}}\tau^3 . $$$Тогда погрешность по всему отрезку [a, b] будет составлять
Другой способ получения погрешности квадратурных формул состоит в следующем. Рассмотрим интеграл по элементарному отрезку
$$\begin{gather*} I_k = \int\limits_{t_k}^{t_{k + 1}}{f(t)}dt, \\ f(t) = f(z) + f^{\prime}(z)(t - z) + \frac{{f^{\prime\prime}}}{2}(t - z)^2 + \ldots + R(t - z), \end{gather*}$$где $$z \in \left[{t_k , t_{k + 1}}\right]$$ — некая опорная точка, тогда для приближенного значения интеграла верно
$$I_k = f(z)\tau_k + \xi\tau^2_k + \eta\tau^3_k + \ldots + \int\limits_{t_k}^{t_{k + 1}} R(t - z)dz, \tau_k = t_{k + 1} - t_k.$$Коэффициенты $$\xi , \eta , ...$$ зависят от производных f'(z), f''(z), ....
С другой стороны, любая из рассмотренных квадратурных формул представима в виде
$$$ \bar {I_k} = \tau_k (af_k + bf_{k + 1/2} + cf_{k + 1}). $$$Заменяя в этой формуле значения функций f в точках fk, fk + 1/2, fk + 1 ее разложением по формуле Тейлора, получим
где $$z \in [t_k, t_{k + 1} ].$$
Сравнивая разложения для $$I_k, \bar {I_k},$$ легко заметить, что вместе с первым слагаемым совпадают и другие слагаемые до (m - 1) - го порядка, так что $$\xi = \xi _{1}, \eta = \eta _{1}, ...$$
Разность же несовпадающих слагаемых будет, очевидно, оценкой погрешности
квадратурной формулы на интервале $$[t_k, t_{k + 1} ]: \varepsilon_k =
\left|I_k - \bar {I_k}\right| \le v \max\limits_{[t_k, t_{k + 1}]}\left|f^ {(m)}\right| \tau^{m + 1}_k,$$ где v — константа.
Если просуммировать локальные погрешности по всем интервалам [ tk, tk + 1 ], то получим оценку погрешности квадратурной формулы по всему отрезку [a, b]:
где $${\tau}= \max\limits_k \tau_k$$ на неравномерной сетке, или $$\tau = (b - a)/N$$ на равномерной. Число m называется порядком точности квадратуры.
Получим теперь погрешность формулы прямоугольников (со средней точкой) для $$\tau _{k} = \tau = const$$:
$$$ \varepsilon \le \frac{b - a}{24} \max\limits_{[a, b ]}\left|{f^{\prime\prime}_{t} (t)}\right|{\tau}^2, $$$погрешность формулы трапеций
$$$ \varepsilon \le \frac{b - a}{12} \max\limits_{[a, b]}\left|{f^{\prime\prime}_{t} (t)}\right|{\tau}^2, $$$погрешность формул Симпсона (с дробными и без дробных индексов соответственно)
$$\begin{multiple} \varepsilon \le \frac{b - a}{2880}\max\limits_{[a, b]} \left|{f_{t}^{(4)} (t)}\right|{\tau}^4, \\ \varepsilon \le \frac{{b - a}}{{180}}\max\limits_{[a, b]} \left|{f_{t}^{(4)} (t)}\right|{\tau}^4 \end{multiple)$$Заметим, что если функция f(t) имеет только три непрерывных
производных, то оценка погрешности формулы Симпсона ухудшается на порядок:
Рассмотрим, для примера двукратный интеграл по прямоугольной области $$\Omega (a\le x\le b, c\le y\le d).$$ Аналогично одномерному случаю, в соответствии с формулой прямоугольников (со средней точкой) заменим функцию
ее значением в точке пересечения диагоналей прямоугольника. В таком случае получим $$$ I = \int\limits_{a}^{b} \int\limits_{c}^{d}f(x, y) dx dy \approx Sf (\frac{a + b}{2}, \frac{c + d}{2}) $,$$ где $$S = \tau _{x}\tau _{y}, \tau _{x} = (b - a), \tau _{y} = (d - c).$$ Если разбить прямоугольник на прямоугольные ячейки Si, то получим аналог формулы прямоугольников со средней точкой в виде суммы по i: $$I \approx \sum\limits_i S_if(x_i, y_i),$$ где Si — площадь i - ой прямоугольной ячейки, xi, yi — координаты точки пересечения ее диагоналей.
Применим теперь формулу Симпсона для вычисления двукратного интеграла путем редукции к методу вычисления одномерного интеграла, зачем представим двойной интеграл как
$$I = \int\limits_{a}^{b}{\int\limits_{c}^{d} {f(x, y) dx dy = \int\limits_{a}^{b}{dx}} } \int\limits_{c}^{d}{f(x, y)} dy.$$Сначала применим формулу Симпсона для вычисления внешнего интеграла:
$$$ I \approx \frac{{\tau_x}}{6}\left[{\int\limits_{c}^{d}{f(a, y)} dy + 4\int\limits_{c}^{d}{f(\frac{{a + b}}{2}, y) dy + \int\limits_{c}^{d}{f(b, y) dy}} }\right] . $$$Теперь применим формулу Симпсона для каждого из трех полученных интервалов:
$$\begin{gather*} I_1 = \int\limits_{c}^{d} f(a, y) dy \approx \frac{\tau_y}{6} \left[f(a, c) + 4f(a, \frac{c + d}{2}) + f(a, d)\right], \\ I_2 = \int\limits_{c}^{d} f(\frac{a + b}{2}, y) dy \approx \frac{\tau_y}{6}\left[f(\frac{a + b}{2}, c) + 4f(\frac{a + b}{2}, \frac{c + d}{2}) + f(\frac{a + b}{2}, d)\right], \\ I_3 = \int\limits_{c}^{d} f(b, y) dy \approx \frac{\tau_y}{6} \left[f(b, c) + 4f(b, \frac{c + d}{2}) + f(b, d)\right] \end{gather*}$$Подставляя приближенные формулы для вычисления интегралов I1, I2, I3 в формулу для вычисления i, получим
Если область интегрирования не является прямоугольной, то ее можно сделать
подобластью большей по площади прямоугольной области, которую в свою очередь разбить на прямоугольные ячейки. Ячейки в этом случае разделяются на внутренние, к ним применяются приведенные формулы численного интегрирования, и граничные, не прямоугольные, площади которых вычисляются по более сложным алгоритмам. При этом приближенное значение интеграла, например, при использовании формулы средних, можно записать как $$I \approx \sum I^{i}_{\mbox{вн}} + \sum I^{j}_{\mbox{гр}},$$ где $$I^{i}_{\mbox{вн}}$$ — приближенное значение интегралов по внутренним ячейкам, $$I^{j}_{\mbox{гр}}$$ — по граничным, i, j — номера ячеек с площадями Si и Sj.
Поскольку 2N. По этой причине обычно используются полиномы степени от нуля до трех (соответственно, формулы прямоугольников со средней точкой, трапеций, Симпсона, 3/8). Вычисление с их помощью интегралов от функций, обладающих высокой степенью гладкости, например, близким к полиномам высокой степени, представляется нерациональным. В выражение для погрешности этих формул входят первая, вторая или четвертая производные. Погрешность определяется низким порядком производной при высокой степени гладкости
интегрируемой функции. Этих недостатков лишены
Формулировка задачи построения квадратурных формул, поставленная Гауссом, такова.
Для заданного количества точек, а именно, для (N + 1) точки, найти такое расположение узлов и такие веса ci, чтобы квадратурная формула
была точной для полиномов как можно более высокой степени, т.е. чтобы rN(t) = 0.
Пояснение. Для некоторых классов функций существуют квадратурные формулы с rN(t) = 0, которые называются точными. Примером такого класса функций являются полиномы
на отрезке [a, b]. Определим на этом отрезке узлы ti, i = 1, … , N и веса ci так, что
Представим PN(t) в виде интерполяционного полинома
при этом остаточный член интерполяции полинома равен нулю: $$P_N^{(N + 1)} (t) = 0.$$
Тогда из предыдущего условия следует
$$$ c_i = \int\limits_{a}^{b}{{\mathop \Pi\limits_{\substack{k \ne i \\ k = 0 }}^{N}} \frac{(t - t_k)}{(t_i - t_k)}} dt,$$$где ci являются базисными функциями полиномов Лагранжа. Квадратурная формула
является точной для любого полинома степени N. Оказывается, эта
формула может быть точной и для полиномов более высокой степени, а именно, 2N + 1, что используется при построении
Пусть формула численного интегрирования имеет вид
$$I = \int\limits_{a}^{b}{f(t)} dt = c_0 f(t_0 ) + c_1 f(t_1 ) + \ldots + c_n f(t_N) + r_N,$$где ci — веса, rN — остаточный член квадратуры.
Положим, что существует многочлен PM(t) степени M > N, для которого квадратурная формула точна, т.е. rN = 0 при f(t) = PM(t):
f(t) = PM(t) = a0 + a1t + a2t2 + ... + aMtM,
где ai — коэффициенты. В этом случае получим
Приравняем выражения в обеих частях равенства при aj:
Получается нелинейная система из M + 1 уравнения с 2(N +
1) неизвестными ci, ti. Отсюда следует, что максимальное значение M есть 2N + 1. Решение этой системы или исследование на его существование и единственность в общем случае затруднительны. Ниже будет рассмотрен пример получения
Гаусс решил эту задачу более простым (в смысле реализации, но не решения!) способом, доказав следующую теорему. Приведем ее без доказательства.
Теорема. Если в качестве узлов ti, i = 0, ..., N в квадратурной формуле используются нули полиномов Лежандра qN + 1(t), а веса ci вычисляются по формулам
то квадратурная формула
$$\int\limits_{- 1}^1 {f(t)dt = \sum\limits_{i = 0}^{N}{c_i f(t_i)} + r_N(t)}$$точна для полиномов степени 2N + 1.
Напомним, что полиномы Лежандра образуют ортогональную систему
функций на отрезке [- 1; 1]:
Первые несколько полиномов Лежандра будут $$$ q_0(t) = 1, q_1 (t) = t, q_2 (t) = \frac{1}{3}(3t^2 - 1), q_4 (t) = \frac{1}{35}(35t ^4 - 30t^2 + 3), \ldots $$$ рекуррентная и общая формулы имеют вид:
$$\begin{gather*} (n + 1)q_{n + 1} (t) = (2n + 1) t q_n (t) - nq_{n - 1} (t), \\ q_n (t) = \frac{1}{{2^{n} (n!)}}\frac{{d^{n}}}{{dt^{n}}}(t^2 - 1)^{n} \end{gather*}$$Заметим, что рекуррентные формулы, связывающие три полинома порядка n
- 1, n и n + 1 уже встречались для полиномов Чебышева. Такие рекуррентные формулы существуют для всех систем ортогональных полиномов.
Погрешность N. Здесь $$$ \alpha_N = \frac{[(N + 1)!]^4}{\{[2(N + )]!\}^3
[2(N + 1) + 1]} $.$$
Пусть требуется вычислить несобственный интеграл
$$\int\limits_{a}^{b}{f(t)dt}$$от функции, обращающейся в бесконечность в некоторой точке $$c \in \left[{a, b}\right].$$ В этом случае интеграл обычно разбивают на два
$$\int\limits_{a}^{b}{f(t)dt = \lim\limits_{\substack{\delta_1 \to 0 \\ \delta_2 \to 0 }} \left\{{\int\limits_{a}^{c - \delta_1 }{f(t)dt + \int\limits_{c + \delta_2 }^{b}{f(t)dt}} }\right\}}.$$Числа $$\delta_1$$ и $$\delta_2$$ выбирают малыми величинами так, чтобы выполнялась оценка:
$$$ \left|{\int\limits_{c - \delta_1 }^{c + \delta_2 }{f(t)dt}}\right| < \frac{1}{2}\varepsilon , $$$где $$\varepsilon$$ — заданное малое положительное число (точность вычисления интеграла). После этого по квадратурным формулам вычисляют определенные интегралы $$I_1 = \int\limits_{a}^{c - \delta_1} f(t)dt$$ и $$I_2 = \int\limits_{c + \delta_2 }^{b}{f(t) dt}$$ с точностью $$\varepsilon /4$$ каждый. После таких вычислений за приближенное значение интеграла с особенностью принимают $$\int\limits_{a}^{b}{f(t)dt \approx I_1 + I_2}$$ (с точностью $$\varepsilon$$ ).
Другой способ вычисления интеграла от функции особенностью, называемый методом Канторовича выделения особенностей, состоит в следующем. Представим подынтегральную функцию в виде суммы:
$$\int\limits_{a}^{b}{f(t)dt = \int\limits_{a}^{b}{g(t)dt} + \int\limits_{a}^{b}{\left[{f(t) - g(t)}\right] dt.}}$$При этом g(t) подбирают так, чтобы она была интегрируемой, а
разность [f(t) - g(t)] — ограниченной.
Пример. Пусть необходимо вычислить $$$ {I} = \int\limits_0^1 {\frac{dt} {{\sqrt {t(1 + t^2 )}}}}. $$$
Представим I как сумму двух интегралов I = I1 + I2, где $$$ I_1 = \int\limits_0^1 \frac{dt}{\sqrt{t}}, I_2 = \int\limits_0^1 [\frac{1}{\sqrt{t(1 + t^2)}} - \frac{1}{\sqrt{t}}]dt . $$$
Интеграл I1 вычисляется аналитически, а I2, поскольку подынтегральная функция ограничена, можно вычислить по квадратным формулам.
Аналогично можно поступить и в следующей задаче:
$$$ \int\limits_0^1 {\frac{e^{- t^2 }}{\sqrt{t}}} dt = \int\limits_0^1 {\frac{1 - t^2 }{\sqrt{t}}dt + \int\limits_0^1 {\frac{e^{- t^2 } - 1 + t^2 }{\sqrt{t}}dt}. $$$Интегрирование быстро осциллирующих функций типа $$I = \int\limits_{a}^{b}{f(t)e^{i\omega t}}dt$$ можно проводить, заменив f(t) на
Метод Монте - Карло используется, как правило, для вычисления кратных интегралов. Рассмотрим задачу вычисления интеграла по многомерному кубу:
$$I = \int\limits_0^1 \int\limits_0^1 \ldots \int\limits_0^1 f(t_1, t_2, \ldots , t_n)dt_{1dt2} \ldots dt_n.$$Для его вычисления можно построить I на
Проблема вычисления подобных интегралов заключается в том, что при росте размерности задачи объем вычисления значительно увеличивается, а задача численного интегрирования превращается из довольно простой в одну из самых сложных и трудоемких. По этой причине приведенные выше квадратурные формулы используются обычно для решения одно - , дву - и трехмерных задач.
Для вычисления интегралов по гиперкубу высокой размерности обычно используется метод Монте - Карло. Суть его состоит в том, что генерируется последовательность случайных точек единичного n - мерного куба $$t_1, t_2, \ldots , t_n \in R^{n}$$ ; очевидно, что чем больше точек участвует в вычислительном процессе, тем больше точность расчета.
Пусть теперь необходимо взять интеграл по области $$\Omega,$$
принадлежащей n - мерному кубу, причем, $$\Omega$$ выделяется неравенствами
Далее генерируется последовательность случайных чисел, равномерно распределенная в единичном гиперкубе, и для всех точек проверяются неравенства $$g_j (t_k) \le 0.$$ Если они выполнены, т.е. $$t_k \in \Omega,$$ то вычисляются значения f(t_k), прибавляющиеся к сумме.
Пусть вычислено M точек, из которых {K} попали в $$\Omega$$ и накоплена сумма $$\sum\limits_{k = 1}^{K}{f(t_k)}.$$
Среднее по объему $$\Omega$$ значение функции f
вычисляется по формуле $$$ \bar{f} = \frac{I_{\Omega}}{S_{\Omega}} $,$$ где $$S_{\Omega } = K/M,$$ $$I_{\Omega }$$ — кратный интеграл по $$\Omega.$$
С другой стороны, это же значение можно приближенно вычислить как сумму:
$$\left[{\sum {f(t_k)}}\right]/K.$$Приравнивая эти выражения, получим:
$$$ \frac{M}{K}I_\Omega \approx \frac{1}{K}\sum\limits_{k = 1}^{K}{f(t_k)}, \mbox{откуда } I_\Omega \approx \frac{1}{M}\sum\limits_{k = 1}^{K}{f(t_k)} . $$$Решение. Локальную оценку погрешности получим, используя формулы для истинного члена интерполяционного полинома:
$$$ \varepsilon_N = \left|{\int\limits_{t_n}^{t_{n + 1}}{R_N (t)dt}}\right| \le \int\limits_{t_n}^{t_{n + 1}}{\frac{{\max \left|{f^{(N + 1)} (\xi )}\right|}}{{(N + 1)!}}} {\tau}+ \max \left|{\prod\limits_{n = 0}^{N}{(t - t_n)}}\right|. $$$В случае формулы трапеции имеем
$$$ \varepsilon_N \le \int\limits_{t_n}^{t_{n + 1}}{\frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{2}\max \left|{(t - t_n)(t - t_{n + 1})}\right|dt} = \frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{{12}}{\tau}^3 . $$$Погрешность для всего отрезка [a, b] будет следующей $$(\tau = (b - a)(N, t_{m} = n\tau )$$:
Решение. В случае двух узлов N = 1 (количество отрезков разбиения), M = 2, N + 1 = 3 (степень полинома). Узлы tn и веса Cn должны удовлетворять следующей системе уравнений:
В данном случае система уравнений будет:
$$\begin{gather*} c_0 + c_1 = \int\limits_{- 1}^1 {1 dt} = 2, \\ c_0 t + c_1 t_1 = \int\limits_{- 1}^1 {t dt} = 0, \\ c_0 t_0^2 + c_1 t_1^2 = \int\limits_{- 1}^1 {t^2 dt} = \frac{2}{3}, \\ c_0 t_0^3 + c_1 t_1^3 = \int\limits_{- 1}^1 {t^3 dt} = 0, \end{gather*} $$откуда получим
$$c_0 = c_1 = 1\quad t_0 = t_1 = \frac{1}{\sqrt{3}} . $$$Эта формула будет точной для полиномов третьей степени.
$$$ I = \int\limits_0^1 {\frac{1}{\sqrt{x}}e^{- x^2 } dx.} $$$
Решение. Представим интеграл I как сумму двух интегралов
Первый интеграл I1 вычисляется аналитически: I1 = 1, 6. Поскольку подынтегральная функция в I2 трижды непрерывно дифференцируема и ограничена, то I2 можно вычислить, например, по формуле прямоугольников с центральной точкой:
h и h/2 ( правило Рунге ).Решение. Интеграл I может быть представлен в виде
где Ip(h) — приближенные значения интеграла,
вычисленные по формуле с порядком точности p с шагом h, $$$ I^{p} \left({\frac{h}{2}}\right) $$$ - значение интеграла, вычисленное по той же формуле с шагом вдвое меньшим. При малых h
константы C и C1 близки. Этот факт тоже необходимо доказывать. Доказательство труда не представляет — пользуясь теоремой Лагранжа о среднем, легко получить, что эти величины отличаются на O(hp). Тогда получим
откуда следует
$$$ C \approx C_1 \approx \frac{{I^{p} \left({\frac{h}{2}}\right) - I^{p}(h)}}{{h^{p} - 2\left({\frac{h}{2}}\right)^{p}}}. $$$Подставив C во вторую формулу для вычисления I(c, h/2), получим:
Во - первых, эта простая формула позволяет относительно дешевым способом
уточнить вычислительное значение интеграла с шагом h/2. Во - вторых, получаем возможность контролировать точность численного интегрирования путем вычисления значения интеграла дважды (с шагами h и h/2 ).
Примечание. Легко получается аналог правила Рунге при вычислении интеграла
для табличной функции. Необходимо лишь с использованием одних и тех же квадратурных формул вычислить интеграл с шагом таблицы h и затем повторить вычисления, выкинув половину точек, с шагом 2h.
$$$ \int\limits_0^1 {\frac{\sin 100x}{1 + x}dx}, \int\limits_1^2\cos 100 x\ln xdx} $$$
(формулы Филона. Подробнее о них в [7.3]).
$$\begin{gather*} \int\limits_0^{1, 5}{\frac{e^{x}}{x^2 }dx}, \int\limits_0^1 {\frac{\arctg x}{x^2 }dx}, \int\limits_0^1 {\frac{\sqrt{x^3 + 1}}{\sqrt {x}}dx}, \\ \int\limits_0^1 {\frac{\cos x - 1}{x^2 }dx}, \int\limits_0^1 {\frac{\cos x}{\sqrt{x}}dx}, \int\limits_0^1 {\frac{x\sin x}{\ln (1 + x)}dx}, \\ \int\limits_1^\infty {\frac{1 - \cos x}{x\sqrt {x}}dx, }\int\limits_1^0 {\frac{\sin x}{x}dx}, \int\limits_0^1 {\frac{\ln (1 + x)}{x}dx} . \end{gather*}$$
$$$ \pi = \int\limits_0^1 {\frac{4}{1 + x^2 }dx, } $$$
используя формулы прямоугольников с центральной точкой, трапеций, Симпсона.
$$\int\limits_0^1 {\sin x^2 dx}$$
по формуле прямоугольников с центральной точкой для достижения точности $$\varepsilon = 10^{ - 4}.$$ Тот же вопрос для подсчета интеграла
$$\int\limits_0^1 {\exp (x^2 )} dx.$$$$$ I = \int\limits_0^1 {\frac{dx}{1 + x^2 }} $$$
по формуле прямоугольников с центральной точкой, трапеций, Симпсона.
[xi, xi - 1] подынтегральная функция аппроксимируется кубическим сплайном:$$$ S_3^{(n)} = a_n + b_n (x - x_n) + \frac{c_n}{2}(x - x_n)^2 + \frac{d_n}{6}(x - x_n)^3 . $$$
Вывести соответствующую формулу сплайн - квадратуры и исследовать ее точность.
В данной лекции будет рассматриваться задача численного интегрирования. Формулы
численного интегрирования функций одного переменного называют квадратурными формулами. Задача приближенного вычисления определенного интеграла (на отрезке или по многомерной области) фактически разбивается на две самостоятельные подзадачи. Первая — это интегрирование
Вторая задача — подсчет значения определенного интеграла от известной
функции. При этом самая ресурсоемкая операция с точки зрения вычислений — подсчет значения функции. Желательно построить численный метод, позволяющий получать как можно более высокую точность при наименьшем количестве
вычислений, при этом выбор узлов квадратурных формул целиком в руках вычислителя. В этом случае наиболее эффективными окажутся
Простейшую квадратурную формулу (формулу численного интегрирования) можно получить следующим образом. Пусть необходимо вычислить интеграл
$$I = \int\limits_{a}^{b}{f(t)}dt.$$Положим, что f(t) на рассматриваемом отрезке [a, b] не изменяется ( $$f(t) \approx const$$ ). Тогда $$I = f(\xi )(b - a), \xi \in \left[{a, b}\right].$$ Если $$$ \xi = \frac{a + b}{2} $,$$ то получим формулу прямоугольников с центральной точкой
Конечно, для константы приведенная выше формула точна — говорят, что
построенная квадратурная формула будет точна на полиномах степени 0. Легко можно доказать, что формула прямоугольников с центральной точкой будет давать точное значение и в случае линейной функции. Для всех других функций эту формулу будем рассматривать как приближенную.
Если предположить, что функция f(t) на отрезке интегрирования [a, b] достаточно близка к линейной, то можно заменить приближенное значение интеграла i площадью трапеции с высотой (b - a) и основаниями f(a) и f(b). Тогда получается формула трапеций
Введем на отрезке интегрирования сетку, определим значения функции в узлах сетки. Узлы в дальнейшем будем именовать узлами квадратурной формулы (или квадратуры). Пусть, как и в задаче интерполяции, имеется совокупность узлов $$\left\{{t_n}\right\}_{n = 0}^{N}, t_n = a + n{\tau}, {\tau}= (b - a)/{N}, t \in \left[{a, b}\right].$$ Пусть также задана таблица $$f_n = \left\{{f(t_n)}\right\}_{n = 0}^{N}.$$ Отрезок [tk, tk + 1] далее иногда будем называть элементарным отрезком.
Заменим подынтегральную функцию ее интерполяционным полиномом в форме Лагранжа. Будем полагать, что
$$\int\limits_{a}^{b}{f(t)} dt \approx \int\limits_{a}^{b}{L_n} (t)dt.$$Рассмотрим некоторые частные случаи.
Формула трапеций. На отрезке [tk, tk + 1] проводим замену подынтегральной функции интерполяционным полиномом первой степени:
после чего, выполнив интегрирование по элементарному отрезку, получим приближенное значение интеграла на [tk, tk + 1]:
После суммирования интегралов по всем элементарным отрезкам [tk, tk + 1] получаем формулу трапеций для отрезка [a, b]:
На равномерной сетке (сетке с равноотстоящими узлами) при $$\tau _{k} = \tau = (b - a)/N$$ полученная формула принимает вид
$$$ I \approx \frac{{\tau}}{2}\sum\limits_{k = 0}^{N - 1}{\left[{f(t_k) + f(t_{k + 1})}\right]} = \frac{{\tau}}{2}\left[{f(t_0 ) + 2f(t_1 ) + \ldots + 2f(t_{N - 1}) + f(t_N)}\right]. $$$Формула Симпсона. Заменим подынтегральную функцию f(t) на отрезке [tk - 1, tk] интерполяционным полиномом (в форме Лагранжа) второй степени. Для простоты положим $$\tau = (b - a)/N = t_{k} - t_{k - 1} = const$$ для всех k — сетка на отрезке интегрирования равномерная. Тогда
После вычисления интеграла от полинома получим приближенное значение интеграла по элементарному отрезку
$$$ I_k \approx \frac{{\tau}}{6}\left[{f_{k - 1} + 4f_{k - 1/2} + f_k}\right]. $$$Суммируя по всем элементарным отрезкам [tk - 1, tk], получим
где $$$ f_k = f(t_k), f_{k + 1/2} = f(\frac{{t_k + t_{k + 1}}}{2}) $.$$ Формулу Симпсона можно также записать, не используя дробных индексов:
$$$ I \approx \frac{{\tau}}{6}(f_0 + 4f_1 + 2f_2 + 4f_3 + \ldots + 2f_{N - 2} + 4f_{N - 1} + f_N), $$$если локальную формулу получать путем интегрирования интерполяционного
полинома второй степени по отрезку [tk - 1, tk + 1]:
где F(t) — [tk - 1, tk + 1] по точкам tk - 1, tk, tk + 1. В этом случае N — число разбиений отрезка на элементарные отрезки — должно быть четным.
Еще одна используемая на практике
Конечно, существуют и формулы интерполяционного типа более высоких порядков. Они не применяются на практике по следующим обстоятельствам. Любая формула интерполяционного типа записывается в виде
$$I \approx \sum\limits_{k = 0}^{M}{\alpha_k f_k}.$$Во всех приведенных выше формулах коэффициенты $$\alpha_k$$ были
положительными. Такие квадратурные формулы называются
Для степени интерполяционного полинома более 7 среди коэффициентов встречаются отрицательные. Д.Пойа показал, что
где $$\alpha_{nk}$$ — веса квадратурной формулы, получающейся
при замене подынтегрального выражения интерполяционным полиномом степени n. Такое увеличение суммы абсолютных значений коэффициентов связано с быстрым ростом
Подробнее о свойствах
Погрешность квадратурных формул может быть оценена, например, с использованием остаточного члена интерполяционного полинома:
$$$ \varepsilon_k = \left|\int\limits_{t_k}^{t_{k + 1}}R_N(t) dt\right| \le \frac{\max\limits_{[t_k, t_{k + 1}]}\left|f^{(N + 1)}(\xi)\right| }{(N + 1)!} {\tau}\max\limits_{[t_k, t_{k + 1}]}\left|{\mathop \Pi\limits_{k = 0}^{N}(t - t_k)}\right|. $$$Из последней формулы следует, что квадратурная формула точна, если
подынтегральная функция является многочленом степени не выше N.
Получим, например, локальную оценку погрешности для формулы трапеций, используя формулу для остаточного члена интерполяционного полинома 1 - го порядка:
$$$ \varepsilon_k \le \int\limits_{t_k}^{t_{k + 1}}{\frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{2}} \max \left|{(t - t_k)(t - t_{k + 1})}\right|dt = \frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{{12}}\tau^3 . $$$Тогда погрешность по всему отрезку [a, b] будет составлять
Другой способ получения погрешности квадратурных формул состоит в следующем. Рассмотрим интеграл по элементарному отрезку
$$\begin{gather*} I_k = \int\limits_{t_k}^{t_{k + 1}}{f(t)}dt, \\ f(t) = f(z) + f^{\prime}(z)(t - z) + \frac{{f^{\prime\prime}}}{2}(t - z)^2 + \ldots + R(t - z), \end{gather*}$$где $$z \in \left[{t_k , t_{k + 1}}\right]$$ — некая опорная точка, тогда для приближенного значения интеграла верно
$$I_k = f(z)\tau_k + \xi\tau^2_k + \eta\tau^3_k + \ldots + \int\limits_{t_k}^{t_{k + 1}} R(t - z)dz, \tau_k = t_{k + 1} - t_k.$$Коэффициенты $$\xi , \eta , ...$$ зависят от производных f'(z), f''(z), ....
С другой стороны, любая из рассмотренных квадратурных формул представима в виде
$$$ \bar {I_k} = \tau_k (af_k + bf_{k + 1/2} + cf_{k + 1}). $$$Заменяя в этой формуле значения функций f в точках fk, fk + 1/2, fk + 1 ее разложением по формуле Тейлора, получим
где $$z \in [t_k, t_{k + 1} ].$$
Сравнивая разложения для $$I_k, \bar {I_k},$$ легко заметить, что вместе с первым слагаемым совпадают и другие слагаемые до (m - 1) - го порядка, так что $$\xi = \xi _{1}, \eta = \eta _{1}, ...$$
Разность же несовпадающих слагаемых будет, очевидно, оценкой погрешности
квадратурной формулы на интервале $$[t_k, t_{k + 1} ]: \varepsilon_k =
\left|I_k - \bar {I_k}\right| \le v \max\limits_{[t_k, t_{k + 1}]}\left|f^ {(m)}\right| \tau^{m + 1}_k,$$ где v — константа.
Если просуммировать локальные погрешности по всем интервалам [ tk, tk + 1 ], то получим оценку погрешности квадратурной формулы по всему отрезку [a, b]:
где $${\tau}= \max\limits_k \tau_k$$ на неравномерной сетке, или $$\tau = (b - a)/N$$ на равномерной. Число m называется порядком точности квадратуры.
Получим теперь погрешность формулы прямоугольников (со средней точкой) для $$\tau _{k} = \tau = const$$:
$$$ \varepsilon \le \frac{b - a}{24} \max\limits_{[a, b ]}\left|{f^{\prime\prime}_{t} (t)}\right|{\tau}^2, $$$погрешность формулы трапеций
$$$ \varepsilon \le \frac{b - a}{12} \max\limits_{[a, b]}\left|{f^{\prime\prime}_{t} (t)}\right|{\tau}^2, $$$погрешность формул Симпсона (с дробными и без дробных индексов соответственно)
$$\begin{multiple} \varepsilon \le \frac{b - a}{2880}\max\limits_{[a, b]} \left|{f_{t}^{(4)} (t)}\right|{\tau}^4, \\ \varepsilon \le \frac{{b - a}}{{180}}\max\limits_{[a, b]} \left|{f_{t}^{(4)} (t)}\right|{\tau}^4 \end{multiple)$$Заметим, что если функция f(t) имеет только три непрерывных
производных, то оценка погрешности формулы Симпсона ухудшается на порядок:
Рассмотрим, для примера двукратный интеграл по прямоугольной области $$\Omega (a\le x\le b, c\le y\le d).$$ Аналогично одномерному случаю, в соответствии с формулой прямоугольников (со средней точкой) заменим функцию
ее значением в точке пересечения диагоналей прямоугольника. В таком случае получим $$$ I = \int\limits_{a}^{b} \int\limits_{c}^{d}f(x, y) dx dy \approx Sf (\frac{a + b}{2}, \frac{c + d}{2}) $,$$ где $$S = \tau _{x}\tau _{y}, \tau _{x} = (b - a), \tau _{y} = (d - c).$$ Если разбить прямоугольник на прямоугольные ячейки Si, то получим аналог формулы прямоугольников со средней точкой в виде суммы по i: $$I \approx \sum\limits_i S_if(x_i, y_i),$$ где Si — площадь i - ой прямоугольной ячейки, xi, yi — координаты точки пересечения ее диагоналей.
Применим теперь формулу Симпсона для вычисления двукратного интеграла путем редукции к методу вычисления одномерного интеграла, зачем представим двойной интеграл как
$$I = \int\limits_{a}^{b}{\int\limits_{c}^{d} {f(x, y) dx dy = \int\limits_{a}^{b}{dx}} } \int\limits_{c}^{d}{f(x, y)} dy.$$Сначала применим формулу Симпсона для вычисления внешнего интеграла:
$$$ I \approx \frac{{\tau_x}}{6}\left[{\int\limits_{c}^{d}{f(a, y)} dy + 4\int\limits_{c}^{d}{f(\frac{{a + b}}{2}, y) dy + \int\limits_{c}^{d}{f(b, y) dy}} }\right] . $$$Теперь применим формулу Симпсона для каждого из трех полученных интервалов:
$$\begin{gather*} I_1 = \int\limits_{c}^{d} f(a, y) dy \approx \frac{\tau_y}{6} \left[f(a, c) + 4f(a, \frac{c + d}{2}) + f(a, d)\right], \\ I_2 = \int\limits_{c}^{d} f(\frac{a + b}{2}, y) dy \approx \frac{\tau_y}{6}\left[f(\frac{a + b}{2}, c) + 4f(\frac{a + b}{2}, \frac{c + d}{2}) + f(\frac{a + b}{2}, d)\right], \\ I_3 = \int\limits_{c}^{d} f(b, y) dy \approx \frac{\tau_y}{6} \left[f(b, c) + 4f(b, \frac{c + d}{2}) + f(b, d)\right] \end{gather*}$$Подставляя приближенные формулы для вычисления интегралов I1, I2, I3 в формулу для вычисления i, получим
Если область интегрирования не является прямоугольной, то ее можно сделать
подобластью большей по площади прямоугольной области, которую в свою очередь разбить на прямоугольные ячейки. Ячейки в этом случае разделяются на внутренние, к ним применяются приведенные формулы численного интегрирования, и граничные, не прямоугольные, площади которых вычисляются по более сложным алгоритмам. При этом приближенное значение интеграла, например, при использовании формулы средних, можно записать как $$I \approx \sum I^{i}_{\mbox{вн}} + \sum I^{j}_{\mbox{гр}},$$ где $$I^{i}_{\mbox{вн}}$$ — приближенное значение интегралов по внутренним ячейкам, $$I^{j}_{\mbox{гр}}$$ — по граничным, i, j — номера ячеек с площадями Si и Sj.
Поскольку 2N. По этой причине обычно используются полиномы степени от нуля до трех (соответственно, формулы прямоугольников со средней точкой, трапеций, Симпсона, 3/8). Вычисление с их помощью интегралов от функций, обладающих высокой степенью гладкости, например, близким к полиномам высокой степени, представляется нерациональным. В выражение для погрешности этих формул входят первая, вторая или четвертая производные. Погрешность определяется низким порядком производной при высокой степени гладкости
интегрируемой функции. Этих недостатков лишены
Формулировка задачи построения квадратурных формул, поставленная Гауссом, такова.
Для заданного количества точек, а именно, для (N + 1) точки, найти такое расположение узлов и такие веса ci, чтобы квадратурная формула
была точной для полиномов как можно более высокой степени, т.е. чтобы rN(t) = 0.
Пояснение. Для некоторых классов функций существуют квадратурные формулы с rN(t) = 0, которые называются точными. Примером такого класса функций являются полиномы
на отрезке [a, b]. Определим на этом отрезке узлы ti, i = 1, … , N и веса ci так, что
Представим PN(t) в виде интерполяционного полинома
при этом остаточный член интерполяции полинома равен нулю: $$P_N^{(N + 1)} (t) = 0.$$
Тогда из предыдущего условия следует
$$$ c_i = \int\limits_{a}^{b}{{\mathop \Pi\limits_{\substack{k \ne i \\ k = 0 }}^{N}} \frac{(t - t_k)}{(t_i - t_k)}} dt,$$$где ci являются базисными функциями полиномов Лагранжа. Квадратурная формула
является точной для любого полинома степени N. Оказывается, эта
формула может быть точной и для полиномов более высокой степени, а именно, 2N + 1, что используется при построении
Пусть формула численного интегрирования имеет вид
$$I = \int\limits_{a}^{b}{f(t)} dt = c_0 f(t_0 ) + c_1 f(t_1 ) + \ldots + c_n f(t_N) + r_N,$$где ci — веса, rN — остаточный член квадратуры.
Положим, что существует многочлен PM(t) степени M > N, для которого квадратурная формула точна, т.е. rN = 0 при f(t) = PM(t):
f(t) = PM(t) = a0 + a1t + a2t2 + ... + aMtM,
где ai — коэффициенты. В этом случае получим
Приравняем выражения в обеих частях равенства при aj:
Получается нелинейная система из M + 1 уравнения с 2(N +
1) неизвестными ci, ti. Отсюда следует, что максимальное значение M есть 2N + 1. Решение этой системы или исследование на его существование и единственность в общем случае затруднительны. Ниже будет рассмотрен пример получения
Гаусс решил эту задачу более простым (в смысле реализации, но не решения!) способом, доказав следующую теорему. Приведем ее без доказательства.
Теорема. Если в качестве узлов ti, i = 0, ..., N в квадратурной формуле используются нули полиномов Лежандра qN + 1(t), а веса ci вычисляются по формулам
то квадратурная формула
$$\int\limits_{- 1}^1 {f(t)dt = \sum\limits_{i = 0}^{N}{c_i f(t_i)} + r_N(t)}$$точна для полиномов степени 2N + 1.
Напомним, что полиномы Лежандра образуют ортогональную систему
функций на отрезке [- 1; 1]:
Первые несколько полиномов Лежандра будут $$$ q_0(t) = 1, q_1 (t) = t, q_2 (t) = \frac{1}{3}(3t^2 - 1), q_4 (t) = \frac{1}{35}(35t ^4 - 30t^2 + 3), \ldots $$$ рекуррентная и общая формулы имеют вид:
$$\begin{gather*} (n + 1)q_{n + 1} (t) = (2n + 1) t q_n (t) - nq_{n - 1} (t), \\ q_n (t) = \frac{1}{{2^{n} (n!)}}\frac{{d^{n}}}{{dt^{n}}}(t^2 - 1)^{n} \end{gather*}$$Заметим, что рекуррентные формулы, связывающие три полинома порядка n
- 1, n и n + 1 уже встречались для полиномов Чебышева. Такие рекуррентные формулы существуют для всех систем ортогональных полиномов.
Погрешность N. Здесь $$$ \alpha_N = \frac{[(N + 1)!]^4}{\{[2(N + )]!\}^3
[2(N + 1) + 1]} $.$$
Пусть требуется вычислить несобственный интеграл
$$\int\limits_{a}^{b}{f(t)dt}$$от функции, обращающейся в бесконечность в некоторой точке $$c \in \left[{a, b}\right].$$ В этом случае интеграл обычно разбивают на два
$$\int\limits_{a}^{b}{f(t)dt = \lim\limits_{\substack{\delta_1 \to 0 \\ \delta_2 \to 0 }} \left\{{\int\limits_{a}^{c - \delta_1 }{f(t)dt + \int\limits_{c + \delta_2 }^{b}{f(t)dt}} }\right\}}.$$Числа $$\delta_1$$ и $$\delta_2$$ выбирают малыми величинами так, чтобы выполнялась оценка:
$$$ \left|{\int\limits_{c - \delta_1 }^{c + \delta_2 }{f(t)dt}}\right| < \frac{1}{2}\varepsilon , $$$где $$\varepsilon$$ — заданное малое положительное число (точность вычисления интеграла). После этого по квадратурным формулам вычисляют определенные интегралы $$I_1 = \int\limits_{a}^{c - \delta_1} f(t)dt$$ и $$I_2 = \int\limits_{c + \delta_2 }^{b}{f(t) dt}$$ с точностью $$\varepsilon /4$$ каждый. После таких вычислений за приближенное значение интеграла с особенностью принимают $$\int\limits_{a}^{b}{f(t)dt \approx I_1 + I_2}$$ (с точностью $$\varepsilon$$ ).
Другой способ вычисления интеграла от функции особенностью, называемый методом Канторовича выделения особенностей, состоит в следующем. Представим подынтегральную функцию в виде суммы:
$$\int\limits_{a}^{b}{f(t)dt = \int\limits_{a}^{b}{g(t)dt} + \int\limits_{a}^{b}{\left[{f(t) - g(t)}\right] dt.}}$$При этом g(t) подбирают так, чтобы она была интегрируемой, а
разность [f(t) - g(t)] — ограниченной.
Пример. Пусть необходимо вычислить $$$ {I} = \int\limits_0^1 {\frac{dt} {{\sqrt {t(1 + t^2 )}}}}. $$$
Представим I как сумму двух интегралов I = I1 + I2, где $$$ I_1 = \int\limits_0^1 \frac{dt}{\sqrt{t}}, I_2 = \int\limits_0^1 [\frac{1}{\sqrt{t(1 + t^2)}} - \frac{1}{\sqrt{t}}]dt . $$$
Интеграл I1 вычисляется аналитически, а I2, поскольку подынтегральная функция ограничена, можно вычислить по квадратным формулам.
Аналогично можно поступить и в следующей задаче:
$$$ \int\limits_0^1 {\frac{e^{- t^2 }}{\sqrt{t}}} dt = \int\limits_0^1 {\frac{1 - t^2 }{\sqrt{t}}dt + \int\limits_0^1 {\frac{e^{- t^2 } - 1 + t^2 }{\sqrt{t}}dt}. $$$Интегрирование быстро осциллирующих функций типа $$I = \int\limits_{a}^{b}{f(t)e^{i\omega t}}dt$$ можно проводить, заменив f(t) на
Метод Монте - Карло используется, как правило, для вычисления кратных интегралов. Рассмотрим задачу вычисления интеграла по многомерному кубу:
$$I = \int\limits_0^1 \int\limits_0^1 \ldots \int\limits_0^1 f(t_1, t_2, \ldots , t_n)dt_{1dt2} \ldots dt_n.$$Для его вычисления можно построить I на
Проблема вычисления подобных интегралов заключается в том, что при росте размерности задачи объем вычисления значительно увеличивается, а задача численного интегрирования превращается из довольно простой в одну из самых сложных и трудоемких. По этой причине приведенные выше квадратурные формулы используются обычно для решения одно - , дву - и трехмерных задач.
Для вычисления интегралов по гиперкубу высокой размерности обычно используется метод Монте - Карло. Суть его состоит в том, что генерируется последовательность случайных точек единичного n - мерного куба $$t_1, t_2, \ldots , t_n \in R^{n}$$ ; очевидно, что чем больше точек участвует в вычислительном процессе, тем больше точность расчета.
Пусть теперь необходимо взять интеграл по области $$\Omega,$$
принадлежащей n - мерному кубу, причем, $$\Omega$$ выделяется неравенствами
Далее генерируется последовательность случайных чисел, равномерно распределенная в единичном гиперкубе, и для всех точек проверяются неравенства $$g_j (t_k) \le 0.$$ Если они выполнены, т.е. $$t_k \in \Omega,$$ то вычисляются значения f(t_k), прибавляющиеся к сумме.
Пусть вычислено M точек, из которых {K} попали в $$\Omega$$ и накоплена сумма $$\sum\limits_{k = 1}^{K}{f(t_k)}.$$
Среднее по объему $$\Omega$$ значение функции f
вычисляется по формуле $$$ \bar{f} = \frac{I_{\Omega}}{S_{\Omega}} $,$$ где $$S_{\Omega } = K/M,$$ $$I_{\Omega }$$ — кратный интеграл по $$\Omega.$$
С другой стороны, это же значение можно приближенно вычислить как сумму:
$$\left[{\sum {f(t_k)}}\right]/K.$$Приравнивая эти выражения, получим:
$$$ \frac{M}{K}I_\Omega \approx \frac{1}{K}\sum\limits_{k = 1}^{K}{f(t_k)}, \mbox{откуда } I_\Omega \approx \frac{1}{M}\sum\limits_{k = 1}^{K}{f(t_k)} . $$$Решение. Локальную оценку погрешности получим, используя формулы для истинного члена интерполяционного полинома:
$$$ \varepsilon_N = \left|{\int\limits_{t_n}^{t_{n + 1}}{R_N (t)dt}}\right| \le \int\limits_{t_n}^{t_{n + 1}}{\frac{{\max \left|{f^{(N + 1)} (\xi )}\right|}}{{(N + 1)!}}} {\tau}+ \max \left|{\prod\limits_{n = 0}^{N}{(t - t_n)}}\right|. $$$В случае формулы трапеции имеем
$$$ \varepsilon_N \le \int\limits_{t_n}^{t_{n + 1}}{\frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{2}\max \left|{(t - t_n)(t - t_{n + 1})}\right|dt} = \frac{{\max \left|{f^{\prime\prime}(\xi )}\right|}}{{12}}{\tau}^3 . $$$Погрешность для всего отрезка [a, b] будет следующей $$(\tau = (b - a)(N, t_{m} = n\tau )$$:
Решение. В случае двух узлов N = 1 (количество отрезков разбиения), M = 2, N + 1 = 3 (степень полинома). Узлы tn и веса Cn должны удовлетворять следующей системе уравнений:
В данном случае система уравнений будет:
$$\begin{gather*} c_0 + c_1 = \int\limits_{- 1}^1 {1 dt} = 2, \\ c_0 t + c_1 t_1 = \int\limits_{- 1}^1 {t dt} = 0, \\ c_0 t_0^2 + c_1 t_1^2 = \int\limits_{- 1}^1 {t^2 dt} = \frac{2}{3}, \\ c_0 t_0^3 + c_1 t_1^3 = \int\limits_{- 1}^1 {t^3 dt} = 0, \end{gather*} $$откуда получим
$$c_0 = c_1 = 1\quad t_0 = t_1 = \frac{1}{\sqrt{3}} . $$$Эта формула будет точной для полиномов третьей степени.
$$$ I = \int\limits_0^1 {\frac{1}{\sqrt{x}}e^{- x^2 } dx.} $$$
Решение. Представим интеграл I как сумму двух интегралов
Первый интеграл I1 вычисляется аналитически: I1 = 1, 6. Поскольку подынтегральная функция в I2 трижды непрерывно дифференцируема и ограничена, то I2 можно вычислить, например, по формуле прямоугольников с центральной точкой:
h и h/2 ( правило Рунге ).Решение. Интеграл I может быть представлен в виде
где Ip(h) — приближенные значения интеграла,
вычисленные по формуле с порядком точности p с шагом h, $$$ I^{p} \left({\frac{h}{2}}\right) $$$ - значение интеграла, вычисленное по той же формуле с шагом вдвое меньшим. При малых h
константы C и C1 близки. Этот факт тоже необходимо доказывать. Доказательство труда не представляет — пользуясь теоремой Лагранжа о среднем, легко получить, что эти величины отличаются на O(hp). Тогда получим
откуда следует
$$$ C \approx C_1 \approx \frac{{I^{p} \left({\frac{h}{2}}\right) - I^{p}(h)}}{{h^{p} - 2\left({\frac{h}{2}}\right)^{p}}}. $$$Подставив C во вторую формулу для вычисления I(c, h/2), получим:
Во - первых, эта простая формула позволяет относительно дешевым способом
уточнить вычислительное значение интеграла с шагом h/2. Во - вторых, получаем возможность контролировать точность численного интегрирования путем вычисления значения интеграла дважды (с шагами h и h/2 ).
Примечание. Легко получается аналог правила Рунге при вычислении интеграла
для табличной функции. Необходимо лишь с использованием одних и тех же квадратурных формул вычислить интеграл с шагом таблицы h и затем повторить вычисления, выкинув половину точек, с шагом 2h.
$$$ \int\limits_0^1 {\frac{\sin 100x}{1 + x}dx}, \int\limits_1^2\cos 100 x\ln xdx} $$$
(формулы Филона. Подробнее о них в [7.3]).
$$\begin{gather*} \int\limits_0^{1, 5}{\frac{e^{x}}{x^2 }dx}, \int\limits_0^1 {\frac{\arctg x}{x^2 }dx}, \int\limits_0^1 {\frac{\sqrt{x^3 + 1}}{\sqrt {x}}dx}, \\ \int\limits_0^1 {\frac{\cos x - 1}{x^2 }dx}, \int\limits_0^1 {\frac{\cos x}{\sqrt{x}}dx}, \int\limits_0^1 {\frac{x\sin x}{\ln (1 + x)}dx}, \\ \int\limits_1^\infty {\frac{1 - \cos x}{x\sqrt {x}}dx, }\int\limits_1^0 {\frac{\sin x}{x}dx}, \int\limits_0^1 {\frac{\ln (1 + x)}{x}dx} . \end{gather*}$$
$$$ \pi = \int\limits_0^1 {\frac{4}{1 + x^2 }dx, } $$$
используя формулы прямоугольников с центральной точкой, трапеций, Симпсона.
$$\int\limits_0^1 {\sin x^2 dx}$$
по формуле прямоугольников с центральной точкой для достижения точности $$\varepsilon = 10^{ - 4}.$$ Тот же вопрос для подсчета интеграла
$$\int\limits_0^1 {\exp (x^2 )} dx.$$$$$ I = \int\limits_0^1 {\frac{dx}{1 + x^2 }} $$$
по формуле прямоугольников с центральной точкой, трапеций, Симпсона.
[xi, xi - 1] подынтегральная функция аппроксимируется кубическим сплайном:$$$ S_3^{(n)} = a_n + b_n (x - x_n) + \frac{c_n}{2}(x - x_n)^2 + \frac{d_n}{6}(x - x_n)^3 . $$$
Вывести соответствующую формулу сплайн - квадратуры и исследовать ее точность.
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.