авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 10 |

«С.П. Шарый Курс ВЫЧИСЛИТЕЛЬНЫХ МЕТОДОВ Курс ВЫЧИСЛИТЕЛЬНЫХ МЕТОДОВ С. П. Шарый Институт вычислительных технологий СО РАН ...»

-- [ Страница 4 ] --

n f (m+1) ((x)) (x xi )Ni, f (x) Hm (x) = (m + 1)! i= где Ni кратности узлов, m = N0 + N1 +... + Nn 1 степень интер поляционного полинома, а (x) некоторая точка из [a, b], зависящая от x. В случае решения задачи (2.116)–(2.117) справедливо f (4) ((x)) a+b f (x) H3 (x) = · (x a) x (x b).

24 Далее, из того, что формула Симпсона точна для полиномов тре тьей степени, и из условий (2.116)–(2.117) следуют равенства b ba a+b · H3 (a) + 4H H3 (x) dx = + H3 (b) 6 a ba a+b (2.118) · f (a) + 4f = + f (b).

6 Отсюда уже нетрудно вывести выражение для погрешности квадра турной формулы Симпсона:

b ba a+b f (x) dx · f (a) + 4f R(f ) = + f (b) 6 a b в силу (2.118) f (x) H3 (x) dx = a b f (4) ((x)) a+b · (x a) x (x b) dx.

= 24 a 150 2. Численные методы анализа Из него следует оценка b f (4) ((x)) a+b |R(f )| = · (x a) x (x b) dx 24 a b f (4) ((x)) a+b · (x a) x (x b) dx 24 a b M4 a+b (2.119) · (x a) x (x b) dx, 24 a коль скоро подинтегральное выражение (x a) x (a + b)/2 (x b) не меняет знак на интервале интегрирования [a, b]. Здесь обозначено M4 = maxx[a,b] |f (4) (x)|.

Для вычисления фигурирующего в (2.119) интеграла сделаем заме ну переменных a+b t=x, тогда b a+b (x a) x (x b) dx a ba ba 2 ba t t = t+ dt 2 ba ba (b a)2 (b a) t2 t2 dt = =.

4 ba Окончательно M4 (b a) |R(f )|.

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

2.12. Численное интегрирование 2.12г Общие интерполяционные квадратурные формулы Квадратурными формулами интерполяционного типа мы назвали формулы, получающиеся в результате замены подинтегральной функ ции f (x) интерполяционным полиномом Pn (x), который построен по некоторой совокупности простых узлов из интервала интегрирования.

Выпишем этот полином в форме Лагранжа:

n Pn (x) = f (xi ) i (x), i= где n (x) i (x) =, (x xi ) n (xi ) базисные полиномы Лагранжа (стр. 50), а функция n (x) была опре делена в (2.12):

n (x) = (x x0 )(x x1 ) · · · (x xn ).

Интерполяционная квадратурная формула должна получаться из приближёного равенства b b f (x) dx Pn (x) dx a a в результате выполнения интегрирования в правой части. Как след ствие, в представлении (2.111) её весовые коэфициенты имеют вид b b n (x) i = 0, 1,..., n, (2.120) ci = i (x) dx = dx, (x xi ) n (xi ) a a и эти значения весов ci, определяемых по узлам x0, x1,..., xn, явля ются отличительным характеристическим признаком именно интерпо ляционной квадратурной формулы. Если для заданного набора узлов у какой-либо квадратурной формулы весовые коэффициенты равны (2.120), то можно считать, что она построена на основе алгебраиче ской интерполяции подинтегральной функции по этим узлам, взятым с единичной кратностью.

152 2. Численные методы анализа Теорема 2.12.1 Для того, чтобы квадратурная формула (2.111), по строенная по (n + 1) попарно различным узлам, была интерполяци онной, необходимо и достаточно, чтобы её алгебраическая степень точности была не меньшей n.

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

Доказательство. Необходимость условий теоремы очевидна: интер поляционная квадратурная формула на n + 1 узлах, конечно же, точна на полиномах степени n, когда подинтегральная функция заменяется на равный ей интерполянт.

Покажем достаточность: если квадратурная формула (2.111), по строенная по (n + 1) узлу, является точной для любого алгебраическо го полинома степени n, то её весовые коэффициенты вычисляются по формулам (2.120), т. е. она является квадратурной формулой интерпо ляционного типа.

В самом деле, для базисных интерполяционных полиномов i (x) выполнено свойство (2.9) при i = j, 0, i (xj ) = ij = при i = j, 1, и они имеют степень n. Следовательно, применяя рассматриваемую квадратурную формулу для вычисления интеграла от i (x), получим n n b i (x) dx = ck i (xk ) = ck ik = ci, i = 0, 1,..., n.

a k=0 k= Иными словами, имеет место равенство (2.120), что и требовалось до казать.

В частности, если в интерполяционной квадратурной формуле вме сто подинтегральной функции взять полином P0 (x) = x0 = 1, то полу чаем равенство n b ba = 1 dx = ck, a k= 2.12. Численное интегрирование сумма весов квадратурной формулы равна длине интервала инте грирования.

Погрешность интерполяционных квадратурных формул равна b R(f ) = Rn (f, x) dx, a где Rn (f, x) остаточный член алгебраической интерполяции. В §2.2д была получена оценка для Rn (f, x) в виде (2.24) f (n+1) ((x)) · n (x), Rn (f, x) = (n + 1)!

где (x) [a, b], и поэтому b f (n+1) ((x)) n (x) dx.

R(f ) = (n + 1)! a Следовательно, справедлива оценка b Mn+ (2.121) |R(f )| |n (x)| dx, (n + 1)! a где Mn+1 = maxx[a,b] |f (n+1) (x)|. Из неё можно ещё раз заключить, что квадратурная формула интерполяционного типа, построенная по (n + 1) узлам, является точной для любого полинома степени не более n, поскольку тогда Mn+1 = 0.

Отметим, что оценка (2.121) в некоторых случаях может оказаться завышенной, как это имеет место, к примеру, для формулы Симпсона.

2.12д Дальнейшие формулы Ньютона-Котеса В §2.12б и §2.12в простейшие квадратурные формулы Ньютона Котеса формулы прямоугольников и трапеций, формула Симпсона были выведены и исследованы средствами, индивидуальными для каждой отдельной формулы. В этом разделе мы взглянем на формулы Ньютона-Котеса с более общих позиций.

Возьмём на интервале интегрирования [a, b] равноотстоящие друг от друга узлы ba (n) xk = a + kh, k = 0, 1,..., n, h=.

n 154 2. Численные методы анализа Для определения весов формул Ньютона-Котеса необходимо вычис лить величины (2.120), которые мы обозначим для рассматриваемого частного случая как b n (x) (n) (2.122) Ak = dx, k = 0, 1,..., n, (n) (n) x xk n xk a где n (x) = (x a)(x a h) · · · (x a nh).

Сделаем в интеграле (2.122) замену переменных x = a + th. Тогда dx = h dt, n (x) = hn+1 t(t 1)(t 2) · · · (t n), (n) x xk = h(t k), (n) (n) (n) (n) (n) (n) (n) (n) · · · xk xk1 xk xk+1 · · · xk x(n) = xk x n xk n = (1)nk hn k!(n k)!, где считается, что 0! = 1. Окончательно n (1)nk (n) t(t 1) · · · (t k + 1)(t k 1) · · · (t n) dt, Ak =h k!(n k)! k = 0, 1,..., n. Чтобы придать результату не зависящий от интервала интегрирования вид, положим (n) (n) = (b a) Bk, Ak где n (1)nk (n) t(t 1) · · · (t k + 1)(t k 1) · · · (t n) dt.

Bk = n k!(n k)! (n) Теперь уже величины Bk не зависят от h и [a, b]. Они называются коэффициентами Котеса.

2.12. Численное интегрирование К примеру, для n = (t 1)2 (1) B0 = (t 1) dt = =, 2 0 t2 (1) B1 = t dt = =.

2 0 Мы вновь получили веса квадратурной формулы трапеций (2.112). Для случая n = t3 t 1 1 (2) (t 1)(t 2) dt = 3 + 2t B0 = =, 4 4 3 2 0 t 1 1 (2) t B1 = t(t 2) dt = =, 2 2 3 0 t3 t 1 1 (2) t(t 1) dt = B1 = =.

4 4 3 2 0 Полученные коэффициенты соответствуют формуле Симпсона (2.115).

И так далее.

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

В Табл. 2.2, заимствованной из книги [16], приведены коэффициенты Котеса для n 10 (см. также [4, 22, 36, 66]).

(n) Можно видеть, что с ростом n значения коэффициентов Котеса Bk в зависимости от номера k начинают всё сильнее и сильнее осцилли ровать (напоминая в чём-то пример Рунге, стр. 80). Результатом этого является то необычное и противоестественное обстоятельство, что сре ди весов формул Ньютона-Котеса при числе узлов n = 8, 10 и бльших о встречаются отрицательные. Это снижает ценность соответствующих формул, так как при интегрировании знакопостоянных функций может приводить к вычитанию близких чисел и потере точности.

К середине XX века выяснилось, что отмеченный недостаток ти пичен для формул Ньютона-Котеса высоких порядков. Р.О. Кузьмин получил в [49] асимптотические формулы для коэффициентов Коте 156 2. Численные методы анализа Таблица 2.2. Коэффициенты Котеса n 1 2 3 4 5 6 7 8 9 1 1 7 19 41 751 989 2857 k=0 1 6 8 90 288 840 17280 28350 89600 4 3 16 25 9 3577 5838 15741 k=1 1 6 8 45 96 35 17280 28350 89600 1 3 2 25 9 1323 928 1080 28350 k=2 6 8 15 144 280 17280 1 16 25 34 2989 10496 19344 k=3 8 45 144 105 17280 28350 89600 7 25 9 2989 4540 k=4 90 96 280 17280 89600 19 9 1323 10496 5778 k=5 288 35 17280 28350 89600 41 3577 928 k=6 840 17280 89600 751 5838 1080 k=7 17280 28350 89600 989 15741 k=8 28350 2857 k=9 89600 k = 10 са13, из которых следует, что сумма их модулей, т. е.

n (n) |Bk |, k= неограниченно возрастает с ростом n. Отсюда вытекает, во-первых, что погрешности вычислений с формулами Ньютона-Котеса могут быть 13 Помимо оригинальной статьи Р.О. Кузьмина [49] эти асимптотические формулы можно также увидеть в учебнике [18].

2.12. Численное интегрирование сколь угодно велики (см. §2.12а). Во-вторых, так как дополнительно n n b 1 (n) (n) Bk = Ak = 1 dx = 1, ba ba a k=0 k= (n) то при больших n среди коэффициентов Bk должны быть как поло жительные, так и отрицательные. Доказательство упрощённого вари анта этого результата можно найти в [26].

Общую теорию квадратурных формул Ньютона-Котеса вместе с тщательным исследованием их погрешностей читатель может увидеть, к примеру, в книгах [4, 18, 53]. Отметим, что формулы Ньютона-Котеса высоких порядков не очень употребительны. Помимо отмеченной выше численной неустойчивости они проигрывают по точности результатов на одинаковом количестве узлов формулам Гаусса (изучаемым далее в §2.13) и другим квадратурным формулам.

Из популярных квадратурных формул Ньютона-Котеса приведём ещё формулу трёх восьмых, которая получается при замене подин тегральной функции интерполяционным полиномом 3-й степени:

b ba 2a + b a + 2b + f (b). (2.123) f (x) dx · f (a)+ 3f + 3f 8 3 a Её погрешность оценивается как M4 (b a) |R(f )|, где M4 = maxx[a,b] |f (4) (x)|. Порядок точности этой формулы такой же, как и у формулы Симпсона. Вообще, можно показать, что форму лы Ньютона-Котеса с нечётным числом узлов, один из которых при ходится на середину интервала интегрирования, имеют (как формула Симпсона) повышенный порядок точности [1, 4, 18].

2.12е Метод неопределённых коэффициентов Пусть требуется вычислить интеграл b f (x) dx, a 158 2. Численные методы анализа для которого мы ищем приближение в виде n b f (x) dx ck f (xk ), a k= с заданными узлами x0, x1,..., xn. Весовые коэффициенты ck можно найти из условия зануления погрешности этого равенства для какого то достаточно представительного набора несложно интегрируемых функций fi (x), i = 1, 2,.... Каждое отдельное равенство является урав нением на неизвестные c0, c1,..., cn, и потому, выписав достаточное число подобных уравнений, мы сможем определить желаемые веса, т. е.

построить квадратурную формулу. В этом суть метода неопределён ных коэффициентов.

В качестве пробных функций fp (x), p = 1, 2,... часто берут алгеб раические полиномы. Для равномерно расположенных узлов при этом получаются знакомые нам квадратурные формулы Ньютона-Котеса.

Продемонстрируем работу метода неопределённых коэффициентов для тригонометрических полиномов 1, sin(px), cos(px), p = 1, 2,....

2.13 Квадратурные формулы Гаусса 2.13а Задача оптимизации квадратур Параметрами квадратурной формулы (2.111) n b f (x) dx ck f (xk ) a k= являются узлы xk и весовые коэффициенты ck, k = 0, 1,..., n. Однако, строя квадратурные формулы Ньютона-Котеса, в частности, формулу трапеций и формулу Симпсона, мы заранее задавали положение узлов и потом по ним находили веса. Таким образом, возможности общей формулы (2.111) были использованы не в полной мере, поскольку для достижения наилучших результатов можно было бы управлять поло жением её узлов. Только в формуле средних прямоугольников положе ние единственного узла было выбрано из соображений симметрии, и 2.13. Квадратурные формулы Гаусса это привело к существеному улучшению точности. Напомним для при мера, что специальное неравномерное расположение узлов интерполя ции по корням полиномов Чебышёва существенно улучшает точность интерполирования (см. §2.3).

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

Как следствие, сформулированную в начале этого параграфа зада чу оптимизации узлов можно поставить, к примеру, следующим об разом: для заданного фиксированого числа узлов из интервала ин тегрирования нужно построить квадратурную формулу, которая име ет наивысшую алгебраическую степень точности, т. е. является точ ной на полиномах наиболее высокой степени. Нетривиальное решение этой задачи действительно существует, и формулы наивысшей алгеб раической степени точности называются квадратурными формулами Гаусса, поскольку впервые они были рассмотрены в начале XIX века К.Ф. Гауссом.

Далее для удобства мы будем записывать квадратурные формулы Гаусса не в виде (2.111), а как n b (2.124) f (x) dx ck f (xk ), a k= нумеруя узлы с k = 1, а не с нуля. Требование точного равенства для любого полинома степени m в этой формуле в силу её линейности экви валентно тому, что формула является точной для одночленов f (x) = xl, l = 0, 1, 2,..., m, т. е.

n b xl dx = ck xl l = 0, 1, 2,..., m.

k a k= или n ck xl = bl+1 al+1, (2.125) l = 0, 1, 2,..., m.

k l+ k= 160 2. Численные методы анализа Это система из (m + 1) нелинейных уравнений с 2n неизвестными ве личинами c1, c2,..., cn, x1, x2,..., xn. Число уравнений совпадает с числом неизвестных при m + 1 = 2n, т. е. m = 2n 1, и это, вообще говоря, есть максимальное возможное значение для m. При бльших о значениях m система уравнений (2.125) переопределена и в случае об щего положения оказывается неразрешимой.

Сделанное заключение можно обосновать строго.

Предложение 2.13.1 Алгебраическая степень точности квадратур ной формулы, построенной по n узлам, не может превосходить 2n1.

Доказательство. Пусть x1, x2,..., xn узлы квадратурной форму лы. Рассмотрим интегрирование по интервалу [a, b] функции g(x) = (x x1 )(x x2 ) · · · (x xn ), которая является полиномом степени 2n. Если квадратурная формула (2.124) точна для g(x), то n n ck · 0 = 0, ck g(xk ) = k=1 k= тогда как значение интеграла от g(x) явно не равно нулю. Подынте гральная функция g(x) всюду на [a, b] положительна за исключением лишь конечного множества точек узлов x1, x2,..., xn, и поэтому b g(x) dx 0.

a Полученное противоречие показывает, что квадратурная формула (2.124) не является точной для полиномов степени 2n.

Итак, наивысшая алгебраическая степень точности квадратурной формулы, построенной по n узлам, в общем случае может быть равна 2n1. Для двух узлов это 3, при трёх узлах имеем 5, и т. д. Для сравне ния напомним, что алгебраические степени точности формул трапеций и Симпсона, построенных по двум и трём узлам соответственно, равны всего 1 и 3. При возрастании числа узлов этот выигрыш в алгебраиче ской степени точности формул Гаусса, достигаемый за счёт разумного расположения узлов, нарастает.

2.13. Квадратурные формулы Гаусса 2.13б Простейшие квадратуры Гаусса Перейдём к построению квадратурных формул Гаусса. При неболь ших n система уравнений (2.125) для узлов и весов может быть решена с помощью несложных аналитических преобразований.

Пусть n = 1, тогда m = 2n 1 = 1, и система уравнений (2.125) принимает вид c1 = b a, a2 ).

c1 x1 = 2 (b Отсюда c1 = b a, (b a2 ) = 1 (a + b).

x1 = 2c Как легко видеть, получающаяся квадратурная формула это фор мула (средних) прямоугольников b a+b f (x) dx (b a) · f.

a Нам в самом деле известно, что она резко выделяется своей точностью среди родственных квадратурных формул.

Пусть n = 2, тогда алгебраическая степень точности соответствую щей квадратурной формулы равна m = 2n 1 = 3. Система уравнений (2.125) для узлов и весов принимает вид c1 + c2 = b a, c x + c x = 1 (b2 a2 ), 11 22 c1 x2 + c2 x2 = a3 ), 3 (b 1 c1 x1 + c2 x3 = a4 ).

4 (b Пользуясь симметричностью этой системы уравнений, можно искать симметричное решение, в котором c1 = c2. Это даёт c1 = c2 = 1 (b a), и из первого и второго уравнений тогда получаем (2.126) x1 + x2 = a + b.

162 2. Численные методы анализа Следовательно, после возведения в квадрат x2 + 2x1 x2 + x2 = a2 + 2ab + b2. (2.127) 1 В то же время, с учётом найденных значений c1 и c2 из третьего урав нения системы следует x2 + x2 = 3 (b2 + ab + a2 ), 1 и, вычитая это равенство из (2.127), будем иметь x1 x2 = 6 (b2 + 4ba + a2 ).

(2.128) Соотношения (2.126) и (2.128) на основе известной из элементарной алгебры теоремы Виета позволяют сделать вывод, что x1 и x2 являются корнями квадратного уравнения x2 (a + b) x + 1 (b2 + 4ba + a2 ) = 0, так что 1 (2.129) x1,2 = (a + b) ± (b a).

2 Удовлетворение полученными решениями четвёртого уравнения систе мы проверяется прямой подстановкой. Кроме того, поскольку 1 11 · =, 2 2 то x1 и x2 лежат на интервале [a, b] и действительно могут служить узлами квадратурной формулы.

В целом мы вывели квадратурную формулу Гаусса b ba (2.130) · f (x1 ) + f (x2 ), f (x) dx = a где узлы x1 и x2 определяются посредством (2.129).

Пример 2.13.1 Вычислим с помощью полученной выше формулы Гаусса с двумя узлами (2.130) интеграл / cos x dx, 2.13. Квадратурные формулы Гаусса точное значение которого согласно формуле Ньютона-Лейбница равно sin(/2) sin 0 = 1. В соответствии с (2.129) и (2.130) имеем / /2 3 cos x dx · cos + cos + 2 4 62 4 = 0.998473.

Формула Ньютона-Котеса с двумя узлами 0 и /2 формула тра пеций даёт для этого интеграла значение / cos x dx · cos 0 + cos = 0.785398, 2 точность которого весьма низка.

Чтобы получить с формулами Ньютона-Котеса точность вычисле ния рассматриваемого интеграла, сравнимую с той, что даёт форму ла Гаусса, приходится брать больше узлов. Так, формула Симпсона (2.115), использующая три узла 0, /4 и /2, приводит к резуль тату / /2 cos x dx · cos 0 + 4 cos + cos 6 4 = 1 + 2 2 = 1.00228, погрешность которого по порядку величины примерно равна погреш ности ответа по формуле Гаусса (2.130), но всё-таки превосходит её в полтора раза.

С ростом n сложность системы уравнений (2.125) для узлов и весов формул Гаусса быстро нарастает, так что в общем случае остаётся неяс ным, будет ли разрешима эта система (2.125) при любом наперёд за данном n. Будут ли её решения вещественными? Будут ли эти решения принадлежать интервалу [a, b], чтобы служить практически удобными узлами квадратурной формулы?

Получение ответов на поставленные вопросы непосредственно из системы уравнений (2.125) представляется громоздким и малоперспек тивным. Обычно следуют другим путём, расчленяя получившуюся за дачу на отдельные подзадачи 164 2. Численные методы анализа 1) исследования выбора узлов и 2) построения весов по формулам (2.120), путём интегрирования коэффициентов интерполяционного полинома Лагранжа.

В последнем случае мы пользуемся тем фактом, что конструируемая квадратурная формула Гаусса оказывается квадратурной формулой интерполяционного типа. Это прямо следует из Теоремы 2.12.1, коль скоро формула Гаусса, построенная по n узлам, является точной для полиномов степени не менее n 1. Другой возможный способ решения подзадачи 2 нахождение весов из системы уравнений (2.125), которая линейна относительно c1, c2,..., cn.

2.13в Выбор узлов для квадратурных формул Гаусса Теорема 2.13.1 Квадратурная формула (2.124) n b f (x) dx ck f (xk ) a k= имеет алгебраическую степень точности (2n 1) тогда и только тогда, когда (1) она является интерполяционной формулой;

(2) её узлы x1, x2,..., xn суть корни такого полинома (x) = (x x1 )(x x2 ) · · · (x xn ), что b (2.131) (x) q(x) dx = a для любого полинома q(x) степени не выше (n 1).

Выражение b (x) q(x) dx a интеграл от произведения двух функций, уже встречалось нам в §2.10в. Мы могли видеть, что на пространстве L2 [a, b] всех интегри руемых с квадратом функций оно задаёт скалярное произведение, т. е.

2.13. Квадратурные формулы Гаусса симметричную билинейную и положительно определённую форму. По этой причине утверждение Теоремы 2.13.1 часто формулируют так: для того, чтобы квадратурная формула n b f (x) dx ck f (xk ) a k= являлась точной на полиномах степени (2n 1), необходимо и доста точно, чтобы эта формула была интерполяционной, а её узлы x1, x2,..., xn являлись корнями полинома (x), который с единичным весом ортогонален на [a, b] любому полиному степени не выше (n 1).

Доказательство. Необходимость. Пусть рассматриваемая квадратур ная формула имеет алгебраическую степень точности (2n 1), т. е. точ на на полиномах степени (2n 1). Таковым является, в частности, по лином (x) q(x), имеющий степень n + (n 1), если степень q(x) не превосходит (n 1). Тогда справедливо точное равенство n b (x) q(x) dx = ck (xk ) q(xk ) = 0, a k= поскольку все (xk ) = 0. Коль скоро этот результат верен для любого полинома q(x) степени не выше n 1, то отсюда следует выполнение условия (2).

Справедливость условия (1) следует из Теоремы 2.12.1: если постро енная по n узлам квадратурная формула (2.124) является точной для любого полинома степени n 1, то она интерполяционная.

Достаточность. Пусть имеется полином (x) степени n, имеющий n различных корней на интервале [a, b] и удовлетворяющий условию ор тогональности (2.131) с любым полиномом q(x) степени не выше (n1).

Покажем, что интерполяционная квадратурная формула, построенная по узлам x1, x2,..., xn, которые являются корнями (x), будет точна на полиномах степени 2n 1.

Пусть f (x) произвольный полином степени 2n 1. Тогда после деления его на (x) получим представление (2.132) f (x) = (x) q(x) + r(x), где q(x) и r(x) соответственно частное и остаток от деления f (x) на (x). При этом полином q(x) имеет степень (2n 1) n = n 1, а 166 2. Численные методы анализа степень полинома-остатка r(x) по определению не превосходит n 1.

Отсюда b b b b (2.133) f (x) dx = (x) q(x) dx + r(x) dx = r(x) dx a a a a в силу сделанного нами предположения об ортогональности (x) всем полиномам степени не выше n 1.

Но по условиям теоремы рассматриваемая квадратурная форму ла является интерполяционной и построена по n узлам. Поэтому она является точной на полиномах степени n 1 (см. Теорему 2.12.1), в частности, на полиноме r(x). Следовательно, n n b r(x) dx = ck r(xk ) = ck (xk ) q(xk ) + r(xk ) a k=1 k= в силу равенств (xk ) = n поскольку имеет место (2.132).

= ck f (xk ), k= Итак, сравнивая результаты этой выкладки с (2.133), будем иметь n b f (x) dx = ck f (xk ), a k= т. е. исследуемая квадратурная формула действительно является точ ной на полиномах степени 2n 1.

Подведём промежуточные итоги. Процедура построения квадратур ных формул Гаусса разделена нами на две отдельные задачи нахожде ния узлов и вычисления весов. В свою очередь, узлы квадратурной формулы, как выясняется, можно взять корнями некоторых специаль ных полиномов (x), удовлетворяющих условиям Теоремы 2.13.1. В этих полиномах легко угадываются знакомые нам из §2.11 ортогональ ные полиномы, которые являются полиномами Лежандра для случая [a, b] = [1, 1] или соответствующим образом преобразованы из них для произвольного интервала интегрирования [a, b].

2.13. Квадратурные формулы Гаусса 2.13г Практическое применение формул Гаусса Отдельное нахождение узлов и весов формул Гаусса для каждого конкретного интервала интегрирования [a, b] является весьма трудоза тратным, и если бы нам нужно было проделывать эту процедуру вся кий раз при смене интервала [a, b], то практическое применение формул Гаусса значительно потеряло бы свою привлекательность. Естествен ная идея состоит в том, чтобы найти узлы и веса формул Гаусса для какого-то одного канонического интервала, а затем получать их для любого другого интервала с помощью несложных преобразований.

В качестве канонического интервала обычно берут [1, 1], т. е. тот интервал, для которого строятся ортогональные полиномы Лежандра.

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

Если x = 1 (a + b) + 2 (b a) y, (2.134) то переменная x будет пробегать интервал [a, b], когда y изменяется в [1, 1]. Обратное преобразование даётся формулой 2x (a + b).

y= ba В частности, если yi, i = 1, 2,..., n, корни полинома Лежандра, кото рые согласно Предложению 2.11.2 все различны и лежат на интервале [1, 1], то узлы квадратурной формулы Гаусса для интервала интегри рования [a, b] суть xi = 1 (a + b) + 2 (b a) yi, (2.135) i = 1, 2,..., n.

Все они также различны и лежат на интервале интегрирования [a, b].

Далее, веса ck любой интерполяционной квадратурной формулы вы ражаются в виде интегралов (2.120). В случае формул Гаусса (когда узлы нумеруются с единицы) они принимают вид b ck = k (x) dx, k = 1, 2,..., n, a 168 2. Численные методы анализа где k (x) k-ый базисный полином Лагранжа (см. стр. 50), построен ный по узлам (2.135):

(x x0 ) · · · (x xi1 )(x xi+1 ) · · · (x xn ) i (x) =.

(xi x0 ) · · · (xi xi1 )(xi xi+1 ) · · · (xi xn ) Тогда, выполняя замену переменных (2.134), получим + b) + 2 (b a) y = 1 (b a) dy, dx = d 2 (a и потому b a) ck = k (x) dx = 2 (b k (y) dy, k = 1, 2,..., n, a где k (y) k-ый базисный полином Лагранжа, построенный по узлам yi, i = 1, 2,..., n, которые являются корнями n-го полинома Лежандра.

Получается, что веса квадратурной формулы Гаусса для произвольно го интервала интегрирования [a, b] вычисляются простым умножением весов для канонического интервала [1, 1] на множитель 2 (b a) радиус интервала интегрирования.

Для интервала [1, 1] узлы квадратурных формул Гаусса (т. е. кор ни полиномов Лежандра) и их веса тщательно затабулированы для первых натуральных чисел n вплоть до нескольких десятков. Обсуж дение вычислительных формул и других деталей численных процедур для их нахождения читатель может найти, к примеру, в книгах [4, 53].

Конкретные значения узлов и весов квадратур Гаусса приводятся в по дробных руководствах по численным методам [2, 4, 10, 16, 17, 56] или в специальных справочниках, например, в [48, 36]. В частности, в учебни ке [4] значения весов и узлов фомул Гаусса приведены для небольших n с 16 значащими цифрами, в книге [17] с 15 значащими цифрами вплоть до n = 16, а в справочниках [48, 36] с 20 значащими цифрами вплоть до n = 48 и n = 96. Таким образом, практическое применение квадратур Гаусса обычно не встречает затруднений.

При небольших значениях n можно дать точные аналитические вы ражения для узлов формул Гаусса, как корней полиномов Лежандра Ln (x), имеющих явные представления (2.106). Так, для n = 5x3 3x = x 5x2 3.

1 L3 (x) = 2 2.13. Квадратурные формулы Гаусса Таблица 2.3. Узлы и веса квадратурных формул Гаусса Узлы Веса n= ±0.57735 02691 89626 1.00000 00000 n= 0.00000 00000 00000 0.88888 88888 ±0.77459 66692 41483 0.55555 55555 n= ±0.33998 10435 84856 0.65214 51548 ±0.86113 63115 94053 0.34785 48451 n= 0.00000 00000 00000 0.56888 88888 ±0.53846 93101 05683 0.47862 86704 ±0.90617 98459 38664 0.23692 68850 Поэтому для канонического интервала интегрирования [1, 1] и n = узлы квадратурной формулы Гаусса суть x1 = = 0.77459 66692 41483..., x2 = 0, x3 = = 0.77459 66692 41483....

Для n = 35x4 30x2 + 3, L4 (x) = и нахождение корней этого биквадратного полинома труда не пред ставляет. Аналогично и для n = 5, когда 63x5 70x3 + 15x = x 63x4 70x2 + 15.

1 L5 (x) = 8 170 2. Численные методы анализа Численные значения узлов и весов квадратурных формул Гаусса для n = 2, 3, 4, 5 сведены в Табл. 2.3. Видно, что узлы располагают ся симметрично относительно середины интервала интегрирования, а равноотстоящие от неё весовые коэффициенты одинаковы.

2.13д Погрешность квадратур Гаусса Для исследования остаточного члена квадратурных формул Гаус са предположим, что подинтегральная функция f (x) имеет достаточно высокую гладкость. Построим для неё интерполяционный многочлен, принимающий в узлах x1, x2,..., xn значения f (x1 ), f (x2 ),..., f (xn ).

Коль скоро квадратурная формула Гаусса точна на полиномах степени 2n 1, то для адекватного учёта этого факта степень полинома, ин терполирующего подинтегральную функцию, также нужно взять рав ной 2n 1. Необходимую степень можно получить, рассматривая, как и для формулы Симпсона, интерполяцию с кратными узлами, напри мер, назначая кратность всех n узлов равной двум и предполагая в них известными фиктивные значения производных f (x1 ), f (x2 ),..., f (xn ).

При этом согласно (2.44) погрешность интерполирования функции f (x) полиномом Эрмита H2n1 (x) равна R2n1 (f, x) = f (x) H2n1 (x) n f (2n) ((x)) (x xi )2, · = (2n)! i= f (2n) ((x)) · (x), = (2n)!

где (x) = (x x1 )(x x2 ) · · · (x xn ) и (x) некоторая точка, за висящая от x, из открытого интервала интерполирования ] a, b [. По условиям интерполяции H2n1 (xi ) = f (xi ), i = 1, 2,..., n, следователь 2.13. Квадратурные формулы Гаусса но, b b f (x) dx = H2n1 (x) + R2n1 (f, x) dx a a b b = H2n1 (x) dx + R2n1 (f, x) dx a a n b = ci H2n1 (xi ) + R2n1 (f, x) dx a i= n b 1 f (2n) ((x)) (x) = ci f (xi ) + dx, (2n)! a i= где ci веса квадратурной формулы Гаусса.

Выражение для второго слагаемого последней суммы, т. е. для оста точного члена квадратуры, можно упростить и далее, приняв во внима ние знакопостоянство множителя (x). Тогда в силу интегральной теоремы о среднем (см., к примеру, [35]) имеем b b 2 f (2n) ((x)) (x) dx = f (2n) () (x) dx a a для некоторой точки ]a, b[. Таким образом, погрешность квадра турной формулы Гаусса, построенной по n узлам x1, x2,..., xn [a, b], равна f (2n) () b R(f ) = (x) dx, (2n)! a где ] a, b [.

Узлы x1, x2,..., xn это корни полинома, полученного из поли нома Лежандра линейной заменой переменных. По этой причине инте грал в полученной формуле для погрешности можно вычислить точно, основываясь на результате Предложения 2.11.1. Это даёт (n!) (b a)2n+1 f (2n) () (2.136) R(f ) = (2n)! (2n + 1) для некоторой промежуточной точки из интервала интегрирования 172 2. Численные методы анализа ] a, b [. Иногда удобнее грубая оценка (n!) M2n (b a)2n+1, |R(f )| (2n)! (2n + 1) где, как обычно, обозначено Mp = maxx[a,b] |f (p) (x)|. В частности, для формулы Гаусса (2.129)–(2.130) с двумя узлами M4 (b a) |R2 (f )|, что даже лучше оценки погрешности для формулы Симпсона. Практи ческое поведение этой погрешности мы могли видеть в Примере 2.13.1.

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

Другое важное наблюдение состоит в том, что в выражении (2.136) числитель (n!)4 с ростом n может быть сделан сколь угодно меньшим знаменателя, превосходящего ((2n)!)4. Как следствие, квадратурные формулы Гаусса дают пример ненасыщаемого численного метода, по рядок точности которого может быть сделан сколь угодно высоким с ростом числа узлов и гладкости функции, которая подаётся ему на вход. В этом квадратуры Гаусса принципиально отличаются, к при меру, от интерполяции с помощью сплайнов, которая сталкивается с ограничением на порядок сходимости, не зависящим от гладкости ис ходных данных (стр. 91).

В заключение темы следует сказать, что на практике нередко тре буется включение во множество узлов квадратурной формулы каких либо фиксированных точек интервала интегрирования, например, его концов (одного или обоих), либо каких-то выделенных внутренних то чек. Основная идея формул Гаусса может быть применена и в этом случае, что приводит к квадратурам Маркова [4, 12, 53, 56], которые называются также квадратурами Лобатто [2, 36].

Построение квадратурных формул Гаусса основывалось на оптими зации алгебраической степени точности квадратур. Эта идея может быть модифицирована и приспособлена к другим ситуациям, когда ал гебраические полиномы и их степень уже не являются наиболее адек ватным мерилом качества квадратурной формулы. Например, мож но развивать квадратуры наивысшего тригонометрического порядка 2.14. Составные квадратурные формулы точности, которые окажутся практичнее при вычислении интегралов от осциллирующих функций [56].

2.14 Составные квадратурные формулы Рассмотренные выше квадратурные формулы дают приемлемую по грешность в случае, когда длина интервала интегрирования [a, b] неве лика и подинтегральная функция имеет на нём ограниченные произ водные. Но если величина (b a) относительно велика или интегрируе мая функция имеет большие производные, то погрешность вычисления интеграла делается значительной или даже большой. Тогда для полу чения требуемой точности вычисления интеграла применяют состав ные квадратурные формулы, основанные на разбиении интервала инте грирования на подинтервалы меньшей длины, по каждому из которых вычисляется значение элементарной квадратуры, а затем искомый интеграл приближается их суммой. Зафиксируем некоторую квадратурную формулу. Будем также счи тать, что её погрешность R(f ) имеет оценку |R(f )| C(b a)p, где C константа, зависящая от типа квадратурной формулы и инте грируемой функции, но в любом случае p 1. К примеру, для формулы средних прямоугольников C = M2 /24 и p = 3, а для формулы Симпсо на C = M4 /2880 и p = 5. Разобьём интервал интегрирования точками r1, r2,..., rN 1 на N равных частей [a, r1 ], [r1, r2 ],..., [rN 1, b] дли ны h = (b a)/N. Тогда, пользуясь аддитивностью интеграла, можно вычислить по рассматриваемой формуле интегралы ri+ i = 0, 1,..., N 1, f (x) dx, ri где r0 = a, b = rN, а затем положить N b ri+ (2.137) f (x) dx f (x) dx.

a ri i= 14 Интересно, что при этом подинтегральная функция интерполируется на всём интервале интегрирования [a, b] при помощи сплайна.

174 2. Численные методы анализа Рис. 2.24. Составная квадратурная формула трапеций На каждом подинтервале [ri, ri+1 ] p ba = Chp, |R(f )| C N а полная погрешность интегрирования R(f ) при использовании пред ставления (2.137) не превосходит суммы погрешностей отдельных сла гаемых, т. е.

p C(b a)p ba = C(b a)hp1.

|R(f )| CN = N p N Как видим, эта погрешность уменьшилась в N p1 раз, и потенциаль но таким способом погрешность вычисления интеграла можно сделать сколь угодно малой.

Число (p1) часто называют порядком точности (составной) квад ратурной формулы, и это понятие очевидно согласуется с данным ранее Определением 2.8.1. Ясно, что основная идея составных квадратурных формул работает и в случае неравномерного разбиения интервала ин тегрирования на более мелкие части, но анализ погрешности проводить тогда труднее.

Для равномерного разбиения интервала интегрирования составные квадратурные формулы выглядят особенно просто. Выпишем их яв ный вид для рассмотренных выше простейших квадратур Ньютона Котеса и разбиения интервала интегрирования [a, b] на N равных ча стей [r0, r1 ], [r1, r2 ],..., [rN 1, rN ] длины h = (b a)/N каждая, в ко тором a = r0 и rN = b.

2.14. Составные квадратурные формулы Составная формула (средних) прямоугольников N b f (x) dx h f (ri1/2 ), a i= где ri1/2 = ri h/2. Её полная погрешность (b a)h |R(f )| M2, т. е. она имеет второй порядок точности. Эта формула, как нетрудно видеть, совпадает с интегральной суммой Римана для интеграла от f (x) по интервалу [a, b].

Составная формула трапеций N b f (ri ) + 1 f (b).

f (x) dx h 2 f (a) + a i= Её полная погрешность (b a)h |R(f )| M2, т. е. порядок точности тоже второй.

Составная формула Симпсона (парабол) N b h f (x) dx f (ri1 ) + 4f (ri1/2 ) + f (ri ), a i= где ri1/2 = ri h/2. Её полная погрешность (b a)h |R(f )| M4, т. е. формула имеет четвёртый порядок точности.

Аналогично можно сконструировать составные квадратурные фор мулы Гаусса, но мы не будем здесь развёртывать детали этого постро ения.

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

Так же растёт и погрешность суммирования результатов интегрирова ния по отдельным подобластям общей области интегрирования. Поэто му эффект увеличения точности составной формулы при возрастании размерности становится всё менее ощутимым.

2.15 Сходимость квадратур С теоретический точки зрения интересен вопрос о сходимости квад ратур при неограниченном возрастании числа узлов. Иными словами, верно ли, что n b ck f (xk ) f (x) dx a k= при n (здесь узлы и веса квадратурных формул нумеруются с нуля)?

Похожий вопрос вставал при исследовании интерполяционного про цесса, и мы обсуждали его в §2.5. Но в случае квадратурных формул помимо бесконечной треугольной матрицы узлов (1) ··· x (2) (2) ··· x x 0 (2.138), (3) (3) (3) ··· x x1 x 0.......

.....

...

(n) (n) (n) таких что xk лежат на интервале интегрирования [a, b] и xi = xj при i = j, необходимо задавать ещё и треугольную матрицу весовых 2.15. Сходимость квадратур коэффициентов квадратурных формул (1) ··· c (2) (2) ··· c c 0 (2.139).

(3) (3) (3) ··· c c1 c 0.......

.....

...

В случае задания бесконечных треугольных матриц (2.138)–(2.139), по которым организуется приближённое вычисление интегралов на по следовательности сеток, будем говорить, что определён квадратурный процесс.

Определение 2.15.1 Квадратурный процесс, задаваемый зависящим от целочисленного параметра n семейством квадратурных формул n b (n) (n) f (x) dx ck f xk, n = 0, 1, 2,..., a k= которые определяются матрицами узлов и весов (2.138)–(2.139), бу дем называть сходящимся для функции f (x) на интервале [a, b], если n b (n) (n) lim ck f xk = f (x) dx, n a k= т. е. если при неограниченном возрастании числа узлов n предел ре зультатов квадратурных формул равен точному интегралу от функ ции f по [a, b].

Необходимо оговориться, что в практическом плане вопрос о сходи мости квадратур решается положительно с помощью составных фор мул, рассмотренных в предыдущем §2.14. Для достаточно общих функ ций путём построения составной квадратурной формулы всегда мож но добиться сходимости приближённого значения интеграла к точно му (для составной формулы прямоугольников это следует из самого определения интегрируемости по Риману). Обсуждаемый ниже круг вопросов относится больше к теоретическим качествам тех или иных чистых квадратурных формул, их предельному поведению.

Весьма общие достаточные условия для сходимости квадратур были сформулированы и обоснованы В.А. Стекловым [61], а впоследствии Д. Пойа [70] доказал также необходимость условий Стеклова.

178 2. Численные методы анализа Теорема 2.15.1 (теорема Стеклова-Пойа) Квадратурный процесс, по рождаемый матрицами узлов и весов (2.138)–(2.139), сходится для любой непрерывной на [a, b] функции тогда и только тогда, когда (1) этот процесс сходится для полиномов, (2) суммы абсолютных значений весов равномерно по n ограничены, т. е. существует такая константа C, что n (n) (2.140) C ck k= для всех n = 0, 1, 2,....

Покажем достаточность условий теоремы Стеклова-Пойа. С этой целью, задавшись каким-то 0, найдём полином PN (x), который рав номерно с погрешностью приближает непрерывную подинтегральную функцию f (x) на рассматриваемом интервале [a, b]. Существование та кого полинома обеспечивается теоремой Вейерштрасса (см. §2.5). Далее преобразуем выражение для остаточного члена квадратурной форму лы:

n b (n) (n) f (x) dx Rn (f ) = ck f xk a k= n b b (n) (n) f (x) PN (x) dx + PN (x) dx = ck f xk a a k= b f (x) PN (x) dx = a n b (n) (n) PN (x) dx + ck PN xk a k= n (n) (n) (n) f xk + ck PN xk.

k= Отдельные слагаемые полученной суммы, расположенные выше в раз личных строчках, оцениваются при достаточно больших номерах n сле 2.15. Сходимость квадратур дующим образом:

b f (x) PN (x) dx (b a), так как PN (x) приближает f (x) равномерно с погрешностью a на интервале [a, b], n b (n) (n), так как квадратуры PN (x) dx ck PN xk сходятся на полиномах, a k= n n (n) (n) (n) (n) в силу (2.140).

f xk C ck PN xk ck k=0 k= Поэтому в целом, если n достаточно велико, имеем |Rn (x)| (b a + 1 + C).

Это и означает сходимость рассматриваемого квадратурного процесса.

Доказательство необходимости условия теоремы Стеклова-Пойа по мимо оригинальной статьи [70] можно найти также в книгах [4, 26].

В формулировке теоремы фигурирует величина n (2.141) ck, k= сумма абсолютных значений весов, которая, как мы видели в §2.12а, является коэффициентом увеличения погрешности в данных и игра ет очень важную роль при оценке качества различных квадратурных формул. В §2.12д уже упоминался результат Р.О. Кузьмина [49] о том, что для формул Ньютона-Котеса величина (2.141) неограниченно уве личивается с ростом числа узлов n. Как следствие, на произвольных непрерывных функциях эти квадратурные формулы сходимостью не обладают.

Для квадратурных формул Гаусса ситуация другая. Справедливо Предложение 2.15.1 Весовые коэффициенты квадратурных формул Гаусса положительны.

180 2. Численные методы анализа ck 0. 0. 0. k 0 48 Рис. 2.25. Зависимость весовых коэффициентов от номера для квадратуры Гаусса 96-го порядка Доказательство. Ранее мы уже выводили для весов интерполяцион ных квадратурных формул выражение (2.120). Зафиксировав индекс i {1, 2,..., n}, дадим другое явное представление для весового коэф фициента ci квадратурной формулы Гаусса, из которого и будет следо вать доказываемое предложение.

Пусть x1, x2,..., xn узлы квадратурной формулы Гаусса на ин тервале интегрирования [a, b]. Так как формулы Гаусса имеют алгеб раическую степень точности 2n 1, то для полинома i (x) = (x x1 ) · · · (x xi1 )(x xi+1 ) · · · (x xn ) степени 2(n 1) должно выполняться точное равенство n b (2.142) i (x) dx = ck i (xk ).

a k= Но i (xk ) = ik по построению полинома i, так что от суммы справа в (2.142) остаётся лишь одно слагаемое ci i (xi ):

b i (x) dx = ci i (xi ).

a Следовательно, b ci = i (x) dx i (xi ).

a Далее, i (x) 0 всюду на интервале интегрирования [a, b] за ис ключением конечного числа точек, и потому положителен интеграл в 2.16. Вычисление интегралов методом Монте-Карло числителе выписанного выражения. Кроме того, i (xi ) 0, откуда можно заключить, что ci 0.

Кроме того, сумма весов формул Гаусса равна длине интервала ин тегрирования (как и для всех интерполяционных квадратурных фор мул, см. §2.12г). Как следствие, величина (2.141) при этом ограничена, и квадратурный процесс по формулам Гаусса всегда сходится.

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

2.16 Вычисление интегралов методом Монте-Карло В методе Монте-Карло, называемом также методом статисти ческого моделирования, искомое решение задачи представляется в виде какой-то характеристики некоторого случайного процесса. Затем этот процесс каким-либо образом моделируется и по его реализации мы вы числяем нужную характеристику, т. е. решение задачи.

В качестве примера рассмотрим задачу вычисления определённого интеграла b f (x) dx.

a Согласно известной из интегрального исчисления теореме о среднем b f (x) dx = (b a) f (c) a для некоторой точки c [a, b]. Смысл средней точки c можно понять глубже с помощью следующего рассуждения. Пусть интервал интегри рования [a, b] разбит на N равных подинтервалов. По определению ин теграла Римана, если xi точки из этих подинтервалов, то N N b ba f (x) dx f (xi ) = (b a) · f (xi ) N N a i=1 i= для достаточно больших N. Сумма в правой части это произведение ширины интервала интегрирования (b a) на среднее арифметическое 182 2. Численные методы анализа значений подинтегральной функции f в точках xi. Таким образом, ин теграл от f (x) по [a, b] есть не что иное, как среднее значение функ ции f (x) на интервале [a, b], умноженное на ширину этого интервала.

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

Для формализации высказанных идей целесообразно привлечь ап парат теории вероятностей. Эта математическая дисциплина исследует случайные явления, которые подчиняются свойству статистической устойчивости, т. е. обнаруживают закономерности поведения в боль ших сериях повторяющихся испытаний. Одними из основных понятий теории вероятностей являются понятия случайной величины и её функ ции распределения. Случайной величиной называется переменная вели чина, значения которой зависят от случая и для которой определена функция распределения вероятностей. Конкретное значение, которое случайная величина принимает в результате опыта, обычно называют реализацией случайной величины.

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

На этом пути мы и приходим к простейшему методу Монте-Карло для вычисления определённого интеграла:

2.16. Вычисление интегралов методом Монте-Карло Рис. 2.26. Вычисление объёма области методом Монте-Карло.

организуем реализации i, i = 1, 2,..., N, для случайной величины, имеющей равномерное распределение на интервале [a, b] ;

(2.143).

N ba искомый интеграл f (i ) ;

· N i= Получение равномерно распределённой случайной величины (как и других случайных распределений) является не вполне тривиальной задачей. Но она удовлетворительно решена на существующем уровне развития вычислительной техники и информатики. Так, практически во всех современных языках программирования имеются средства для моделирования простейших случайных величин, в частности, равно мерного распределения на интервале.

Рассмотрим теперь задачу определения площади фигуры с криво линейными границами (Рис. 2.26). Погрузим её в прямоугольник со сторонами, параллельными координатным осям, имеющий известные размеры, и станем случайным образом раскидывать точки внутри это го прямоугольника. Ясно, что при равномерном распределении случай ных бросаний вероятность попадания точки в рассматриваемую фигу ру равна отношению площадей этой фигуры и объемлющего её пря 184 2. Численные методы анализа моугольника. С другой стороны, это отношение будет приблизительно равно относительной доле количества точек, которые попали в фигу ру. Оно может быть вычислено в достаточно длинной серии случайных бросаний точек в прямоугольник.

На основе сформулированной выше идеи можно реализовать ещё один способ вычисления интеграла от функции одной переменной. По мещаем криволинейную трапецию, ограниченную графиком интегри руемой функции, в прямоугольник на плоскости 0xy. Затем организуем равномерное случайное бросание точек в этом прямоугольнике и под считываем относительную частоту точек, попадающих ниже графика интегрируемой функции. Искомый интеграл равен её произведению на площадь большого прямоугольника (см. Рис. 2.27).

x Рис. 2.27. Один из способов приближённого вычисления определённого интеграла методом Монте-Карло Результаты вычислений по методу Монте-Карло сами являются слу чайной величиной, и два результата различных решений одной и той же задачи, вообще говоря, могут отличаться друг от друга. Можно показать, что второй (геометрический) способ вычисления интеграла методом Монте-Карло, вообще говоря, уступает по качеству результа тов первому способу, основанному на нахождении среднего значения функции, так как дисперсия (среднеквадратичный разброс) получае мых оценок у него больше [30].

Сформулированные выше идеи и основанные на них алгоритмы в действительности применимы для интегрирования функций от произ вольного количества переменных. Более того, вероятностные оценки погрешности, пропорциональные N 1/2, также не зависят от размер 2.17. Правило Рунге для оценки погрешности ности n пространства, в котором берётся интеграл, тогда как для тра диционных детерминистских методов интегрирования они ухудшаются с ростом n. Начиная c 7–8 переменных методы Монте-Карло уже пре восходят по своей эффективности классические кубатурные формулы и являются главным методом вычисления многомерных интегралов.


В заключение параграфа краткий исторический очерк. Идея мо делирования случайных явлений очень стара. В современной истории науки использование статистического моделировния для решения кон кретных практических задач можно отсчитывать с конца XVIII века, когда Ж.-Л. Бюффоном (в 1777 году) был предложен способ определе ния числа с помощью случайных бросаний иглы на бумагу, разграф лённую параллельными линиями.15 Тем не менее, идея использования случайности при решении различных задач не получила большого раз вития вплоть до Второй мировой войны, т. е. до середины XX века.

В 1944 году в связи с работами по созданию атомной бомбы в США, поставившими ряд очень больших и сложных задач, С. Улам и Дж. фон Нейман предложили широко использовать для их решения статисти ческое моделирование и аппарат теории вероятностей.16 Этому спо собствовало появление к тому времени электронных вычислительных машин, позволивших быстро выполнять многократные статистические испытания (Дж. фон Нейман также принимал активное участие в со здании первых цифровых ЭВМ). С конца 40-х годов XX века начинает ся широкое развитие метода Монте-Карло и методов статистического моделирования во всём мире. В настоящее время их успешно применя ют для решения самых разнообразных задач практики (см., к примеру, [30, 55] и цитированную там литературу).

2.17 Правило Рунге для оценки погрешности Предположим, что нам необходимо численно найти интеграл или производную функции, либо решение дифференциального или инте грального уравнения, т. е. решить какую-либо задачу, где фигурирует сетка на некотором интервале вещественной оси. Пусть для решения 15 Наиболее известная докомпьютерная реализация метода Бюффона была осу ществлена американским астрономом А. Холлом [68].

16 Интересно, что примерно в те же самые годы в СССР решение аналогичных задач советского атомного проекта было успешно выполнено другими методами.

186 2. Численные методы анализа этой задачи применяется численный метод порядка p, так что главный член его погрешности равен Chp, где h шаг рассматриваемой сетки, аC величина, напрямую от h не зависящая. Как правило, значение C не известно точно и его нахождение непосредственно из исходных данных задачи является делом трудным и малоперспективным. Мы могли видеть, к примеру, что для задач интерполирования и числен ного интегрирования выражение для этой константы вовлекает оценки для производных высоких порядков от подинтегральной функции. Во многих случаях их практическое вычисление не представляется воз можным, так что оценки эти носят, главным образом, теоретический характер. Аналогична ситуация и с другими постановками задач и их погрешностями решения.

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

Предположим для простоты анализа, что численные решения рас сматриваемой задачи расчитаны на сетках с шагом h и h/2 и равны соответственно Ih и Ih/2, а точное решение есть I. Тогда Ih I = Chp, p hp h Ih/2 I = C =C.

2p Вычитая второе равенство из первого, получим hp 2p Ih Ih/2 = Chp C = Chp, p 2p так что 2p Ih Ih/ · C=.

2p 1 hp Зная константу C, можно уже находить оценку погрешности расчитан ных решений Ih, Ih/2 или любых других Правило Рунге работает плохо, если главный член погрешности Chp не доминирует над последующими членами её разложения, которые соответствуют (p + 1)-ой и более высоким степеням шага сетки h. Это происходит, как правило, для сильно меняющихся решений.

Литература к главе Литература к главе Основная [1] Барахнин В.Б., Шапеев В.П. Введение в численный анализ. – Санкт Петербург–Москва– Краснодар: Лань, 2005.

[2] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. – Москва: Бином, 2003, а также другие издания этой книги.

[3] Бахвалов Н.С., Корнев А.А., Чижонков Е.В. Численные методы. Реше ния задач и упражнения. – Москва: Дрофа, 2008.

[4] Березин И.С., Жидков Н.П. Методы вычислений. Т. 1–2. – Москва: Наука, 1966.

[5] Вержбицкий В.М. Численные методы. Части 1–2. – Москва: Оникс 21 век, 2005.

[6] Волков Е.А. Численные методы. – Москва: Наука, 1987.

[7] Гельфонд А.О. Исчисление конечных разностей. – Москва: Наука, 1967.

[8] Гончаров В.Л. Теория интерполирования и приближения функций. – Москва: ГИТТЛ, 1954.

[9] Даугавет И.К. Введение в теорию приближения функций. – Ленинград: Из дательство Ленинградского университета, 1977.

[10] Демидович Б.П., Марон А.А. Основы вычислительной математики. – Москва: Наука, 1970.

[11] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн функций. – Москва: Наука, 1980.

[12] Калиткин Н.Н. Численные методы. – Москва: Наука, 1978.

[13] Кобков В.В., Шокин Ю.И. Сплайн-функции в численном анализе. – Ново сибирск: Издательство НГУ, 1983.

[14] Коллатц Л. Функциональный анализ и вычислительная математика. – Москва: Мир, 1969.

[15] Кострикин А.Н. Введение в алгебру. Часть 1. Основы алгебры. – Москва:

Физматлит, 2001.

[16] Крылов А.Н. Лекции о приближённых вычислениях. – Москва: ГИТТЛ, 1954, а также более ранние издания.

[17] Крылов В.И. Приближённое вычисление интегралов. – Москва: Наука, 1967.

[18] Крылов В.И., Бобков В.В., Монастырный П.И. Вычислительные методы.

Т. 1–2. – Москва: Наука, 1976.

[19] Кунц К.С. Численный анализ. – Киев: Техника, 1964.

[20] Люстерник Л.А., Червоненкис О.А., Янпольский А.Р. Математический анализ. Вычисление элементарных функций. – Москва: ГИФМЛ, 1963.

[21] Мак-Кракен Д., Дорн У. Численные методы и программирование на ФОР Москва: Мир, 1977.

ТРАНе.

188 2. Численные методы анализа [22] Марков А.А. Исчисление конечных разностей. Одесса: Mathesis, 1910.

[23] Мацокин А.М., Сорокин С.Б. Численные методы. Часть 1. Численный ана лиз. – Новосибирск: НГУ, 2006.

[24] Миньков С.Л., Миньков Л.Л. Основы численных методов. – Томск: Изда тельство научно-технической литературы, 2005.

[25] Мысовских И.П. Интерполяционные кубатурные формулы. – Москва: Нау ка, 1981.

[26] Натансон И.П. Конструктивная теория функций. – Москва–Ленинград:

ГИТТЛ, 1949.

[27] Никольский С.М. Квадратурные формулы. – Москва: Наука, 1988.

[28] Ортега Дж., Пул У. Введение в численные методы решения дифференци альных уравнений. – Москва: Наука, 1986.

[29] Самарский А.А., Гулин А.В. Численные методы. – Москва:

Наука, 1989.

[30] Соболь И.М. Численные методы Монте-Карло. – Москва: Наука, 1973.

[31] Стечкин С.Б., Субботин Ю.Н. Сплайны в вычислительной математике. – Москва: Наука, 1976.

[32] Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. – Москва: Наука, 1974.

[33] Тыртышников Е.Е. Матричный анализ и линейная алгебра. – Москва: Физ матлит, 2007.

[34] Тыртышников Е.Е. Методы численного анализа. – Москва: Академия, 2007.

[35] Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления.

Т. 1, 2. – Москва: Наука, 1966.

Дополнительная [36] Абрамовиц М., Стиган И. Таблицы специальных функций. – Москва: Наука, 1979.

[37] Алберг Дж., Нильсон Э., Уолш Дж. Теория сплайнов и её приложения. – Москва: Мир, 1972.

[38] Бабенко К.И. Основы численного анализа. – Москва: Наука, 1986.

[39] Геронимус Я.Л. Теория ортогональных многочленов. – Москва: Госуд. изд-во технико-теоретической литературы, 1950.

[40] Гурвиц А., Курант Р. Теория функций. – Москва: Наука, Физматлит, 1968.

[41] Демиденко Е.З. Оптимизация и регрессия. – Москва: Наука, 1989.

[42] Дробышевич В.И., Дымников В.П., Ривин Г.С. Задачи по вычислитель ной математике. – Москва: Наука, 1980.

[43] Карлин С., Стадден В. Чебышевские системы и их применение в анализе и статистике. – Москва: Наука, 1976.

Литература к главе [44] Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспе чение. – Москва: Мир, 1998.

[45] Квасов Б.И. Методы изогеометрической аппроксимации сплайнами. – Москва: Физматлит, 2006.

[46] Коллатц Л., Крабс В. Теория приближений. Чебышёвские приближения и их приложения. – Москва: Наука, 1978.

[47] Кронрод А.С. Узлы и веса квадратурных формул. Шестнадцатизначные таб лицы. – Москва: Наука, 1964.

[48] Крылов В.И., Шульгина Л.Т. Справочная книга по численному интегриро ванию. – Москва: Наука, 1966.

[49] Кузьмин Р.О. К теории механических квадратур // Известия Ленинградско го политехнического института. Отделение техн. естеств. и матем. 1931. – Т. 33.

– С. 5–14.

[50] Локуциевский О.В., Гавриков М.Б. Начала численного анализа. – Москва:

ТОО Янус, 1994.

[51] Лоран Ж.-П. Аппроксимация и оптимизация. – Москва: Мир, 1975.

[52] Меньшиков Г.Г. Локализующие вычисления. Конспект лекций. – Санкт Петербург: СПбГУ, Факультет прикладной математики–процессов управле ния, 2003.

[53] Микеладзе Ш.Е. Численные методы математического анализа. – Москва:

ГИТТЛ, 1953.

[54] Милн В.Э. Численный анализ. – Москва: Издательство иностранной литера туры, 1951.

[55] Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование.

Методы Монте-Карло. – Москва: Изд. центр Академия, 2006.

[56] Мысовских И.П. Лекции по методам вычислений. – Санкт-Петербург: Изда тельство Санкт-Петербургского университета, 1998.

[57] Погорелов А.И. Дифференциальная геометрия. – Москва: Наука, 1974.


[58] Ремез Е.Я. Основы численных методов чебышёвского приближения. – Киев:

Наукова думка, 1969.

[59] Сегё Г. Ортогональные многочлены. – Москва: Физматлит, 1962.

[60] Соболев С.Л. Введение в теорию кубатурных формул. – Москва: Наука, 1974.

[61] Стеклов В.А. О приближённом вычислении определённых интегралов // Из вестия Академии Наук. – 1916. – Т. 10, №6. – С. 169–186.

[62] Суетин П.К. Классические ортогональные многочлены. – Москва: Наука, 1979.

[63] Стефенсен И.Ф. Теория интерполяции. – Москва: Объединённое научно техническое издательство НКТП СССР, 1935.

[64] Хансен Э., Уолстер Дж.У. Глобальная оптимизация с помощью методов интервального анализа. – Москва-Ижевск: Издательство РХД, 2012.

190 2. Численные методы анализа [65] Хаусхолдер А.С. Основы численного анализа. – Москва: Издательство ино странной литературы, 1956.

[66] Хемминг Р.В. Численные методы. – Москва: Наука, 1972.

[67] Aberth O. Precise numerical methods using C++. – San Diego: Academic Press, 1998.

[68] Hall A. On an experimental determination of // Messenger of Mathematics. – 1873. – Vol. 2. – P. 113-114.

[69] Moore R.E., Kearfott R.B., Cloud M. Introduction to interval analysis. – Philadelphia: SIAM, 2009.

[70] Polya G. Uber Konvergenz von Quadraturverfahren // Mathematische Zeitschrift.

– 1933. – Bd. 37. – S. 264–286.

[71] Schoenberg I.J Contributions to the problem of approximation of equidistant data by analytic functions. Part A: On the problem of smoothing or graduation. A rst class of analytic approximation formulae. Part B: On the problem of osculatory interpolation. A second class of analytic approximation formulae // Quart. Appl.

Math. – 1946. – Vol. 4. – P. 45–99, 112–141.

Глава Численные методы линейной алгебры 3.1 Задачи вычислительной линейной алгебры Численные методы линейной алгебры это один из классических разделов вычислительной математики, который в середине XX века вычленился даже в отдельное научное направление в связи с бурным развитием математических вычислений на ЭВМ.1 Традиционный, ис торически сложившийся список задач вычислительной линейной ал гебры по состоянию на 50–60-е годы прошлого века можно найти в капитальной книге Д.К. Фаддеева и В.Н. Фаддеевой [44]. Он включал • решение систем линейных алгебраических уравнений, • вычисление определителей матриц, • нахождение обратной матрицы, • нахождение собственных значений и собственных векторов матриц, а также многочисленные разновидности этих задач.

1 В англоязычной учебной и научной литературе часто используют термин мат ричные вычисления, который уже по объёму, не охватывая, к примеру, такую часть современной вычислительной линейной алгебры как тензорные вычисления.

192 3. Численные методы линейной алгебры Но всё течёт, всё меняется. По мере развития науки и технологий в фокусе развития вычислительной линейной алгебры оказались новые задачи. Вот как формулирует список важнейших задач в 2001 году американский специалист Дж. Деммель в книге [13]:

• решение систем линейных алгебраических уравнений;

• линейная задача о наименьших квадратах:

найти вектор x, минимизирующий Ax b, Ax b для заданных mn-матрицы A и m-вектора b;

• нахождение собственных значений и собственных векторов матриц;

• нахождение сингулярных чисел и сингулярных векторов матриц.

Постановку последней задачи мы будем обсуждать ниже в §3.2г.

Вторая задача из этого списка линейная задача о наименьших квад ратах является одним из вариантов дискретной задачи о наилучшем среднеквадратичном приближении. Она возникает обычно в связи с ре шением переопределённых систем линейных алгебраических уравнений (СЛАУ), которые, к примеру, получаются при обработке эксперимен тальных данных.

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

С точки зрения классических разделов математики решение выпи санных задач даётся вполне конструктивными способами и как будто не встречает затруднений:

– решение квадратной СЛАУ получается покомпонентно по фор муле Крамера, как частное двух определителей, которые, в свою очередь, могут быть вычислены по явной формуле;

– для вычисления собственных значений матрицы A нужно выпи сать её характеристическое (вековое) уравнение det(A I) = и найти его корни.

И так далее. Но практическая реализация этих теоретических рецептов наталкивается на почти непреодолимые трудности.

3.1. Задачи вычислительной линейной алгебры К примеру, явная формула для определителя n n-матрицы вы ражает его как сумму n! слагаемых, каждое из которых есть произ ведение n элементов из разных строк и столбцов матрицы. Раскрытие определителя по этой формуле требует n! (n 1) умножений и (n! 1) сложений, т. е. всего примерно n! n арифметических операций, и потому из-за взрывного роста факториала2 решение СЛАУ по правилу Краме ра при n 20–30 делается невозможным даже на самых современных ЭВМ.

Производительность современных ЭВМ принято выражать в так называемых флопах (сокращение от английской фразы oating point operation), и 1 флоп это одна усреднённая арифметическая опера ция в арифметике с плавающей точкой в секунду (см. §1.2). Для наи более мощных на сегодняшний день ЭВМ скорость работы измеряется так называемым петафлопами, 1015 операций с плавающей точкой в секунду. Для круглого счёта можно даже взять производительность нашего гипотетического компьютера равной 1 экзафлоп = 1018 опера ций с плавающей точкой в секунду. Решение на такой вычислительной машине системы линейных алгебраических уравнений размера по правилу Крамера, с раскрытием определителей по явной комбина торной формуле, потребует времени 30 · 30! операций 30 компонент решения · 1018 флоп · 3600 час · 24 сутки · 365 сутки сек час год лет, т. е. примерно 7.57 · 109 лет. Для сравнения, возраст Земли в насто ящее время оценивается в 4.5 · 109 лет.

Обращаясь к задаче вычисления собственных значений матрицы, напомним известную из алгебры теорема Абеля-Руффини3 : общее ал гебраическое уравнение степени выше четвёртой неразрешимо в ради калах, т. е. не существует конечная формула, выражающая решения такого уравнения через коэффициенты с помощью арифметических операций и взятия корней произвольной степени. Таким образом, для матриц размера 5 5 и более мы по необходимости должны развивать для нахождения собственных значений какие-то численные методы. К этому добавляются трудности в раскрытии определителя, который вхо дит в характеристическое уравнение матрицы.

2 Напомним в этой связи известную в математическом анализе асимптотическую n! 2n (n/e)n, где e = 2.7182818...

формулу Стирлинга 3 Иногда её называют просто теоремой Абеля (см., к примеру, [60]).

194 3. Численные методы линейной алгебры Помимо неприемлемой трудоёмкости ещё одной причиной непри годности для реальных вычислений некоторых широко известных ал горитмов из чистой математики является сильное влияние на их ре зультаты неизбежных погрешностей счёта и ввода данных. Например, очень неустойчиво к ошибкам решение СЛАУ по правилу Крамера.

3.2 Теоретическое введение 3.2а Основные понятия Вектор это упорядоченный кортеж из чисел либо объектов какой то другой природы, расположенный вертикально (вектор-столбец) или горизонтально (вектор-строка). Таким образом, если a1, a2,..., an некоторые числа, то a a это вектор-столбец, a=.

..

an а это вектор-строка.

a = a1, a2, · · ·, an По умолчанию, если не оговорено противное, условимся считать, что векторами являются вектор-столбцы. Множество векторов, компо ненты которых принадлежат R или C, мы будем обозначать через Rn или Cn. При этом нулевые векторы, т. е. векторы, все компоненты ко торых суть нули, мы традиционно обозначаем через 0.

Ненулевые векторы a и b называются коллинеарными, если a = b для некоторого скаляра. Иногда различают сонаправленные колли неарные векторы, отвечающие случаю 0, и противоположно на правленные, для которых 0. Нулевой вектор по определению кол линеарен любому вектору.

Линейной оболочкой векторов v2,..., vk называют множество все возможных линейных комбинаций этих векторов, т. е. наименьшее ли нейное подпространство, содержащее эти векторы v1,..., vk. Мы будем обозначать линейную оболочку посредством lin {v1,..., vn }, так что n i vi i R.

lin {v1,..., vn } = i= 3.2. Теоретическое введение На линейных пространствах Rn и Cn можно задать скалярные про изведения. Напомним, что в вещественном случае это положительно определённая симметричная и билинейная форма, а в комплексном положительно определённая эрмитова форма. Обычно они задаются в следующем стандартном виде n для a, b Rn (3.1) a, b = ai b i i= или n для a, b Cn. (3.2) a, b = ai bi i= Наличие скалярного произведения позволяет измерять углы между векторами и ввести очень важное понятие ортогональности векторов.

При этом пространства Rn и Cn становятся евклидовыми простран ствами, и для них справедливы многие красивые и важные свойства, существенно упрощающие математические рассуждения.

Матрица это прямоугольная таблица, составленная из чисел или каких-нибудь других объектов. Если она имеет m строк и n столбцов, то обычно её записывают в виде a11 a12... a1n a21 a22... a2n A :=..,...

...

.

...

am1 am2... amn называя aij элементами матрицы A = ( aij ). При этом мы можем отождествлять векторы с матрицами размера n 1 (вектор-столбцы) либо 1 n (вектор-строки).

Матрицы используются в современной математике для самых раз нообразных целей. В частности, если дана система, составленная из конечного числа объектов (подсистем), то взаимодействие i-го объек та с j-ым можно описывать матрицей, элементы которой суть aij. В простейшем случае эти элементы принимают значения 1 или 0, соот ветствующие ситуациям связь есть и никак не связано. Для нашего курса особенно важно, что с помощью матриц, как это показывается в линейной алгебре, даётся удобное представление для линейных отоб ражений конечномерных пространств.

196 3. Численные методы линейной алгебры Ведущей подматрицей некоторой матрицы называется матрица, со ставленная из строк и столбцов с первыми номерами. Ведущий минор матрицы это определитель ведущей подматрицы.

Транспонированной к m n-матрице A = ( aij ) называется n m матрица A, в которой ij-ым элементом является aji. Иными словами, a11 a21... an a12 a22... an A :=.

..

...

...

.

...

a1m a2m... anm Для числовых матриц определены сумма, разность и произведение.

Напомним, что сумма (разность) двух матриц одинакового размера есть матрица того же размера, образованная поэлементными суммами (разностями) операндов. Если A = ( aij ) m l-матрица и B = ( bij ) l n-матрица, то произведение матриц A и B есть такая m n-матрица C = ( cij ), что l cij := aik bkj.

k= Строчным рангом матрицы (или рангом по строкам) называется, как известно, количество её линейно независимых строк. Столбцовым рангом матрицы (или рангом по столбцам) называется максимальное количество её линейно независимых столбцов. В курсах линейной ал гебры показывается, что строчный и столбцовый ранги матрицы совпа дают друг с другом и равны максимальному размеру ненулевого мино ра, порождённого этой матрицей. Как следствие, мы можем говорить просто о ранге матрицы.

Квадратная матрица, все строки которой (или столбцы) линейно независимы, называется неособенной (регулярной, невырожденной). Её ранг равен, таким образом, её порядку. В противном случае матрица называется особенной (вырожденной).

В случае, когда нулевые и ненулевые элементы в матрице A струк турированы определённым образом, по отношению к A будут употреб 3.2. Теоретическое введение Рис. 3.1. Наглядные образы нижней треугольной и верхней треугольной матриц.

ляться дополнительные определяющие термины. Например, a11 a12... a1n a a22... a2n a21 a и...

....

...

..

...

ann an1 an2... ann это верхняя треугольная и нижняя треугольная матрицы соответ ственно. Равносильные термины правая треугольная и левая тре угольная матрицы. Выбор того или иного варианта названия обычно диктуется контекстом или сложившейся традицией.

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

Блочными называются матрицы вида A11 A12... A1n A21 A22... A2n.,....

...

.

...

Am1 Am2... Amn элементы Aij которых, в свою очередь, являются матрицами. Матрицы 198 3. Численные методы линейной алгебры вида A11 A11 A12... A1n A22 A22... A2n и.

.....

...

Ann Ann где внедиагональные блоки или же блоки ниже главной диагонали яв ляются нулевыми, назовём соответственно блочно-диагональными или верхними блочно треугольными (правыми блочно треугольными), см.

Рис. 3.2. Аналогичным образом определяются нижние блочно треуголь ные (левые блочно треугольные) матрицы.

0 Рис. 3.2. Наглядные образы блочно-диагональной и верхней блочно-треугольной матриц.

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

Линейная алгебра и её численные методы в некоторых ситуациях по существу требуют выхода в поле комплексных чисел C, алгебраически пополняющее вещественную ось R. Это необходимо, в частности, в свя зи с понятиями собственных чисел и собственных векторов матриц, но может также диктоваться самой содержательной постановкой задачи.

Например, привлечение комплексных чисел бывает необходимым при исследовании колебательных режимов в различных системах, так как 3.2. Теоретическое введение в силу известной из математического анализа формулы Эйлера гар монические колебания с угловой частотой обычно представляются в виде комплексной экспоненты exp(it).

Эрмитово-сопряжённой к mn-матрице A = ( aij ) называют nm матрица A, в которой ij-ым элементом является комплексно-сопря жённый aji. Иными словами, a11 a21... an a12 a22... an A :=.

.,...

...

.

...

a1m a2m... anm и эрмитово сопряжение матрицы есть композиция транспонирования и взятия комплексного сопряжения элементов.

или Рис. 3.3. Наглядные образы симметричной матрицы.

В линейной алгебре её приложениях широко используются специ альные типы матриц эрмитовы, симметричные, косоэрмитовы, ко сосимметричные, унитарные, ортогональные и т. п. Напомним, что эр митовыми матрицами называются такие комплексные матрицы A, что A = A. Симметричными матрицами4 называют матрицы, совпада ющие со своими транспонированными. Матрица Q называется уни тарной, если Q Q = I. Матрица Q называется ортогональной, если Q Q = I.

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

200 3. Численные методы линейной алгебры другими элементами системы. Это происходит, например, если связи между элементами системы носят локальный характер. В противопо ложность этому, плотно заполненными называют матрицы, которые не являются разреженными. Иными словами, в плотно заполненных матрицах большинство элементов не равны нулю.

0 0 Рис. 3.4. Наглядные образы некоторых ленточных матриц.

Ленточными матрицами называют матрицы, у которых ненулевые элементы образуют выраженную ленту вокруг главной диагонали. В формальных терминах, матрица A = (aij ) называется ленточной, если существуют такие натуральные числа p и q, что aij = 0 при j i p и i j q. В этом случае велична p + q + 1 называется шириной ленты. Простейшими и важнейшими из ленточных матриц являются трёхдиагональные матрицы, для которых p = q = 1, и двухдиагональ ные матрицы, для которых p = 0 и q = 1 или p = 1 и q = 0. Такие матрицы встретятся нам в §3.8.

3.2б Собственные числа и собственные векторы матрицы Как должно быть известно читателю, для квадратных веществен ных или комплексных матриц большую роль в теории и приложениях играют собственные значения и собственные векторы. Если обозна чить посредством собственное значение n n-матрицы A, а x, x = 0, её собственный вектор, то они удовлетворяют матричному уравне нию (3.3) Ax = x.

Содержательный смысл этого равенства состоит в том, что на одномер ном линейном подпространстве в Rn или Cn, порождённом собствен 3.2. Теоретическое введение ным вектором x, задаваемое матрицей A линейное преобразование дей ствует как умножение на скаляр, т. е. как растяжение или сжатие.

Собственные значения являются корнями так называемого характери стического уравнения матрицы, которое имеет вид det(AI) = 0. Для n n-матрицы это алгебраическое уравнение n-ой степени, так что для его разрешимости по существу требуется привлечение поля комплекс ных чисел C. Совокупность собственных чисел матрицы называется её спектром, и в общем случае спектр подмножество комплексной плоскости.

Наконец, широко известный факт: собственные значения эрмито вых и симметричных матриц вещественны.

Предложение 3.2.1 Пусть A mn-матрица, B nm-матрица, так что одновременно определены произведения AB и BA. Спектры матриц AB и BA могут различаться только нулём.

Доказательство. Пусть какое-нибудь ненулевое собственное зна чение матрицы AB, так что (3.4) ABu = u с некоторым вектором u = 0. Умножая это равенство слева на матрицу B, получим B(ABu) = B(u), или BA(Bu) = (Bu), причём Bu = 0, так как иначе в исходном соотношении (3.4) необходи мо должно быть = 0. Сказанное означает, что вектор Bu является собственным вектором матрицы BA, отвечающим такому же собствен ному значению.

И наоборот, если ненулевое µ есть собственное значение для BA, то, домножая слева равенство BAv = µv на матрицу A, получим ABAv = AB(Av) = µ(Av), причём Av = 0. Как следствие, Av есть собственный вектор матрицы AB, отвечающий собственному значению µ. Иными словами, ненулевые 202 3. Численные методы линейной алгебры собственные числа матриц AB и BA находятся во взаимнооднозначном соответствии друг с другом.

Другой вывод этого результата можно найти, к примеру, в [2, 42].

Собственные векторы x, являющиеся решеними (3.3), называют так же правыми собственными векторами, поскольку они соответствуют умножению на матрицу справа. Но нередко возникает необходимость рассмотрения левых собственных векторов, обладающих свойством, аналогичным (3.3), но при умножении на матрицу слева. Очевидно, это должны быть собственные вектор-строки, но, имея пространство вектор-столбцов Cn, нам будет удобно записать условие на левые соб ственные векторы в виде y A = µ y, для y Cn и некоторого µ C. Применяя к этому соотношению эрми тово сопряжение, получим A y = µy, т. е. левые собственные векторы матрицы A являются правыми соб ственными векторами эрмитово сопряжённой матрицы A (что объяс няет редкость самостоятельного использования понятий левого и пра вого собственных векторов). Ясно, что при этом det(A µI) = 0.



Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 10 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.