Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3792 / 1087 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Специальности: Программист, Математик
Лекция 8:

Численное интегрирование

Пояснение. Для некоторых классов функций существуют квадратурные формулы с rN(t) = 0, которые называются точными. Примером такого класса функций являются полиномы

P_N (t) = \sum\limits_{k = 0}^{N}{a_k}t^{k}

на отрезке [a, b]. Определим на этом отрезке узлы ti, i = 1, … , N и веса ci так, что

\int\limits_{a}^{b}{P_N (t)dt =  \sum\limits_{i = 0}^{N}{c_i P_N (t_i)}}.

Представим PN(t) в виде интерполяционного полинома

$  
P_N (t) = \sum\limits_{i = 0}^{N}{P_N(t_i) {\mathop \Pi\limits_{\substack{k \ne i \\
k = 0 }}^{N}} \frac{{(t - t_k)}}{{(t_i - t_k)}}},  $

при этом остаточный член интерполяции полинома равен нулю: 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 являются базисными функциями полиномов Лагранжа. Квадратурная формула

\int\limits_{a}^{b}{P_N (t)dt = \sum\limits_{i = 0}^{N}{c_i P_N(t_i)}}

является точной для любого полинома степени 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коэффициенты. В этом случае получим

\begin{gather*}
a_0 \int {dt} + a_1 \int {tdt} + a_2 \int {t^2 dt} + \ldots + a_M \int {t^{M} dt} = \\ 
= c_0 (a_0 + a_1 t_0 + a_2 t_0^2 + \ldots +  a_M t_0^{M}) + \\ 
+ c_1 (a_0 + a_1 t_1 + a_2 t_1^2 + \ldots +  a_M t_1^{M}) + \ldots + \\ 
+ c_N (a_0 + a_1 t_N + a_2 t_N^2 + \ldots + a_M t_N^{M})
\end{gather*}

Приравняем выражения в обеих частях равенства при aj:

\left. \begin{array}{l}
c_0 + c_1 + \ldots + c_N = \int\limits_{a}^{b}{dt = I_0 }, \\ 
c_0 t_0 + c_1 t_1 + \ldots + c_N t_N = \int\limits_{a}^{b}{tdt = I_1 }, \\ 
 \ldots , \\ 
c_0 t_0^{M} + c_1 t_1^{M} + \ldots + c_N t_N^{M} = \int\limits_{a}^{b}{t^{M} dt = I_M .}
\end{array} \right.

Получается нелинейная система из M + 1 уравнения с 2(N + 1) неизвестными ci, ti. Отсюда следует, что максимальное значение M есть 2N + 1. Решение этой системы или исследование на его существование и единственность в общем случае затруднительны. Ниже будет рассмотрен пример получения квадратурной формулы Гаусса таким путем для двух узлов.

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

Теорема. Если в качестве узлов ti, i = 0, ..., N в квадратурной формуле используются нули полиномов Лежандра qN + 1(t), а веса ci вычисляются по формулам

$  c_i = \int\limits_{- 1}^1 {{\mathop \Pi\limits_{\substack{k \ne i \\ 
k = 0}}^{N}} \frac{{(t - t_k)}}{{(t_i - t_k)}}}dt, $

то квадратурная формула

\int\limits_{- 1}^1 {f(t)dt = \sum\limits_{i = 0}^{N}{c_i f(t_i)} + r_N(t)}

точна для полиномов степени 2N + 1.

Напомним, что полиномы Лежандра образуют ортогональную систему функций на отрезке [- 1; 1]:

\int\limits_{- 1}^1 {q_i (t)q_j (t) dt = 0}\quad \mbox{при} \quad i \ne j; \int\limits_{- 1}^1 {q_i (t) q_j (t)dt \ne 0} \quad \mbox{при}\quad  i = j.

Первые несколько полиномов Лежандра будут

$  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 уже встречались для полиномов Чебышева. Такие рекуррентные формулы существуют для всех систем ортогональных полиномов.

Погрешность квадратурной формулы Гаусса на отрезке будет r_N(t) = 2^{2(N + 1) + 1}\alpha_Nf^{(2(N + 1))}(\xi), при этом \xi  \in [- 1, 1]. Для \xi \in [a, b] формула остаточного члена будет r_N (t) = (b - a)^{2(N + 1) + 1} \alpha_N f^{(2(N + 1))} (\xi ), причем коэффициент \alpha_N быстро убывает с ростом N. Здесь

$  \alpha_N = \frac{[(N + 1)!]^4}{\{[2(N + )]!\}^3
[2(N + 1) + 1]}  $.

Формулы Гаусса обеспечивают высокую точность уже при небольшом количестве узлов (от 4 до 10) В этом случае \alpha_2   \approx  5 \cdot 10^{- 7}, \alpha_3 \approx 6 \cdot 10^{- 10}, \alpha_4 \approx 4 \cdot 10^{- 13} . В практических же вычислениях число узлов составляет от нескольких сотен до нескольких тысяч. Отметим также, что веса квадратур Гаусса всегда положительны, что обеспечивает устойчивость алгоритма вычисления сумм \sum\limits_{i = 0}^{N}c_if(t_i) .