5.2.2. Определенный интеграл

Для вычисления определенных интегралов в конечных пределах стандартно используется приближенное представление площади плос-кой области между графиком y = f(x), x  [a, b] и осью абсцисс.

На каждом промежутке [t _i , t _{i + 1}], t _i = a + i  (i = 0, ... , n) фун-

кция представляется постоянной и тогда для h = имеем:

I =  f(x) dx ≈ h  f(a + (i + 0.5)h)

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

0  x ₁ < x ₂ < ... < x _k  1, то можно определить формулу

 f(x) dx ≈  a _i f(x _i) , (3)

где коэффициенты {a _i} выбираются так, чтобы формула была то-чна для многочленов наиболее высокой степени. При k = 1, x ₁ = 0.5,

a ₁ = 1 имеем формулу прямоугольников, при k = 2, x ₁ = 0, x ₂ = 1, a ₁ =

= a ₂ = 0.5 получим формулу трапеций и т. д.

Наиболее популярны так называемые формулы Гаусса типа (3) с оптимальным выбором не только коэффициентов a _i , но и узлов {x _i }. Можно показать, что в качестве {x _i} следует брать корни мно-гочленов Гаусса, ортогональных на [0, 1] (многочлены P_m (x) и P_n (x), m ≠ n, ортогональны на [0, 1], если  P_m (x) P_n (x) dx = 0). Например, первые два такие многочлена суть P₁ (x) = 2х – 1, P₂ (x) = 6х ² - 6х + + 1. Первый имеет корень х = 0.5 (так что приведенная формула пря-моугольников является простейшей формулой Гаусса), второй - два корня x _1,2 = и таким образом формула

 f(x) dx ≈ [f( ) + f( )]

оказывается точной для многочленов 0, 1, 2, 3-й степени. Полагая s = = , t = 1 - s, h = , можно теперь записать

 f(x) dx ≈ [f(a + (i + s)h) + f(a + (i + t)h)] . (4)

Формулы Гаусса - "открытого" типа (не используют значений функ-ции в узлах сетки) и могут применяться к несобственным интегралам и интегралам типа Коши. Например, для  = ln(3) (=1.0986122…)

при а = -1, b = 3, n = 40 по формуле (4) получим значение с восемью

точными знаками (MathCad этот интеграл не считает). Аналогично мо-жно построить формулы вида:

 f(x) e dx ≈  a _i f(x _i) . (5)

Например, при k = 3, х ₁ = 0.19055415, х ₂ = 0.848251867, х ₃ =

= 1.799776578, а ₁ = 0.44602977, а ₂ = 0.396468267, а ₃ = 0.043728888 и для функции f(x) = (2 + x ^-2)ln(1 + x) - по формуле (5) получа-ем 0.99991525 при точном значении 1.

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

5.2.3. Функции на дискретном множестве

Если функция задана конечным числом своих значений, анали-тическое меню Mcad оказывается бесполезным, поскольку классичес-кие операции анализа на таких функциях неопределены. Продолжение таких функций на числовую область называется интерполяцией (вос-полнением) или аппроксимацией (приближением) в зависимости от метода и цели задачи. Простейшим инструментом здесь является ме-тод Лагранжа: если функция задана n своими значениями y ₀, y ₁, ..., y _{n - 1} в точках соответственно x ₀ , x ₁, ..., x _{n - 1}, то многочлен (Лагран-жа)

L _n (t) = y _i

определен для всех t и L _n (x _i) = y _i (i = 0, ..., n - 1). В терминах MathCad этот многочлен можно определить формулой

L(t) y _i  if(i = j, 1, ).

C помощью аппарата разделенных разностей этот многочлен можно записать в ином виде.

Более гибким инструментом для решения задач интерполяции

являются сплайны, представляющие собой агрегат, гладко склеенный

из многочленов фиксированной степени (в данном случае три). Встро-енные функции Lspline, Pspline, Cspline позволяют найти коэффици-енты (с подходящим выбором краевых условий в виде многочленов первой, второй и третьей степени соответственно) по данным x = {x _i}, y = {y _i}. Сам процесс интерполяции осуществляется функцией Interp: сначала полагаем vs cspline(x, y), например, затем можно опреде-лить функцию f(t) interp(vs, x, y, t), интерполирующую массивы данных x, y кубическим сплайном в точке t.

Пример. Пусть функция задана в точках {x _i} = {0, 2, 3.5, 4.5, 7} значениями {y _i} = {0, 0, 1, 0, 1}. Интерполируем ее значения на равно-мерную сетку {z _j = j - 1, j = 0, ..., 9} с помощью многочлена Лагранжа (функция L(t)), сплайна с параболическими краевыми условиями (функция f(t)) и продифференцируем полученные функции. В резуль-тате получилась таблица:

Таблица 4

z	- 1	0	1	2	3	4	5	6	7	8
L(z)	7.426	0	-1.17	0	0.987	0.647	-0.78	-1.67	1	11.61
f(z)	1.288	0	-0.43	0	0.917	0.584	-0.38	-0.26	1	3.406
L′(z)	-12.2	-3.53	0.54	1.385	0.398	-1.04	-1.54	0.273	5.790	16.40
f ′(z)	-1.72	-0.86	0	0.859	0.603	-1.15	-0.46	0.687	1.833	2.979

Можно для той же цели использовать подходящую систему фун-кций { _i(t)} и подобрать в представлении R(t) =  a _i   _i (t) коэффи-циенты a _i наилучшим (в требуемом смысле) образом. Это линейная задача и она легко решается методами п. 2.1.1. Если в данном приме-ре взять r(t) = (t ² + 1) ^{- 1} и положить  _i (t) = r(t - x _i) (трансляции функ-ции r), a коэффициенты a _i (i) выбрать из условий интерполяции: R(x _i) = y _i (i), то в результате получается аналитическая кривая (см. рис. 4). Вектор коэффициентов получается из системы линейных ура-внений Aa = y c матрицей A _{i, j} = r(x _i - x _j ), i, j. Можно положить a = A ^{- 1}y (или a = Lsolve(A,y)) и определить R(t) = a  r(t - x) (скалярное произведение) либо указанной выше суммой.

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

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

(см., например, п. 5.2.6).

Рис. 4

Наконец, на равномерной сетке x _{i + 1} - x _i = h (i) можно пост-роить формулы приближенного дифференцирования в узлах x _i вида y =  a _i y _i , где y _i = y(x _i),  i . Используя формулу Тэйлора, из равенств  a _i  y = y + o(h ^{k + 1}) нетрудно получить для коэффициентов {a _i} линейную систему уравнений a _j (j - m) ⁱ = (i = 0, ..., k) с A _i,j = (j - m) ⁱ , b _i = . Отсюда, например, следуют популярные формулы вычисления второй производной: y₀′′h ² ≈ 2y₀ - 5y ₁ + 4y ₂ - y ₃, y_i′′h ² ≈ y _{i - 1} -2y _i + y _{i + 1}, y _n′′h ² ≈ 2y _n - 5y _{n - 1} + 4y _{n - 2} - y _{n - 3} и др. Приведем таблицу неко-торых часто используемых формул приближенного дифференцирова-

ния (табл. 5 – 6).

Таблица 5

Формулы для вычисления y h∙k! порядка k
m	k = 2	k = 3	k = 4
0	-3y0 +4y1 –y2	-11y0 +18y1 -9y2 +2y3	-25y0+48y1-36y2+16y3-3y4
1	-y0 +y2	-2y0 -3y1 +6y2 -y3	-3y0-10y1+18y2-6y3+y4
2	y0 –4y1 +3y2	y0 -6y1 +3y2 +2y3	y;\s\do2( 0-8y1 +8y3-y4
3		-2y0 +9y1 -18y2 +11y3	-y0+6y1-18y2+10y3+3y4
4			3y0-16y1+36y2-48y3+25y4

Таблица 6

Формулы для вычисления y h2 порядка k
m	k = 2	k = 3	k = 4
0	y;\s\do2( 0-2y1+y2	2y0-5y1+4y2-y3	(35y0-104y1+114y2-56y3+11y4) / 12
1	y;\s\do2( 0-2y1+y2	y;\s\do2( 0-2y1+y2	(11y0-20y1+6y2+4y3-y4) / 12
2	y;\s\do2( 0-2y1+y2	y;\s\do2( 1-2y2+y3	(-y;\s\do2( 0+16y1-30y2+16y3-y4) / 12
3		-y0+4y1-5y2+2y3	(-y0+4y1+6y2-20y3+11y4) / 12
4			(11y0-56y1+114y2-104y3+35y4) / 12