close

Вход

Забыли?

вход по аккаунту

?

Butenina Strepetov Vychisl mat

код для вставкиСкачать
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Государственное образовательное учреждение
высшего профессионального образования
САНКТПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
Д. В. Бутенина, А. В. Стрепетов
ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА
Учебное пособие
СанктПетербург
2007
УДК 519.6
ББК 22.19
Б93
Рецензенты:
Научноисследовательский институт
прикладной математики и кибернетики
Нижегородского государственного университета;
доктор физикоматематических наук,
профессор СЗТУ И. А. Бригаднов
Утверждено редакционноиздательским советом университета
в качестве учебного пособия
Б93
Бутенина Д. В., Стрепетов А. В.
Вычислительная математика: учебное пособие / Д. В. Буте
нина, А. В. Стрепетов; ГУАП. – СПб., 2007. – 87 с.
ISBN 9785808802612
Излагаются основные вопросы курса «Вычислительная математи
ка», связанные с изучением численных методов решения линейных
и нелинейных алгебраических уравнений и систем, а также методов
приближения функций.
Пособие предназначено для студентов, обучающихся по специ
альностям 230101 «Вычислительные машины, комплексы, системы
и сети»; 230105 «Программное обеспечение вычислительной техни
ки и автоматизированных систем» заочной и очнозаочной форм обу
чения.
УДК 519.6
ББК 22.19
ISBN 9785808802612
2
© ГУАП, 2007
© Д. В. Бутенина, А. В. Стрепетов, 2007
СОДЕРЖАНИЕ
Предисловие ....................................................................
1. Аналитические основы численных методов .......................
1.1. Погрешности численного решения ............................
1.2. Метрическое пространство .......................................
1.3. Линейное нормированное пространство .....................
1.4. Вычислительные погрешности .................................
1.5. Пространство со скалярным произведением ................
1.6. Полнота метрических пространств ............................
2. Численные методы линейной алгебры ...............................
2.1. Норма и число обусловленности матрицы ...................
2.2. Прямые методы решения систем алгебраических уравне
ний .............................................................................
2.3. Итерационные методы решения систем линейных алге
браических уравнений ...................................................
2.4. Решение систем уравнений с матрицами специального вида
(трехдиагональные матрицы) .........................................
2.5. Методы решения проблемы собственных значений и соб
ственных векторов матриц .............................................
3. Численные методы решения алгебраических уравнений выс
ших порядков и трансцендентных уравнений ........................
3.1. Устойчивые решения, число обусловленности функции
одной переменной .........................................................
3.2. Отделение корней ...................................................
3.3. Методы уточнения корней уравнения ........................
4. Численные методы приближения функций ........................
4.1. Интерполяция функций ..........................................
4.2. Аппроксимация функций в среднеквадратичной метрике
Библиографический список ................................................
4
5
5
6
7
9
11
13
17
18
24
34
37
39
47
48
51
52
56
58
71
86
3
ПРЕДИСЛОВИЕ
Данное учебное пособие написано на основе курса лекций по вы
числительной математике, читаемого авторами в течение ряда лет
студентам различных специальностей университета. В пособии в сжа
том виде приводятся необходимые сведения о численных методах
решения прикладных задач, возникающих при разработке и проек
тировании вычислительных систем. При отборе материала авторы
в первую очередь руководствовались требованиями, предусмотрен
ными программой курса «Вычислительная математика». Для пони
мания содержания всех разделов пособия достаточно знания основ
ного курса высшей математики, читаемого в течение первых трех се
местров обучения. Некоторые дополнительные понятия, используе
мые в основном тексте, приведены в первом разделе, носящем
вспомогательный характер и содержащем необходимые для более
ясного понимания численных методов сведения из математического
и функционального анализа.
Как и всякое другое пособие по численным методам, эта книга
содержит изложение основных положений теории, в данном случае
относящихся к решению алгебраических уравнений и систем, а так
же методам приближения функций. Второй раздел пособия целиком
посвящен методам решения задач линейной алгебры. Способы реше
ния нелинейных уравнений приведены в третьем разделе; в четвер
том обсуждаются наиболее распространенные методы интерполяции
и аппроксимации функций.
Следует иметь в виду, что материал пособия включает отдельные,
вызывающие наибольшую трудность при обучении, вопросы програм
мы, и, тем самым, не исчерпывает полностью содержание курса «Вы
числительная математика». Среди существующей учебной литера
туры наибольшей полнотой изложения вычислительных методов
отличается двухтомник [1], достаточно подробное описание содер
жится в стандартных учебниках [2–3]. С отдельными вопросами кур
са можно более детально ознакомиться по специализированным учеб
ным пособиям [4–8], каждое из которых целиком посвящено тому
или иному разделу вычислительной математики.
Предлагаемое пособие должно служить основой для последующе
го более глубокого изучения численных методов и их применения
в различных технических задачах.
4
1. АНАЛИТИЧЕСКИЕ ОСНОВЫ ЧИСЛЕННЫХ МЕТОДОВ
Прежде чем переходить к непосредственному изложению, попыта
емся определить место теории численных методов при решении опреде
ленной научноисследовательской задачи. Как правило, в реальной си
туации явное (аналитическое) решение задачи найти невозможно – не
потому, что мы не умеем его искать, а потому, что ответ не выражается
через привычные функции – элементарные или известные специальные.
При этом развитие одной лишь вычислительной техники не снимает
проблемы: любая ЭВМ может выполнить определенную последователь
ность простейших арифметических и логических операций и не более
того. Численные методы – это и есть методы решения задач, сводящие
ся в конечном итоге к действиям, которые выполняет ЭВМ. В зависимо
сти от конкретной ситуации исследователю приходится выбирать тот
или иной численный метод из ряда имеющихся или разрабатывать но
вый, учитывающий особую специфику задачи. Выбирая (или предла
гая новый) численный метод, нужно знать о его преимуществах и недо
статках по сравнению с другими; убедиться в его сходимости и устойчи
вости; определить число необходимых операций и размеры требующей
ся памяти ЭВМ. Как правило, ответы на эти вопросы, являющиеся
частью вычислительной математики, находятся аналитически. Тем
самым вычислительная математика по своему построению и использу
емому аппарату ближе к традиционному математическому анализу
и алгебре, чем, например, к программированию. Кроме этого, любой
численный метод не зависит от того, на каком из языков программи
рования и на какой ЭВМ будет реализован алгоритм (за исключением,
быть может, развивающихся в последнее время методов параллель
ных вычислений, предусматривающих одновременную работу несколь
ких процессоров). Поэтому в данном пособии основное внимание уделя
ется описанию метода, доказательству его корректности, анализу рабо
ты метода и априорному оцениванию возникающих погрешностей.
1.1. Погрешности численного решения
Решение, полученное численным способом, всегда является при
ближенным. Источников погрешности приближенного решения не
сколько: это, вопервых, несоответствие математической модели изу
чаемому реальному явлению и погрешности исходных данных; во
вторых, погрешности численного метода решения и, втретьих, по
грешности округлений при выполнении арифметических и других
действий над числами.
5
Погрешность, обусловленную неадекватностью математической
модели и погрешностью входных параметров, называют неустрани
мой погрешностью. Оценку влияния такой погрешности можно по
лучить, решая задачу с различными начальными данными в преде
лах погрешности и сравнивая полученные результаты.
Численные методы сами по себе являются в большинстве случаев
приближенными. Даже при точных начальных данных и сколь угод
но точных вычислениях решение исходной задачи получается с не
которой погрешностью, называемой погрешностью численного ме
тода. Такая погрешность появляется, когда метод основывается на
некотором бесконечном сходящемся процессе, который при вычис
лениях обрывается на некотором шаге, либо когда численным мето
дом решается не исходная задача, а более простая, ее аппроксимирую
щая. Анализ или оценку погрешности метода иногда удается провес
ти аналитическими средствами, выразив ее через известные величины
(например, число итераций или шаг при использовании сеточных ме
тодов). Это позволяет при вычислениях эффективно управлять погреш
ностью метода, делая ее меньшей, чем неустранимая погрешность.
Намного более сложным является учет погрешности округления
в арифметических действиях – вычислительной погрешности. Ана
лиз вычислительной погрешности необходим, когда решение задачи
на ЭВМ требует большого числа операций, например при решении
уравнений с частными производными. Причем в этой ситуации ока
зывается, что учитывать влияние погрешности округления в каж
дом действии нереально, и, кроме того, даже если такой учет и прове
сти, получится слишком завышенная оценка погрешности. К сожале
нию, эффективно применяемые оценки вычислительной погрешности
на основе статистического анализа влияния округлений требуют глу
бокого понимания теории систем случайных величин, и поэтому их
описание далеко выходит за рамки настоящего пособия.
Для успешного применения численного метода необходимо, что
бы его собственная погрешность была в несколько раз меньше неуст
ранимой погрешности, а вычислительная погрешность, в свою оче
редь, в несколько раз меньше погрешности метода.
1.2. Метрическое пространство
В результате решения задачи численным методом мы хотим полу
чить результат, который был бы близок к точному ответу. Что такое
«близко»? Очевидно, для двух чисел x1 и x2 надо требовать малости
6
|x2 – x1|; близость же двух векторов или функций можно определить
разными способами. Эти вопросы рассматриваются в функциональ
ном анализе, некоторые основные понятия которого будут сейчас
изложены.
Множество M элементов произвольной природы (векторов, функ
ций и т. п.) называется метрическим пространством, если в нем опре
делена числовая функция ρ(x, y) для каждой пары элементов x и y из
M, называемая расстоянием или метрикой и удовлетворяющая сле
дующим аксиомам
1. Положительность: ρ(x, y) ≥ 0, причем ρ(x, y) = 0 тогда и только
тогда, когда x = y.
2. Симметричность:
ρ(x, y) = ρ(y, x).
(1.1)
3. Неравенство треугольника: ρ(x, y) ≤ ρ(x, z) + ρ(y, z), ∀ x, y, z ∈ M.
Эти аксиомы выполнены, разумеется, для обычного евклидова
расстояния в пространстве nмерных векторов Rn:
ρ(x, y ) =
n
∑ (xk − yk )2 .
(1.2)
k=1
Соответственно, любое подмножество пространства Rn, снабжен
ное метрикой (1.2), будет метрическим пространством. Например,
замкнутый шар Br(a) = {x: ρ(x, a) ≤ r} часто будет рассматриваться
как метрическое пространство. Другие примеры будут приведены
в следующем подразделе.
1.3. Линейное нормированное пространство
Часто множество M, в котором ищется возможное решение задачи,
является линейным пространством, когда в M определены сложение
элементов и умножение их на действительные числа, не выводящие за
пределы M. Более строго, M – линейное пространство, если выполне
ны следующие условия (здесь α, β – числа; f, g, h – элементы M):
f + g ∈ M;
α f ∈ M;
f + g = g + f;
7
(f + g ) + h = f + ( g + h);
существует нулевой элемент 0 ∈ M такой, что
f + 0 = f;
0 ⋅ f = 0;
(α + β)f = αf + βf;
α( f + g ) = αf + αg ;
α(βf ) = (αβ)f;
1 ⋅ f = f.
(1.3)
Из функциональных пространств чаще всего мы будем иметь дело
с линейным пространством непрерывных на интервале [a, b] функ
ций, которое обозначается C[a,b]. Функция f(x) принадлежит C[a,b],
если она определена на отрезке [a, b], и существует функция f*(x),
заданная на (a – ε1, b + ε2) при ε1,2 > 0, непрерывная на этом интервале
и совпадающая с f(x) на [a, b]. Выполнение аксиом (1.3) для про
странства C[a,b] достаточно очевидно; роль нулевого элемента, разу
меется, играет функция, тождественно равная нулю на интервале
[a, b]. Аналогичным образом определяется Ck[a,b] – пространство k раз
непрерывно дифференцируемых функций на интервале [a, b].
Линейное пространство называется нормированным простран
ством, если каждому элементу f поставлено в соответствие действи
тельное число f , которое называется нормой f и должно удовлетво
рять аксиомам нормы:
1. f ≥ 0, причем f = 0 тогда и только тогда, когда f = 0, т. е. f
является нулевым элементом линейного пространства.
2. αf = α f .
3. f + g ≤ f + g – неравенство треугольника для нормы.
В одном и том же линейном пространстве норму можно вводить
разными способами. Рассмотрим это на примере пространства C[a,b].
В качестве нормы в C[a,b] можно взять максимальное отклонение
функции от нуля:
(1.4)
f C = max f (x) .
[a,b]
Выполнение первых двух аксиом для нормы (1.4) очевидно, а не
равенство треугольника следует из теоремы Вейерштрасса для непре
рывных функций о существовании точки x*, в которой достигается
максимум функции, в частности, max|f(x) + g(x)|. Тогда
f+g
C
= max f (x) + g (x) = f (x* ) + g (x* ) ≤ f (x* ) + g (x* ) ≤
[a,b ]
≤ max f (x) + max g (x) = f
[a,b]
8
[a,b]
C
+ g C.
Другой способ задания нормы в C[a,b] получается, если взять сред
неквадратичное отклонение функции на интервале [a, b]:
1
⎛b
⎞ 2
2
(1.5)
f L = ⎜ ∫ f (x) dx ⎟ .
⎜
⎟
2
⎝a
⎠
Выполнение аксиом нормы в этом случае следует из соответствую
щих свойств интеграла. Выбор той или иной нормы часто определя
ется конкретными условиями задачи. Например, при рассмотрении
равномерного приближения функции необходимо пользоваться нор
мой (1.4), а при приближении методом наименьших квадратов ра
зумно взять норму (1.5). Отметим, что любое линейное нормирован
ное пространство является метрическим, поскольку в качестве рас
стояния можно взять норму разности:
ρ(f, g ) = f − g .
(1.6)
В этом случае говорят, что метрика индуцируется нормой. Аксио
мы расстояния 1 и 2 выполняются в силу аксиом 1, 2 нормы, а нера
венство треугольника для расстояния следует из неравенства треу
гольника для нормы:
ρ(f, g ) = f − g = (f − h) + (h − g ) ≤ f − h + h − g = ρ(f, h) + ρ(h, g ).
В соответствии с этим пространство C[a,b] является метрическим
пространством с равномерной ρC (f, g ) = f − g
ной ρL2 (f, g ) = f − g
L2
C
либо среднеквадратич
метрикой.
1.4. Вычислительные погрешности
Пусть x – точное, вообще говоря, неизвестное решение задачи
в некотором линейном нормированном пространстве; x* – найденное
приближенным методом численное решение. Величина
Δ (x* ) = x − x*
называется абсолютной погрешностью приближенного решения x*;
а величина
Δ (x* )
δ(x* ) =
x*
называется его относительной погрешностью. Точное значение по
грешности, разумеется, также неизвестно, как и точное решение x.
9
Можно ввести предельную погрешность, которая оценивала бы на
верняка погрешность вычислений. Любые два числа Δ(x* ) и δ(x * ),
которые оценивают, соответственно, погрешности Δ(x*) и δ(x*):
Δ(x* ) ≤ Δ(x* ), δ(x* ) ≤ δ(x* ),
называются предельной абсолютной и относительной погрешностью.
В простейших ситуациях, например, когда речь идет об арифмети
ческих операциях над числами, учет погрешности можно провести
в явной форме. Очевидно, если c = a + b, то
Δ(c* ) ≤ Δ(a* ) + Δ(b* )
и, следовательно, в качестве предельной погрешности можно взять
сумму
Δ(c* ) = Δ(a* ) + Δ(b* ).
Таким образом, при сложении (и, аналогично, вычитании) двух
чисел их абсолютные предельные погрешности складываются.
Пусть c = a ⋅ b, c* = a* ⋅ b* . Тогда
Δ(c* ) = c − c* = a ⋅ b − a* ⋅ b* = a ⋅ b − a ⋅ b* + a ⋅ b* − a* ⋅ b* ≤
≤ a ⋅ b − b* + b* ⋅ a − a* = a − a* + a* ⋅ b − b* + b* ⋅ a − a* ≤
≤ a* ⋅ b − b* + b* ⋅ a − a* + a − a* ⋅ b − b* .
Поскольку обычно погрешность значительно меньше абсолютно
го значения числа, последним слагаемым можно пренебречь, записав
Δ(c* ) ≈ a* ⋅ Δ(b* ) + b* ⋅ Δ(a* ).
В случае произведения двух приближенных чисел удобнее перей
ти к относительной погрешности
δ(c ) ≈
*
a* ⋅ Δ(b* ) + b* ⋅ Δ(a* )
a* ⋅ b*
≈
Δ(a* )
a*
+
Δ(b* )
b*
= δ(a* ) + δ(b* ).
Аналогичный результат получается и в случае деления. Итак, при
умножении и делении складываются предельные относительные по
грешности приближенных чисел. Получение оценок погрешности
10
в более сложных ситуациях не столь очевидно и требует применения
специальных методов, различных для различных задач. Они будут
изложены в соответствующих разделах пособия при анализе уже кон
кретных вычислительных процедур.
1.5. Пространство со скалярным произведением
Рассмотрим линейное пространство E. Говорят, что в линейном
пространстве введено скалярное произведение, если каждой паре f
и g элементов из E поставлено в соответствие действительное число
〈f, g〉, называемое скалярным произведением элементов f и g, которое
удовлетворяет следующим аксиомам скалярного произведения:
1. 〈f, g〉 = 〈g, f〉.
2. 〈f, f〉 ≥ 0, причем 〈f, f〉 = 0 тогда и только тогда, когда f = 0.
3. 〈αf, g〉 = α〈f, g〉.
4. 〈f1 + f2, g〉 = 〈f1, g〉 + 〈f2, g〉.
Линейное пространство со скалярным произведением называется
евклидовым пространством. В дальнейшем мы будем пользоваться
двумя конкретными евклидовыми пространствами. Первое – про
странство непрерывных на отрезке [a, b] функций со скалярным про
изведением
b
f, g = ∫ f (x) g (x)dx.
(1.7)
a
Второе – линейное пространство En+1 функций, заданных на диск
ретном множестве точек x0, x1, …, xn некоторого отрезка [a, b] со ска
лярным произведением
f, g =
1 n
∑ f (xk )g(xk ).
n + 1 k= 0
(1.8)
В последнем случае функцию f ∈ En+1 можно считать (n + 1)мер
ным вектором, соответственно пространство En+1 эквивалентно про
странству векторов. Убедимся, что (1.7) является скалярным произ
ведением в C[a,b]. Выполнение аксиомы 1 очевидно, проверка аксиом
3 и 4 проводится с помощью свойства линейности интеграла (1.7).
Проверим вторую аксиому.
b
Предположим, что 〈f, f〉 = 0, т. е. ∫ f 2 (x)dx = 0 и допустим, что f ≠ 0
a
(функция f(x) не равна тождественно нулю на [a, b]). Тогда в силу
11
непрерывности f(x) найдется отрезок [α, β] ⊂ [a, b], на котором f2(x) > 0.
Обозначим m = min f 2 (x), m > 0. Теперь можно получить оценку
x∈[α,β]
β
b
2
2
∫ f (x)dx ≥ ∫ f (x)dx ≥ m(β − α) > 0,
α
a
которая противоречит допущению. Следовательно, f(x) = 0.
Отметим, что всякое евклидово пространство всегда является ли
нейным нормированным пространством с нормой
f =
f, f
(1.9)
и, следовательно, метрическим пространством с расстоянием
ρ(f, g ) = f − g =
f − g, f − g .
(1.10)
Для нормы, заданной по формуле (1.9) через скалярное произве
дение, выполнены все три аксиомы нормы. Проверка первых двух
аксиом не вызывает затруднений. Для того чтобы убедиться в вы
полнении неравенства треугольника, необходимы предварительные
рассуждения. Любое скалярное произведение удовлетворяет соотно
шению
f, g ≤
f , f ⋅ g, g ,
(1.11)
которое следует из положительности. Действительно, рассмотрим
элемент f + αg с произвольным числом α и вычислим
f + α g, f + α g = f, f + 2α f, g + α2 g, g .
(1.12)
Получившееся выражение не может быть отрицательным ни при
каких α; следовательно, если рассматривать (1.12) как квадратич
ный трехчлен по α, его дискриминант должен быть отрицательным:
f, g
2
− f, f ⋅ g, g < 0,
откуда (1.11) немедленно следует. Теперь можно непосредственным
вычислением проверить неравенство треугольника:
f+g =
≤
f
2
f + g, f + g =
+2 f ⋅ g + g
2
f, f + 2 f, g + g, g ≤
=
(f
+ g
)
2
= f + g.
Таким образом, (1.9) действительно задает норму в пространстве
со скалярным произведением.
12
1.6. Полнота метрических пространств
При решении многих задач численный ответ ищется итерацион
ными способами, когда искомый результат x (элемент некоторого
метрического пространства M) представляется в виде предела после
довательных приближений xn. Поскольку осуществить явно предель
ный переход, как правило, не удается, вычисления прерывают при
достижении требуемой точности, когда ρ(x, xn) < ε. При таком методе
решения необходимо быть уверенным, что, вопервых, последова
тельность xn действительно сходится к некоторому элементу x и, во
вторых, что этот элемент x принадлежит исходному пространству
M, в котором ищется решение. Будем называть последовательность
элементов {xn} метрического пространства M сходящейся к элементу
x, если x ∈ M и (как в обычном анализе) ρ(xn , x) → 0. Последова
n→∞
тельность {xn}, удовлетворяющую условию
ρ(xn , xm ) → 0,
n,m→∞
будем называть последовательностью, сходящейся в себе, или пос
ледовательностью Коши. Как известно из курса высшей математи
ки, всякая числовая последовательность Коши является сходящей
ся. К сожалению, последовательность Коши элементов произволь
ного метрического пространства не всегда является сходящейся. Это
связано с тем, что предельный элемент, к которому сходится после
довательность, может не принадлежать исходному пространству M.
Так, в пространстве непрерывных функций C[–1,1] со среднеквадра
1
тичной метрикой ρ2 (x, y) =
довательность
∫ x(t) − y(t)
2
dt сходящаяся в себе после
−1
1
⎧
−1 ≤ t ≤ −
⎪ 0,
n
⎪
1
⎪
xn (t) = ⎨1 + nt, − < t ≤ 0
n
⎪
⎪ 1,
0 < t ≤1
⎪
⎩
стремится к функции, имеющей разрыв в точке t = 0, т. е. не являю
щейся непрерывной.
Определение. Метрическое пространство, в котором всякая пос
ледовательность Коши имеет предел, называется полным метричес
ким пространством.
13
Как ясно из приведенного выше примера, пространство C[a,b] с мет
рикой ρ2 не является полным метрическим пространством. Чтобы
можно было рассматривать всевозможные сходящиеся последователь
ности непрерывных функций, пространство C[a,b] необходимо допол
нить всеми функциями, которые являются пределами сходящихся
в себе последовательностей. Такое полное пространство функций обо
значается L2(a, b) и состоит из всех функций (не обязательно непре
рывных), для которых конечен интеграл
b
∫ f (x)
2
dx < +∞.
(1.13)
a
В пространстве L2 функции f1(x) и f2(x) не различаются, если
b
∫ f1 (x) − f2 (x)
2
dx = 0,
(1.14)
a
тем самым элемент пространства L2 – это класс эквивалентных в смыс
ле (1.14) функций, из которого можно выбирать любого представи
теля, т. е. некоторую функцию. Отметим, что пространство непре
рывных функций C[a,b] с метрикой ρc = max x(t) − y(t) является пол
t
ным пространством. Действительно, сходимость в пространстве (C,
ρc) представляет собой равномерную сходимость:
lim max x(t) − xn (t) = 0 означает
n→∞
t
∀ε > 0 ∃Nε : ∀n > Nε , ∀ t ∈ [a,b] ⇒ x(t) − xn (t) < ε;
а равномерно сходящаяся последовательность непрерывных функ
ций имеет своим пределом непрерывную функцию.
Полнота функционального пространства особенно важна для сходи
мости вычислительной процедуры. Основным инструментом доказатель
ства сходимости во многих случаях является следующая теорема.
Принцип сжимающих отображений (теорема Банаха).
Пусть M – полное метрическое пространство; ϕ – отображение M
в себя и, кроме того, ϕ – сжимающее отображение с показателем α < 1:
ρ ( ϕ(x), ϕ(y) ) ≤ αρ ( x, y ).
(1.15)
Тогда у отображения ϕ(x) есть единственная неподвижная точ
ка x*: x* = ϕ(x*), к которой сходятся последовательные приближе
14
ния xk+1 = ϕ(xk), начиная с любого x0 ∈ M, и выполнены следующие
оценки:
(
)
(
)
ρ xk , x* ≤ αkρ x0 , x* ;
(
)
ρ xk , x* ≤
α
ρ ( xk , xk−1 ).
1−α
(1.16)
(1.17)
Доказательство теоремы Банаха проведем в несколько этапов. Вна
чале убедимся в единственности неподвижной точки. Предположим,
что есть две неподвижные точки x = ϕ(x) и y = ϕ(y). Тогда ρ(x, y) =
= ρ(ϕ(x), ϕ(y)) ≤ αρ(x, y). Учитывая условия α < 1 и ρ(x, y) ≥ 0, заклю
чаем, что ρ(x, y) = 0. Таким образом, x = y. Для доказательства суще
∞
ствования неподвижной точки покажем, что {xk }k=0 , xk = ϕ(xk–1) –
последовательность Коши:
ρ(xk , xk +1 ) = ρ ( ϕ(xk−1 ), ϕ(xk ) ) ≤
≤ αρ(xk−1 , xk ) ≤ α2ρ(xk−2 , xk−1 ) ≤ 1 ≤ α kρ(x0 , x1 ).
Таким образом, при n > m и m → ∞
ρ(xn , xm ) ≤
n
∑
j =m+1
ρ(xj , xj −1 ) ≤
αm
ρ(x0 , x1 ) → 0.
1− α
Поскольку M – полное пространство, последовательность Коши
{xn} имеет предел x*, для которого выполнено соотношение
x* = lim xn = lim ϕ(xn−1 ) = ϕ(lim xn ) = ϕ(x* ),
n→∞
n→∞
n→∞
т. е. x* = ϕ(x*) – неподвижная точка отображения ϕ.
Первое из написанных неравенств (1.16) получается из условия
сжимаемости отображения ϕ с учетом неподвижности точки x*:
(
)
ρ(xk , x* ) = ρ ϕ(xk−1 ), ϕ(x* ) ≤ αρ(xk−1, x* ) ≤
≤ α ρ(xk−2 , x ) ≤ 1 ≤ αkρ(x0 , x* ).
2
*
Для проверки неравенства (1.17) воспользуемся третьей аксио
мой нормы
ρ(xk , x* ) ≤ ρ(xk , xk+1 ) + ρ(xk+1, x* ) ≤ α ρ(xk , xk−1 ) + αρ(xk , x* ).
15
Окончательно получаем
(1 − α)ρ(xk , x* ) ≤ α ρ(xk , xk−1 ),
откуда (1.17) немедленно следует.
Замечания к теореме Банаха.
1. Для применения теоремы следует учитывать, что M – именно
метрическое (не обязательно линейное) пространство. В ряде случаев
M может быть замкнутым множеством, например шаром Br(a).
2. Для получения оценки погрешности в процессе вычислений от
дают предпочтение неравенству (1.16), правая часть которого выра
жается только через известные уже на kм шаге величины. Оно явля
ется достаточно точным, поскольку автоматически учитывает воз
можность убывания погрешности на некоторых промежуточных ите
рациях быстрее, чем в α раз.
3. Даже в тех случаях, когда отображение является сжимающим
во всем линейном пространстве, полезно выделить замкнутое мно
жество (например, шар радиусом r), которое переходит в себя после
отображения ϕ. Тогда можно оценить расстояние от нулевого при
ближения до искомого решения
ρ(x0 , x* ) ≤ 2r
и сделать эффективной оценку (1.16):
ρ(xk , x* ) ≤ αk ⋅ 2r;
это позволяет приближенно оценить число итераций, необходимых
для достижения требуемой точности решения.
16
2. ЧИСЛЕННЫЕ МЕТОДЫ ЛИНЕЙНОЙ АЛГЕБРЫ
Задачи линейной алгебры являются наиболее распространенны
ми вычислительными задачами. Вопервых, к таким задачам напря
мую сводится анализ многих вычислительных систем в линейном
приближении. Вовторых, численное решение других собственно
математических задач, таких как численное интегрирование уравне
ний в частных производных или задач приближения функций, вклю
чает в себя решение алгебраических задач.
К численным методам линейной алгебры относятся численные
методы решения систем линейных алгебраических уравнений вида
Ax = b; обращения матриц, вычисления определителей detA и на
хождения всех или отдельных собственных значений λk и собствен
ных векторов ek симметричных матриц Aek = λkek.
Методы, применяемые при решении задач линейной алгебры, мож
но разделить на две группы. К первой группе принадлежат так назы
ваемые точные или прямые методы. Алгоритмы прямых методов по
зволяют найти решение системы за конечное число арифметических
действий. При условии идеально точных вычислений и начальных
данных найденное решение будет точным, но реально, разумеется,
даже «точный» метод дает приближенное решение, поэтому правиль
нее называть его прямым. Из прямых методов мы рассмотрим метод
Гаусса (метод LUразложений) и метод прогонки.
Вторую группу методов составляют приближенные методы, в час
тности, итерационные методы решения систем, когда результат
получается в виде последовательных приближений, сходящихся
в пределе к точному решению. Решение, найденное таким способом,
заведомо не является точным, но за счет увеличения числа итераций
можно достичь достаточно хорошего приближения. Наиболее рас
пространенными приближенными методами являются метод простой
итерации и метод Зейделя.
При выборе того или иного метода следует обязательно учитывать
специфику конкретной задачи. Так, помимо общих методов для про
извольных матриц существует большое число разнообразных мето
дов для матриц специального вида. Например, для ленточных мат
риц целесообразно применять метод прогонки; специальные методы
разработаны для матриц с теплицевой структурой.
Требования, предъявляемые к алгоритмам численного решения
задач линейной алгебры, разумеется, такие же, как и для любых дру
гих вычислительных процедур. Но при работе с матрицами, особен
но больших размеров, необходимо уделять внимание числу выпол
17
няемых арифметических операций и размеру необходимой памяти
ЭВМ, а также оценке погрешности, возникающей при округлении
вычислений. Так, в ЭВМ никогда не применяют правило Крамера,
поскольку, вопервых, оно требует значительно большего числа опе
раций по сравнению со всеми другими методами, и, вовторых, вы
числение определителя по формуле Крамера содержит вычитания
очень больших по абсолютной величине, но близких чисел, что при
водит к недопустимой погрешности.
Оценка числа необходимых арифметических операций получает
ся, как правило, простым вычислением и будет приведена во всех
рассматриваемых методах. Несколько сложнее дело обстоит с опре
делением погрешности. Оценку влияния неустранимой погрешно
сти на конечный результат мы получим в следующем подразделе, где
рассматривается число обусловленности – мера чувствительности си
стемы к возмущению начальных данных. Что же касается вычисли
тельной погрешности, то строгий учет предельных абсолютных по
грешностей на каждом шаге дает заведомо завышенную грубую оцен
ку. Мы лишь упомянем существование вероятностных методов опре
деления погрешности при матричных вычислениях, поскольку их
подробное изложение требует применения аппарата математической
статистики.
2.1. Норма и число обусловленности матрицы
В линейном пространстве Rn векторов x = (x1, x2, …, xn) можно
вводить различные нормы (см. разд. 1). Наиболее употребительны из
них две:
x 1 = max xj ;
1≤ j ≤n
x
=
2
x, x =
(2.1)
n
∑ xi2 .
(2.2)
i =1
Нормы (2.1) и (2.2) являются эквивалентными:
x 1 ≤ x 2 ≤ n x 1.
В пространстве квадратных матриц размером n×n также можно
различными способами определить норму, удовлетворяющую требу
емым аксиомам. Остановимся на одном из способов, когда норма мат
18
рицы вводится так называемым согласованным образом. А именно,
определим норму ||A||:
A = sup
Ax
x
x ≠0
(2.3)
,
где x , Ax – одна из введенных норм (2.1), (2.2) в пространстве век
торов. Пусть A, B – две квадратные числовые матрицы размером n×n,
A + B – их сумма. Тогда
A + B = sup
(A + B)x
x
x ≠0
≤ sup
x ≠0
Ax
x
≤ sup
Ax + Bx
x
x ≠0
+ sup
x ≠0
Bx
x
≤
= A + B.
Таким образом, неравенство треугольника для согласованной мат
ричной нормы (2.3) выполнено. Справедливость остальных матрич
ных норм очевидна.
Рассмотрим первую из введенных норм A 1 = sup Ax 1 x 1 . Лег
x ≠0
ко вычислить, что
⎛
A 1 = sup ⎜ max
x ≠0 ⎜
⎝ i
n
⎞
k=1
⎠
∑ aikxk ⎟⎟
max xi ≤
i
⎛ n
⎞
⎛ n
⎞
≤ max ⎜ ∑ aik ⎟ ⋅ max xk max xi = max ⎜ ∑ aik ⎟.
i
k
i
i
⎝ k=1
⎠
⎝ k=1
⎠
(2.4)
С другой стороны, если в качестве вектора x взять вектор с компонен
тами xk = sign(aik) (т. е. с компонентами, равными ±1 в зависимости от
знака aik), то (2.4) обратится в равенство. Тем самым доказано, что
⎛ n
⎞
A 1 = max ⎜ ∑ aik ⎟.
i
⎝ k=1
⎠
(2.5)
Вторая из введенных норм также допускает оценку
12
A
2
⎛ n n 2⎞
≤ ⎜ ∑∑ aik
⎟
⎝ i=1 k=1 ⎠
,
(2.6)
но вектора x, на котором (2.6) обращалось бы в равенство, не суще
ствует, поэтому явного выражения для нормы ||A||2 в случае произ
19
вольной матрицы нет. Лишь в случае, когда матрица A симметрич
ная (AT = A), норма ||A||2 совпадает с максимальным по модулю соб
ственным числом:
A
= λ max (A).
2
(2.7)
По определению нормы, для любого x имеем
Ax ≤ A ⋅ x ,
(2.8)
поэтому для любых матриц A и B выполнено
ABx ≤ A ⋅ Bx ≤ A ⋅ B ⋅ x
и, следовательно:
2
k
A2 x ≤ A ⋅ x , 1 , A k x ≤ A ⋅ x ,
т. е.
k
Ak ≤ A .
(2.9)
Перейдем к анализу влияния погрешности начальных данных на
погрешность решения. Рассмотрим систему линейных алгебраичес
ких уравнений с невырожденной матрицей
A x = b,
(2.10)
где detA ≠ 0, b ≠ 0. Система (2.10) имеет единственное решение. По
пробуем оценить, как погрешность правой части Δb повлияет на по
грешность решения Δx в предположении, что матрица A известна
точно, т. е. наряду с (2.10) рассмотрим систему
A (x + Δx) = b + Δb.
(2.11)
Поскольку AΔx = Δb или Δx = A–1Δb, имеем
δ(x) Δx
=
δ(b) Δb
≤
20
A ⋅x
x
x
b
⋅
=
b
x
⋅
A −1 Δb
Δb
Δx
Δb
=
Ax
x
⋅
A −1Δb
Δb
= A ⋅ A −1 = ν ( A ).
≤
Таким образом, отношение относительной погрешности решения
к относительной погрешности правой части оценивается через число
ν(A) = ||A||⋅||A–1||,
которое называется числом обусловленности матрицы A. Число обус
ловленности есть максимально возможный коэффициент усиления
относительной погрешности от правой части к решению системы:
δ ( x ) ≤ ν ( A ) ⋅δ ( b ).
Отметим основные свойства числа обусловленности:
1) ν ( E ) = 1, где E – единичная матрица;
2) ν ( A ) ≥ 1;
3) ν ( A ) = ν A −1 ;
4) ν ( AB ) ≤ ν ( A ) ν ( B );
5) ν ( αA ) = ν ( A ), α ≠ 0.
Доказательство первого свойства очевидно:
(
)
ν ( E ) = E ⋅ E −1 = E ⋅ E = 1;
при проверке второго воспользуемся неравенством для нормы матриц
1 = E = A ⋅ A −1 ≤ A ⋅ A −1 = ν ( A ).
Третье свойство следует непосредственно из определения числа обус
ловленности. Четвертое есть вновь следствие неравенства для норм
ν ( AB ) = AB ⋅ ( AB )
−1
≤ A ⋅ B ⋅ B −1 ⋅ A −1 = ν ( A ) ⋅ ν ( B ).
И, наконец, пятое справедливо изза однородности нормы:
ν ( α A ) = αA ⋅ ( αA )
−1
= α ⋅ A ⋅ α −1 ⋅ A −1 = ν ( A ).
Точное значение числа обусловленности зависит от того, какая выб
рана норма матриц, первая или вторая, но качественный вывод о поведе
нии погрешности можно сделать при любом способе определения нормы.
Насколько важно заранее оценить меру близости матрицы, пока
зывает следующий пример плохо обусловленной системы уравнений.
Рассмотрим систему
⎛ 1 −1 1 −1 ⎞
⎛ −1 ⎞
⎜
⎟
⎜ ⎟
⎜ 0 1 1 −1 ⎟ ⋅ x = ⎜ −1 ⎟.
⎜⋅ ⋅
⎜ 2 ⎟
⋅
⋅ ⎟
⎜⎜
⎟⎟
⎜⎜ ⎟⎟
⎝0 0 1 1 ⎠
⎝1⎠
(2.12)
21
Решение такой системы очевидно и находится без труда:
x = ( 0 0 1 1) .
T
Но если допустить, что правая часть системы (2.12) имела по
грешность Δb = (0 0 … ε)T, даже достаточно малую, окажется, что
погрешность решения будет весьма существенной при больших раз
мерностях задачи (когда n >> 1):
(
Δx = 2n−2 ε 2n−3 ε 1 ε
)
T
.
Например, при n = 100 максимальное отклонение будет 298⋅|ε| >
> 1029⋅|ε|. Число обусловленности такой системы очень велико:
ν(A) ≥
Δx
x
⋅
b
Δb
=
2n−2 ε 1
⋅ = 2n−2,
1
ε
с чем и связана неустойчивость решения.
Для дальнейших рассуждений нам понадобится следующее утвер
ждение, которое не только носит вспомогательный характер, но
и представляет самостоятельный интерес.
Лемма. Пусть матрица B ограничена и имеет норму ||B|| ≤ α < 1,
тогда у матрицы (E–B) существует обратная матрица (E–B)–1 с огра
ниченной нормой
( E − B )−1
≤
1
.
1− α
Доказательство сводится к представлению обратной матрицы
в виде сходящегося ряда
( E − B )−1 = E + B + B2 + 1 + Bn + 1.
Оценка нормы тогда получается при использовании соответству
ющих неравенств для норм:
(E − B)
−1
= E + B + B2 + 1 + B n + 1 ≤
2
≤ E + B + B + 1 ≤ 1 + α + α2 + 1 =
1
.
1− α
Из этой оценки немедленно следует сходимость ряда и, следова
тельно, существование обратной матрицы (E–B)–1. Для матрицы
22
(E + B)–1 лемма также справедлива и доказывается полностью ана
логично.
Оценка влияния погрешности всех начальных данных на погреш
ность решения также может быть сделана при помощи числа обус
ловленности. А именно, пусть матрица A, так же как и правая часть
b, известна с некоторой погрешностью ΔA. В такой ситуации эффек
тивная оценка погрешности решения возможна только при условии
1 − ν ( A ) ⋅δ ( A ) > 0,
(2.13)
которое мы будем предполагать выполненным. Система уравнений,
которая реально решается вместо системы (2.10), имеет вид
( A + ΔA ) (x + Δx) = b + Δb.
(2.14)
С учетом (2.10) она превращается в уравнение относительно Δx:
( A + ΔA ) Δx = Δb − ΔA x.
(2.15)
Как следует из только что доказанной леммы, при выполнении усло
вия (2.13) матрица системы (2.15) имеет обратную (E + A–1ΔA)–1A–1,
норма которой допускает оценку
(E + A
−1
ΔA
)
−1
A
−1
≤
A −1
1 − A −1ΔA
≤
A −1
1 − ν ( A ) ⋅δ ( A )
.
Тогда можно записать
Δx
x
(
≤ E + A −1ΔA
≤
)
−1
⎛ Δb
ΔA ⋅ x ⎞
+
A −1 ⋅ ⎜
⎟⎟ ≤
⎜ x
x
⎝
⎠
A −1 ⋅ A
⎛ Δb
ΔA
⋅⎜
+
⎜
1 − ν(A) ⋅ δ(A) ⎝ x ⋅ A
A
⎞
⎟⎟.
⎠
С учетом очевидного неравенства ||b|| ≤ ||A||⋅||x|| окончательно полу
чаем
δ(x ) =
ν(A)
1 − ν ( A ) ⋅δ ( A )
⋅ ( δ ( b ) + δ ( A ) ).
(2.16)
Таким образом, мы получили формулу для оценивания влияния
неустранимой погрешности начальных данных, в которую входят
23
лишь число обусловленности матрицы и относительные погрешнос
ти матрицы A и правой части b. Формула (2.16) является также по
лезной и при определении вычислительной погрешности, поскольку
погрешность округлений при арифметических операциях с матрич
ными элементами часто можно пересчитать в погрешность началь
ных данных.
В заключение отметим, что иногда можно найти число обуслов
ленности ν(A), не вычисляя обратной матрицы и ее нормы. Так, для
симметричных матриц число обусловленности, определенное по нор
ме ||A||2, есть отношение максимального к минимальному по модулю
собственному числу матрицы A:
ν(A) =
λ max ( A )
λ min ( A )
.
Способ определения только этих двух чисел будет изложен в раз
деле, посвященном решению частичной проблемы собственных зна
чений. Численно оценить меру обусловленности можно, решая не
сколько систем с близкими правыми частями и сравнивая получен
ные решения. Такой способ, разумеется, является не строгим.
2.2. Прямые методы решения систем
алгебраических уравнений
Способы решения систем линейных уравнений в основном разде
ляются на две группы. К первой группе относятся точные методы,
представляющие собой конечные алгоритмы для вычисления корней
системы (таковы, например, правило Крамера, метод Гаусса, метод
главных элементов, разложение и др.). Ко второй – итерационные
методы, позволяющие получать корни системы с заданной точнос
тью путем сходящихся бесконечных процессов (к их числу относят
ся метод простых итераций, метод Зейделя и др.). Вследствие неиз
бежных округлений результаты даже точных методов являются при
ближенными, причем оценка погрешностей корней в общем случае
затруднительна.
Точные (прямые) методы решения систем алгебраических уравне
ний используются обычно для сравнительно небольших (n < 200)
систем с плотно заполненной матрицей и неблизким к нулю опре
делителем. Эффективность прямого метода зависит от числа ари
фметических операций, необходимых для решения поставленных
24
задач. Так, в общем случае вычисление определителя порядка n
(без использования специальных приемов) требует (n – 1)n! умноже
ний и n! – 1 сложений, т. е. общее число арифметических операций
равно
n ⋅ n ! − 1.
Общее число операций для вычисления обратной матрицы равно
n2 ⋅ n ! − 1.
Рассмотрим задачу нахождения решения системы алгебраических
уравнений общего вида, т. е. нахождение x1, x2, …, xn по заданным aij
и bi (i, j = 1, n) :
a11x1 + a12x2 + 1 + a1n xn = b1 ⎫
⎪
a21x1 + a22x2 + 1 + a2n xn = b2 ⎪
⎬.
1111111111111 ⎪
an1x1 + an2x2 + 1 + ann xn = bn ⎪⎭
(2.17)
Систему (2.17) запишем в матричном виде
Ax = b,
(2.18)
n×n
где A = {aij }i,j =1,n ∈ R
– данная матица; x – искомый вектор (x1,
x2, …, xn)T; b = (b1, b2, …, bn)T – вектор с компонентами bi.
Рассмотрим векторы в пространстве Rn:
⎛1⎞
⎛0⎞
⎛0⎞
⎛0⎞
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜0⎟
⎜1 ⎟
⎜0⎟
⎜0⎟
⎜1 ⎟
⎜1 ⎟
⎜1 ⎟
⎜1 ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
e1 = ⎜ 0 ⎟, e2 = ⎜ 0 ⎟, 2, ek = ⎜ 1 ⎟ − kя строка, en = ⎜ 0 ⎟.
⎜1 ⎟
⎜1 ⎟
⎜1 ⎟
⎜1 ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜0⎟
⎜0⎟
⎜0⎟
⎜0⎟
⎜0⎟
⎜0⎟
⎜0⎟
⎜1⎟
⎝ ⎠
⎝ ⎠
⎝ ⎠
⎝ ⎠
Система векторов {ek } ∈ R n образует базис в пространстве Rn, так
как это ортонормированная система линейно независимых векторов.
25
n
Действительно, из равенства 0 = ∑ γ i ei = ( γ1, 1, γ n ) следует тотчас γi ≡
T
i =1
≡ 0 для любого i = 1, n. Любой вектор x может быть представлен в виде
линейной комбинации базисных векторов ek:
⎛ x1 ⎞
⎜ ⎟
x
x = ∑ xkek = ⎜ 2 ⎟.
⎜
1 ⎟
k=1
⎜⎜ x ⎟⎟
⎝ n⎠
n
Прежде чем применить прямые методы к матрице общего вида,
рассмотрим решение уравнения Ax = b с матрицей специального вида:
с нижней и верхней треугольной матрицами.
Прямая подстановка. Рассмотрим уравнение вида (2.18) Lx = b,
но матрица L = {lij}; lij = 0, j > i – нижняя треугольная, т. е.
= b1 ⎫
⎪
l21x1 + l22x2
= b2 ⎪
⎬.
. . . . . . .
. .⎪
ln1x1 + ln2x2 + ... + lnn xn = bn ⎪⎭
l11x1
(2.19)
Для того чтобы решить систему Lx = b с произвольной правой
частью, необходимо и достаточно, чтобы определитель матрицы L,
detL ≠ 0, т. е. lii ≠ 0, i = 1, n.
Итак, решение системы (2.19) имеет вид
b1
b −l x ⎫
, x2 = 2 21 1 , ⎪
l11
l22
⎪
................................ ⎪
⎪
i −1
⎪
bi − ∑ lij xj
⎪
j =1
⎪
xi =
,
⎬
lii
⎪
.................................⎪
⎪
n−1
⎪
bn − ∑ lnj xj
⎪
j =1
⎪
xn =
.
⎪⎭
lnn
x1 =
26
(2.20)
Обратная подстановка. Аналогично решается система с верхней
диагональной матрицей U = {uij}; uij = 0, j > i:
Ux = b.
Решение может быть записано в виде x = U–1b или
n
xn =
b − un −1,1 xn
bn
, xn −1 = n −1
, ..., x1 =
unn
un −1, n −1
b1 − ∑ u1j xj
j =2
u11
.
(2.21)
Перейдем непосредственно к методу LUразложений матрицы A.
LUразложение матрицы. Рассмотрим матрицу
⎛ 1 0 1 0 ⎞ ⎛ α1 ⎞
⎜
⎟ ⎜ ⎟
0 1 1 0⎟ ⎜ 2 ⎟
E − αeTk = ⎜
−
( 0, 0,
⎜1 1 1 1 ⎟ ⎜ 2 ⎟
⎜⎜
⎟⎟ ⎜⎜ ⎟⎟
⎝ 0 0 1 1 ⎠ ⎝ αn ⎠
⎛ 1 0 1 −α1 1
⎜
⎜ 0 1 1 −α2 1
⎜2 2
2
=⎜
1
−
αk 1
0
0
1
⎜
⎜2 2
2
⎜⎜
⎝ 0 0 1 −αn 1
1, 1, 1, 0 ) =
0⎞
⎟
0⎟
2⎟
⎟.
0⎟
2⎟
⎟
1 ⎟⎠
↑
k й столбец
Если вместо α выбрать вектор α (k):
α (k)
⎛ 0 ⎞
⎜
⎟
⎜ 0 ⎟
⎜ 1 ⎟
=⎜
⎟,
⎜ lk +1,k ⎟
⎜ 1 ⎟
⎜
⎟
⎜ l
⎟
⎝ n,k ⎠
то матрица E − α(k) eTk станет единичной нижней треугольной матри
цей вида
27
⎛1
⎜
⎜2
⎜0
Mk = ⎜
⎜0
⎜2
⎜⎜
⎝0
1
0
2
1
1
1 −lk+1,k
2
1
−lnk
1 0⎞
⎟
2⎟
2 0⎟
⎟.
1 0⎟
2⎟
⎟
1 1 ⎟⎠
(2.22)
↑
k й столбец
Матрица Mk называется матрицей Гаусса; вектор α(k) – вектором
Гаусса; коэффициенты lk+1,k, …, lnk – множителями Гаусса. Отметим,
что определитель detMk = 1 и
Mk−1 = E + α(k) eTk .
Теорема. Пусть M1, M2, …, Mk – матрицы Гаусса (2.22), k ≤ n, тогда
k
k
i =1
i =1
∏ Mi−1 = M1−1 ⋅ M2−1 ⋅1 ⋅ Mk−1 = E + ∑ α(i) еTi =
0
⎛1
⎜
⎜ l21 1
⎜1 1
=⎜
⎜ lk1 lk2
⎜1 1
⎜⎜
⎝ ln1 ln2
1 0 1 0⎞
⎟
1 0 1 0⎟
1 1 1 1⎟
⎟.
1 1 1 0⎟
1 1 1 1⎟
⎟
1 lnk 1 1 ⎟⎠
Доказательство теоремы проверяется непосредственным перемно
жением Mi−1, i = 1, 2, 1, k.
Рассмотрим далее вектор
(
)
Mk x = E − α(k) eTk x = x − α(k) eTk x = x − xk α(k)
28
⎛ x1
⎞
⎜
⎟
⎜1
⎟
⎜ xk
⎟
=⎜
⎟.
⎜ xk +1 − xk ⋅ lk +1,k ⎟
⎜1
⎟
⎜⎜
⎟⎟
⎝ xn − xk ⋅ lnk
⎠
Нетрудно увидеть, что если xk ≠ 0, то можно подобрать lij (j > i, j ≠ k)
так, чтобы
⎛ x1 ⎞
⎜ ⎟
⎜1 ⎟
⎜x ⎟
x
Mkx = ⎜ k ⎟, lik = i , i = k + 1, 2, n.
xk
⎜0 ⎟
⎜1 ⎟
⎜⎜ ⎟⎟
⎝0 ⎠
(2.23)
Итак, мы доказали следующее утверждение.
Пусть x ∈ R n, тогда, если:
1) xk = 0, то Mk x = x;
2) xk ≠ 0, то существует матрица Mk (т. е. коэффициенты lik) такая,
что справедливо (2.23).
Пусть теперь матрица A = {aij } = ( a1, a2 , 1, aq ) ∈ R n×q , Mk – матри
ца Гаусса. Тогда произведение матриц
Mk A = ( M k a1, Mk a2 , 1, Mk aq ) =
⎛ a11
⎜
⎜2
⎜a
k1
=⎜
⎜ ak+1,1 − ak1 lk+1,1
⎜
⎜2
⎜ an,1 − ak1 lnk
⎝
1 a1q
⎞
⎟
2
⎟
⎟
1 akq
⎟.
1 ak+1,q − akq lk+1,k ⎟
⎟
2
⎟
⎟
1 anq − akq lnk
⎠
Перейдем непосредственно к самому LUразложению. Пусть мат
(0)
рица A = A(0) = a1(0) , 1, a (0)
n . Если a11 ≠ 0, то выберем матрицу M1
так, чтобы
(
)
M1a1(0)
(0) ⎞
⎛ a11
⎜
⎟
a(0)
0 ⎟
1
=⎜
, li1 = i(0)
, i > 1.
⎜ 1 ⎟
a11
⎜⎜
⎟⎟
⎝ 0 ⎠
(2.24а)
29
Умножая матрицу M1 на матрицу A(0), получим матрицу A(1) =
= M1A(0) = M1A вида
(
A(1) = a1(1) , 1, a (1)
n
)
(0)
⎛ a11
⎜
⎜ 0
= M1A = ⎜
⎜ 2
⎜ 0
⎝
(0)
⎞
1 a1(0)
a12
n
⎟
(1)
(1)
a11 1 a2n ⎟
⎟.
2
2 ⎟
(1) ⎟
an(1)2 1 ann
⎠
Коэффициенты матрицы A(1) вычисляются по формулам
aij(1) = aij(0) − li1a1(0)
j , j > 1.
(2.25а)
(1)
Далее, если a22
≠ 0, то M1 выбираем так, чтобы
M2a2(1)
(0) ⎞
⎛ a12
⎜
⎟
(1)
⎜ a22
⎟
a(1)
2
= ⎜ 0 ⎟, li2 = i(1)
, i > 2.
⎜
⎟
a22
⎜ 1 ⎟
⎜⎜
⎟⎟
⎝ 0 ⎠
(2.24б)
(1)
= 0, то a1(2) = M2a1(1) = a1(1) . Следовательно:
Поскольку a12
A
(2)
=
(
a1(2) , 1, a(2)
n
)
(0)
⎛ a11
⎜
⎜ 0
⎜
= M2M1A = ⎜ 0
⎜ 2
⎜
⎜ 0
⎝
(0)
a12
(1)
a22
0
2
0
(0)
⎞
a13
1 a1(0)
n
⎟
(1)
a23
1 a2(1)
n ⎟
⎟
(2)
.
a33
1 a3(2)
n ⎟
2
2 ⎟
⎟
(2) ⎟
1
an(2)
a
nn ⎠
3
Элементы матрицы A(2) вычисляются аналогично коэффициентам
матрицы A(1):
aij(2) = aij(1) − li2a2(1)j , j > 2.
(2.25б)
Последний шаг в этой схеме: матрица A(n–1) = Mn–1Mn…M1A – это
верхняя треугольная матрица, коэффициенты которой вычисляют
ся так:
aij(n−1) = aij(n−1) − li,n−1an(n−−1,2)j , j > n − 1,
30
(2.25в)
li,n−1 =
ai(,nn−−2)
1
an(n−−1,2)n−1
, i > n − 1.
(2.24в)
Обозначим матрицу A(n–1) через U = Mn–1Mn…M1A с коэффициен
тами uij = aij(i−1) , i = 1, 1, n − 1; j > i. Умножая слева U на матрицу L =
= M1−1 ⋅ M2−1 1Mn−1−1, получим
L ⋅ U = M1−1 ⋅ M2−1 1Mn−1−1 ⋅ Mn−1 ⋅ Mn 1M1 ⋅ A = A.
Окончательно имеем M1−1 ⋅ M2−1 1Mn−1−1 ⋅ U = L ⋅ Mn−1 ⋅ Mn 1M1 ⋅ A = A.
Таким образом:
A = LU,
где матрицы L = {lij} и U = {uij} имеют соответственно вид нижней
и верхней треугольных матриц:
⎛ u11 u12 1 u1n ⎞
0 1 0⎞
⎛1
⎜
⎟
⎜
⎟
l21 1 1 0 ⎟
⎜ 0 u11 1 u2n ⎟
⎜
; U=⎜
L=
⎟.
⎜ 2
2
2⎟
2
2
2 ⎟
⎜
⎜⎜
⎟⎟
⎜ 0
⎝ ln1 ln2 1 1 ⎠
0 1 unn ⎟⎠
⎝
Коэффициенты матрицы L = {lij }, lij = 0, j > i вычисляются по фор
мулам (2.24), lii = 1, а uij ≠ 0, j ≥ i определяются по формулам (2.25),
если uij = aij(i−1) , i = 1, 1, n − 1; j > i.
(k−1)
Итак, если на каждом шаге akk ≠ 0, k = 1, n, то можно получить
такие верхнюю U и нижнюю L треугольные матрицы, что A = LU,
причем
(0)
(1)
( n−1)
det(A) = u11 ⋅ u22 ⋅1 ⋅ unn = a11
⋅ a22
⋅1 ⋅ ann
.
(2.26)
Решение системы линейных алгебраических уравнений методом
LUразложений.
Рассмотрим решение уравнения (2.18) методом LUразложений.
Первый этап: находим LUразложение матрицы A по схеме (2.24),
(2.25).
Второй этап: запишем уравнение (2.18) в виде Ax = LUx = b; обо
значим Ux = y, тогда решение (2.18) эквивалентно решению системы
⎧Ly = b
⎨
⎩Ux = y
(2.27)
с треугольными матрицами L и U. По формулам (2.21) найдем вектор
y, а затем вектор x по формулам (2.20).
31
Итак, методом LUразложений получено решение линейной сис
темы уравнений Ax = b за число арифметических операций порядка
O(n3). Одновременно вычислен определитель det матрицы A [см. фор
мулу (2.26)].
Метод Гаусса с выбором главного (ведущего) элемента.
Пусть дана система линейных уравнений Ax = b, причем у матри
цы A некоторые диагональные элементы либо малы, либо равны
нулю. Тогда метод LUразложений неприменим в том виде, который
изложен выше, и следует подправить схему вычислений так, чтобы
все оставалось законным. Рассмотрим расширенную прямоугольную
матрицу (A, b). Выберем ненулевой, как правило, наибольший по
модулю, не принадлежащий столбцу свободных членов элемент apq.
Выбор максимального по модулю элемента можно проводить несколь
кими способами: построчным, постолбцовым выбором либо поэле
ментно. Такой элемент apq матрицы A = (A, b) = (aij, ai,n+1), 1 ≤ i, j ≤ n
называется главным элементом.
Вычислим множители Гаусса
lip =
aiq
a pq
, i ≠ p.
(2.28)
Строка с номером p называется главной. Далее произведем опера
цию: от каждой неглавной строки вычтем главную, умноженную на
соответствующий множитель lip для этой строки. В результате полу
чим новую матрицу, у которой qстолбец состоит из нулей. Отбрасы
вая этот столбец и главную pстроку, получим матрицу с меньшим
числом строк и столбцов на единицу. Над матрицей производим те же
операции с выбором главного элемента и строим матрицу A(2) и т. д.:
A, A(1),…, A(n–1).
Последняя из этой последовательности матриц представляет со
бой двучленную матрицустроку; ее также считаем главной строкой.
Для определения неизвестного xi объединяем в систему все главные
строки, начиная с последней, входящей в матрицу A(n–1). После над
лежащего изменения нумерации неизвестных получается система
с треугольной матрицей, из которой легко шаг за шагом найти неиз
вестные данной системы Ax = b. Метод главных элементов всегда
применим, если det(A) ≠ 0. Сам метод Гаусса является частным слу
чаем метода главных элементов, а схема метода Гаусса получается,
если за главный элемент всегда выбирать левый верхний элемент со
32
ответствующей матрицы. Решение уравнения Ax = b сводится к ре
шению уравнения Ux = U1b, где U – верхняя треугольная матрица:
⎛ u11 u12 1 u1n ⎞
⎜
⎟
⎜ 0 u11 1 u2n ⎟
U=⎜
⎟,
2
2 ⎟
⎜ 2
⎜ 0
0 1 unn ⎟⎠
⎝
(
)
T
а U1b = b1(0) ,b2(1) , 1, bn(n–1) , где bi(i−1) = bi(i−2) − li,i−1bi(−i1−2) .
Метод Гаусса вычисления обратной матрицы. Пусть дана не
особая матрица A = {aij }; i, j = 1, n. Для нахождения ее обратной мат
рицы A–1 = {xij} используем основное соотношение
A ⋅ A −1 = A ⋅ {xij } = E.
(2.29)
Перемножая матрицы A⋅A–1, будем иметь n систем уравнений от
носительно n2 неизвестных xij:
n
∑ aikxkj = δij ,
i, j = 1, 2, 1, n,
(2.30)
k=1
где
⎧1, когда i = j,
δij = ⎨
⎩0, когда i ≠ j.
Полученные n систем линейных уравнений для j = 1, 2, …, n, име
ющих одну и ту же матрицу A и различные свободные члены, одно
временно можно решить методом Гаусса. Применяя метод Гаусса,
получим, например, для компонент jстолбца матрицы A–1
⎫
δ1(0)
j
⎪
(1)
=
δ2(1)j ⎪⎪
a22
x2 j + 1 + a2(1)
n xnj
⎬,
11111111 1
1⎪
(n−1)
(n−1) ⎪
ann
xnj = δnj
⎪⎭
(0)
(0)
a11
x1j + a12
x2 j + 1 + a1(0)
n xnj
=
(2.31)
где δij(i−1) вычислены по формулам (2.25). Количество операций по
рядка O(n3).
33
2.3. Итерационные методы решения
систем линейных алгебраических уравнений
Решения, получаемые с помощью прямых методов, обычно содер
жат погрешности, вызванные округлениями при выполнении опера
ций над числами. Рассмотрим в связи с этим методы, позволяющие
уточнить решение, полученное с помощью прямого метода, либо са
мостоятельно получить это решение.
Метод простой итерации. Найдем решение системы линейных
уравнений (2.17), предполагая, что все диагональные коэффициен
ты матрицы aii ≠ 0.
Приведем систему Ax = b к виду
x = Cx + f,
(2.32)
где C – некоторая матрица; f – векторстолбец.
Исходя из произвольного вектора x(0) строим итерационный про
цесс
x(k+1) = Cx(k) + f, k = 0,1,1
или в развернутом виде
x1(k+1) = c11x1(k) + c12x2( k) + 1 + c1n xn(k) + f1,
11111111111111111
xn(k+1) = cn1x1(k) + cn2x2(k) + 1 + cnn xn(k) + fn .
(2.33)
Произведя итерации, получим последовательность векторов x(1),
x(2), …, x(n). Можно доказать, что если элементы матрицы C удовлет
воряют одному из условий
n
∑ cij
≤ α < 1, i = 1, 2, 1, n
(2.34)
≤ β < 1, j = 1, 2, 1, n,
(2.35)
j =1
или
n
∑ cij
i =1
то процесс итерации сходится к точному решению системы x при лю
бом x(0), т. е. x = lim x(k) .
k→∞
34
Оценка погрешности этого приближенного решения x(k) дается
одной из следующих формул:
xi − xi(k) ≤
α
⋅ max xj(k) − xj(k−1)
1 − α j =1,n
(2.36)
β n (k)
⋅ ∑ xj − xj(k−1) .
1 − β j =1
(2.37)
либо
xi − xi(k) ≤
Иногда за x(0) берут f, хотя наиболее целесообразно в качестве ком
понент x(0) взять приближенные значения неизвестных, полученные
грубой прикидкой либо прямыми методами.
Приведение системы (2.18) к виду (2.32) можно осуществить раз
личными способами. Если диагональные элементы матрицы A не рав
ны нулю, т. е. aii ≠ 0, i = 1, n, то систему Ax = b можно записать в виде
x1 =
1
( b1 − a12x2 − 1 − a1nxn ),
a11
x2 =
1
( b2 − a21x1 − 1 − a2nxn ),
a22
11111111111111
xn =
1
( bn − an1x1 − an2x2 1 − ann−1xn−1 ).
ann
(2.36)
В этом случае элементы матрицы C определяются следующим об
разом:
aij
cij = − (i ≠ j), cii = 0,
aii
и тогда условия (2.34), (2.35) соответственно приобретают вид
n
a
∑ aij
j =1
n
a
∑ aij
i =1
≤ α < 1, i = 1, 2, 1, n;
ii
≤ β < 1, j = 1, 2, 1, n.
ii
35
Эти неравенства будут выполнены, если диагональные элементы
матрицы A удовлетворяют условию
aii > ∑ aij , i = 1, n,
j ≠i
т. е. если модули диагональных коэффициентов для каждого уравне
ния системы больше суммы модулей всех остальных коэффициентов
(не считая свободных членов).
Метод итераций, если он сходится, дает преимущества по сравне
нию с прямыми методами.
1. Если итерации сходятся достаточно быстро, т. е. если для ре
шения системы требуется менее n итераций, то получаем выигрыш
во времени, так как число арифметических операций, необходимых
для одной итерации, пропорционально n2, а общее число арифмети
ческих действий в методе Гаусса, например, пропорционально n3.
2. Погрешности округления в методе итераций сказываются зна
чительно меньше, чем в методе Гаусса. Кроме того, метод итера
ций является самоисправляющимся, т. е. отдельная ошибка, до
пущенная в вычислениях, не отражается на окончательном резуль
тате, так как ошибочное приближение можно рассматривать как
новый начальный вектор. Последнее обстоятельство часто исполь
зуется для уточнения значений неизвестных, полученных методом
Гаусса.
3. Метод итераций становится особенно выгодным при реше
нии систем, у которых значительное число коэффициентов равно
нулю.
4. Процесс итераций приводит к выполнению однообразных опе
раций и сравнительно легко программируется на ЭВМ.
Метод итераций Зейделя. Метод Зейделя является модификаци
ей метода простой итерации. Он заключается в том, что при вычисле
нии (k + 1)го приближения неизвестного xi при i > 1 используются
уже вычисленные ранее соответствующие приближения неизвестных
x1, x2, …, xi–1. Таким образом, для системы (2.32) вычисления по
методу Зейделя ведутся по формулам
x1(k+1) = c11x1(k) + c12x2(k) + 1 + c1n xn(k) + f1,
x2(k+1) = c21x1(k+1) + c22x2(k) + 1 + c2n xn(k) + f2,
11111111111111111
xn(k+1)
36
= cn1x1(k+1) + cn2x2(k+1) + 1 + cnnxn(k) + fn .
(2.37)
Условия сходимости для метода простой итерации остаются
верными и для метода Зейделя. Обычно метод Зейделя дает лучшую
сходимость, чем метод простой итерации. Кроме того, метод Зей
деля может оказаться более удобным при программировании, так
как при вычислении xi(k+1) нет необходимости хранить значения
x1(k), x2(k), 1, xi(−k1).
2.4. Решение систем уравнений с матрицами специального вида
(трехдиагональные матрицы)
Рассмотрим частный случай «разреженной» матрицы – трехдиа
гональную матрицу.
Определение. Квадратная матрица A размером n×n называется k
диагональной (k – конкретное, нечетное, целое, положительное чис
ло), если ее элементы удовлетворяют соотношениям
aij = 0, i − j > (k − 1) 2.
В соответствии с этим определением обычная диагональная мат
рица является однодиагональной, у трехдиагоналольной матрицы
отличны от нуля элементы, стоящие на главной диагонали и на двух
прилежащих к главной диагоналях.
Метод прогонки решения систем алгебраических уравнений. Ме
тод прогонки является модификацией прямого метода Гаусса для
частного случая разреженной матрицы – трехдиагональной матри
цы. Такие системы уравнений возникают при построении сплайн
функций (см. подразд. 4.4), при конечноразностных аппроксима
циях краевых задач для обыкновенных дифференциальных уравне
ний и в уравнениях с частными производными.
Рассмотрим систему n линейных алгебраических уравнений с n
неизвестными, имеющую следующий специальный вид:
b1x1 + c1x2
= d1,
a2 x1 + b2 x2 + c2x3
= d2 ,
⋅
a3 x2 + b3 x3 + c3 x4
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅
⋅
⋅
⋅
⋅
= d3 ,
⋅
⋅
an −1xn −1 + bn −1xn −1 + cn −1xn = d n −1 ,
an xn −1 + bn xn = dn . (2.38)
37
На главной диагонали матрицы этой системы стоят элементы b1,
b2, …, bn, над ней – элементы c1, c2, …, cn, под ней – a1, a2, …, an. При
этом обычно все коэффициенты bi не равны нулю.
Метод прогонки состоит из двух этапов – прямой прогонки (ана
лога прямого хода метода LUразложений) и обратной прогонки (ана
лога обратного хода). Прямая прогонка состоит в том, что каждое
неизвестное xi выражается через xi+1 с помощью прогоночных коэф
фициентов Ai, Bi:
xi = Ai xi+1 + Bi , i = 1, 2, 1, n − 1.
(2.39)
Из первого уравнения системы (2.37) найдем
x1 = −
c1
d
x2 + 1 ⋅
b1
b1
С другой стороны, по формуле (2.39) x1 = A1x2 + B1. Приравнивая
коэффициенты в обоих выражениях для x1, получим
A1 = −
c1
d
, B1 = 1 ⋅
b1
b1
(2.40)
Из второго уравнения системы (2.38) выразим x2 через x3, заме
няя x1 по формуле (2.39):
a2 ( A1x2 + B1 ) + b2x2 + c2x3 = d2.
Отсюда
x2 =
−c2x3 + d2 − a2 B1
= A2x3 + B2,
a2 A1 + b2
где
A2 = −
c2
d −a B
, B2 = 2 2 1 , e2 = a2 A1 + b2 .
e2
e2
Аналогично можно вычислить прогоночные коэффициенты для
любого номера
Ai = −
ci
d −a B
, Bi = i i i−1 , ei = ai Ai−1 + bi , i = 2, 3, 1, n − 1. (2.41)
ei
ei
Обратная прогонка состоит в последовательном вычислении не
известных xi. Сначала нужно найти xn. Для этого воспользуемся
38
выражением (2.39) при i = n – 1 и последним уравнением системы
(2.38). Запишем их:
xn−1 = An−1xn + Bn−1,
an xn−1 + bnxn = dn .
Откуда, исключая xn–1, находим
xn =
dn − an Bn−1
= Bn .
bn + an An−1
Далее, используя формулу (2.39) и выражения для прогоночных
коэффициентов (2.40), (2.41), последовательно вычисляем все неиз
вестные xn−1, xn−2, 1, x1.
При анализе алгоритма метода прогонки надо учитывать возмож
ность деления на нуль в формулах (2.41). Можно показать, что при
выполнении условий
bj ≥ aj + cj ,
причем хотя бы для одного значения j имеет место строгое неравен
ство, деления на нуль не возникает и система (2.39) имеет единствен
ное решение. Приведенное условие преобладания диагональных эле
ментов обеспечивает также устойчивость метода прогонки относи
тельно погрешности округления. Последнее обстоятельство позво
ляет использовать метод прогонки для решения больших систем
уравнений. Заметим, что данное условие устойчивости прогонки яв
ляется достаточным, но не необходимым. В ряде случаев для хорошо
обусловленных систем (2.39) метод прогонки оказывается устойчи
вым при нарушении условия преобладания диагональных элементов.
2.5. Методы решения проблемы собственных значений
и собственных векторов матриц
Собственные значения и векторы матриц. Пусть матрица A – сим
метричная (A = AT) квадратная матрица ∈ R n×n . Если для некоторого
ненулевого вектора q и некоторого числа λ выполнено равенство
Aq = λq,
(2.42)
то λ и q называются соответственно собственным значением (или чис
лом) и собственным вектором матрицы A. В равенстве (2.42) неизве
39
стны и число λ и вектор q. Чтобы найти их, перепишем равенство
(2.42) в виде однородной системы линейных уравнений
( A − λ E n ) q = Оn ,
(2.43)
где En = diag{1, …, 1}; On – нулевой элементвектор ∈ R n. Система
(2.43) имеет нетривиальное (отличное от нуля) решение лишь в том
случае, если
det ( A − λEn ) = 0.
(2.44)
Уравнение (2.44) называется характеристическим уравнением
матрицы A и может быть расписано в виде
( a11 − λ ) x1 + a12x2 + 1 + a1nxn = 0,
a21x1 + ( a22 − λ ) x2 + 1 + a2nxn = 0,
..........................................
an1x1 + an2x2 + 1 + ( ann − λ ) xn = 0.
Определитель системы (2.44) является многочленом степени n
относительно λ:
Pn ( λ ) = λ n + c1λn−1 + 1 + cn−1λ + cn ,
(2.45)
и называется характеристическим многочленом матрицы A. Таким
образом, собственные значения матрицы A – корни характеристи
ческого многочлена: Pn ( λi ) = 0, i = 1, n. Подставляя последовательно
собственные числа λ1, λ2, …, λn в однородную систему (2.43), находим
собственные векторы qi матрицы A:
Aqi = λi qi , i = 1, n.
(2.46)
Считаем для простоты, что все λi , i = 1, n различны. Совокупность
n векторных равенств (2.46) можно записать в виде одного матрич
ного равенства. Введем для этого матрицу Q, столбцы которой есть
собственные векторы матрицы A: Q = {q1, q2, …, qn}, и диагональную
матрицу Λ с собственными значениями {λi }i=1 матрицы A на диаго
нали: Λ = diag{λ1, λ2, …, λn}. Тогда n равенств (2.46) эквивалентны
одному матричному A(q1, q2, …, qn) = (λ1q1, λq2, …, λnqn) или
n
AQ = QΛ
Λ.
40
(2.47)
Если матрица A вещественная симметричная, т. е. AT=A, то соб
ственные векторы, отвечающие собственным числам, ортогональны:
2
(q i , q j ) = δij qi .
Поскольку собственные векторы qi определены с точностью до про
извольного множителя (уравнение (2.42) и слева, и справа можно
умножить на произвольное число α и переобозначить αq через q), то
выберем длину каждого собственного вектора равной единице. Тогда
матрицу Q, имеющую ортогональные столбцы, длины которых рав
ны единице, называют ортогональной: Q–1 = QT. Домножив обе части
матричного равенства (2.47) слева на QT = Q–1, получим так называ
емую диагональную форму матрицы A
Λ = QT AQ.
(2.48)
С другой стороны, домножив (2.47) справа на QT, находим
A = QΛQT,
(2.49)
что позволяет вернуться от диагональной формы матрицы к самой
матрице A.
Для невырожденных матриц (ABC)–1 = C–1B–1A–1, и тогда линей
ную систему Ax = b перепишем:
(
x = QΛ Q T
)
−1
( )
b = QT
−1
Λ−1Q −1b = QΛ−1QT b,
и так как
x = A–1b,
получаем
A −1 = QΛ−1QT
(2.50)
[сравним с формулой (2.49)].
Формула (2.50) является частным случаем формулы
ϕ(A) = Qϕ(Λ )QT ,
(2.51)
определяющей любую функцию ϕ от квадратной матрицы A:
ϕ( Λ ) = diag {ϕ(λ1 ), ϕ(λ2 ),..., ϕ(λn )}.
41
Подытоживая результаты для симметричных матриц с различ
ными ненулевыми собственными значениями, можно сформулиро
вать следующий теоретический факт о спектральном разложении
матрицы.
Справедливо следующее утверждение: для матрицы A, симмет
ричной и имеющей n различных собственных значений λi, имеет мес
то равенство
T
A = λ1q1q1T + λ2q2qT
2 + 1 + λ nq nq n .
(2.52)
По существу, это иная форма записи равенства (2.49).
Задача нахождения всех собственных значений и собственных век
торов числовой матрицы называется полной проблемой собственных
значений. Эта задача в общем случае достаточно сложная. Мы оста
новимся только на некоторых важных частных задачах.
Нахождение максимального по модулю собственного значения.
Для решения частичной проблемы собственных значений, со
стоящей в определении одного или нескольких собственных значе
ний и соответствующих им собственных векторов, обычно использу
ются итерационные методы.
Предположим, что квадратная матрица A порядка n имеет пол
ную систему нормированных собственных векторов e1, e2, …, en, т. е.
Aei = λiei, ||ei|| = ||ei||2 = 〈ei, ei〉1/2 = 1,
(2.53)
где 〈x, y〉 – скалярное произведение в nмерном векторном простран
стве; λi – собственные значения матрицы A, отвечающие собствен
ным векторам ei, i = 1, 2, …, n. Такая полная система собственных
векторов существует у симметричных матриц. Допустим, что
λ1 > λ2 > λ3 > 1 > λ n .
(2.54)
Выберем произвольный вектор x(0) ≠ 0. Имеем
x(0) = c1e1 + c2e2 +…+ cnen,
где c1, c2, …, cn – координаты вектора x(0) в базисе e1, e2, …, en.
Предположим, что с1 ≠ 0.
Последовательно находим векторы
x(k) = Ax(k–1), k = 1, 2, ….
42
(2.55)
Тогда, согласно (2.53):
x(1) = Ax(0) = A ( c1e1 + c2 e2 + ... + cn en ) =
n
(
= ∑ ci λi ei = λ1 c1e1 + η(1)
i =1
)
и вообще:
(
)
x(k) = λ1(k) c1e1 + η(k) ,
k
(2.56)
k
⎛λ ⎞
⎛λ ⎞
где η(k) = c2 ⎜ 2 ⎟ e2 + 1 + cn ⎜ n ⎟ en , причем из (2.54) имеем
λ
⎝ 1⎠
⎝ λ1 ⎠
⎛λ
η( k) = O ⎜ 2
⎜ λ1
⎝
k⎞
⎟ → 0.
⎟ k→∞
⎠
Здесь и ниже в данном подразделе верхний индекс у векторов, на
пример x(k), указывает номер приближения, у скаляров, например,
λk, означает степень, а в качестве нормы векторов принята вторая
норма, т. е. ||x|| = ||x||2 = 〈x, x〉1/2. Далее ограничимся рассмотрением
класса симметричных матриц. Учитывая, что собственные векторы
симметричной матрицы ортогональны, из формулы (2.56) получим
следующее выражение для скалярного произведения:
x(k) , x(k−1) = λ12k−1 c1e1 + η( k) , c1e1 + η(k−1) =
= λ12k−1 ⎡ c12 e1, e1 + c1 e1, η( k−1) + c1 η(k) , e1 + η(k) , η(k−1) ⎤ =
⎣
⎦
= λ12k−1 ⎡c12
⎣
(k)
(k−1)
+ η ,η
⎡
⎛
⎤ = λ12k−1 ⎢c12 + O ⎜ λ2
⎦
⎜λ
⎢
⎝ 1
⎣
2k−1 ⎞ ⎤
⎟⎥,
⎟⎥
⎠⎦
поскольку ||e1|| = 1, вектор η(k) ортогонален вектору e1 при любом k =
= 1, 2,… и
2k−1
⎛λ ⎞
η(k) , η(k−1) = c22 ⎜ 2 ⎟
⎝ λ1 ⎠
2k−1
⎛λ ⎞
+ 1 + cn2 ⎜ n ⎟
⎝ λ1 ⎠
.
43
Аналогично получаем
x
(k−1)
,x
⎡
(k−1)
= λ12k−2 ⎢ c12
⎢
⎣
⎛λ
+ O⎜ 2
⎜λ
⎝ 1
2k−2 ⎞ ⎤
⎟⎥ .
⎟⎥
⎠⎦
Следовательно:
λ1( k)
=
x(k) , x(k−1)
x(k−1) , x(k−1)
x ( k ) = x( k ) , x ( k )
e1(k) =
x (k)
x
(k)
1/2
⎛λ
= λ1 + O ⎜ 2
⎜λ
⎝ 1
2k−1 ⎞
⎛
⎛λ
k
= λ1 ⎜ c1 + O ⎜ 2
⎜
⎜ λ1
⎝
⎝
⎟;
⎟
⎠
k ⎞⎞
⎟ ⎟;
⎟⎟
⎠⎠
= ( signλ1 ) e11 + r (k) ,
⎛λ
где e11 = ( signc1 ) e1, r (k) = O ⎜ 2
⎜ λ1
⎝
k
k
(2.57)
⎞
⎧–1, x < 0,
⎟, signx = ⎨
⎟
⎩ 1, x > 0.
⎠
Таким образом, при условии (2.54) итерационный процесс (2.55)
позволяет найти с любой точностью максимальное по модулю соб
ственное значение λ1 (ибо λ1(k) → λ1 при k → ∞) и соответствующий ему
собственный вектор e1, так как ||r(k)|| → 0 при k → ∞.
Рассмотрим некоторые замечания.
1. Если |λ1| > 1, то ||x(k)|| → ∞, а если |λ1| < 1, то ||x(k)|| → 0 при k → ∞. То
и другое явление при счете на машине нежелательно (в одном случае
произойдет переполнение, в другом – появляется машинный ноль).
Поэтому следует итерации вести по формулам
e1(0) =
x(0)
x(0)
,
x(k) = Ae1( k−1) , λ1(k) = x(k) , e1(k−1) , e1(k) =
x (k)
x (k)
,
не имеющим указанного недостатка и дающим тот же результат (те
же e1(k) , λ1( k) ), что и формулы (2.55), (2.56).
44
2. Если при выборе x(0) будет выполнено c1 = 0, то за счет ошибок
округлений через несколько итераций появится ненулевая компо
нента, отвечающая e1. Тогда с некоторым запаздыванием итераци
онный процесс выйдет на первое собственное значение.
Рассмотрим теперь задачи, решение которых сводится к отыска
нию максимального по модулю собственного значения некоторой
матрицы B = g(A), такой, что это собственное значение соответствует
отыскиваемому собственному значению матрицы A. В качестве та
кой матрицы B часто берется матрица вида
B = a0E + a1A + a2A2 +…+ amAm,
(2.58)
собственные числа λi(B) которой, согласно (2.51), связаны с собствен
ными числами λi(A) матрицы A соотношением
λi ( B ) = a0 + a1λi ( A ) + a2λ2i ( A ) + 1 + amλim ( A ), i = 1, 2, 1, n. (2.59)
Перейдем к решению конкретных задач.
Нахождение минимального и максимального собственного зна
чения.
Пусть найдено наибольшее по модулю собственное значение λ1 мат
рицы A. Возможны два случая.
1. λ1 > 0, т. е. λ1 = λmax(A).
Для нахождения λmin(A) возьмем матрицу
B = λ1E – A,
собственными значениями которой являются числа λi(B) = λ1 – λi(A) ≥
≥ 0 для любого i = 1, 2, …, n. Тогда λmax(B) = λ1 – λmin(A). Причем
λ max ( B ) = max λi ( B ) в силу неотрицательности спектра матрицы B.
1≤i≤n
Определив максимальное по модулю значение λmax(B) итерационным
методом (2.57), находим минимальное собственное значение матри
цы A:
λmin(A) = λ1 – λmax(B).
2. λ1 < 0, т. е. λ1 = λmin(A).
В этом случае, проведя аналогичные рассуждения для матрицы B =
= A – λ1E, максимальное собственное значение матрицы A находим
по формуле
λmax(A) = λ1 + λmax(B).
45
Нахождение расстояния ρ0 от заданной точки λ = λ0 до ближай
шего собственного значения матрицы.
Рассмотрим случай, представляющий наибольший интерес:
min λi ( A ) ≤ λ0 ≤ max λi ( A ).
1≤i≤n
1≤i≤n
Обозначим
ρ0 = min λi ( A ) − λ0 ,
{
1≤i ≤n
}
l = max max λi ( A ) − λ0 , λ0 − min λi ( A ) .
1≤i ≤n
1≤i≤n
(2.60)
Докажем, что
ρ0 = l 1 − λ ( B ),
где λ ( B ) – максимальное по модулю собственное значение матрицы
(
)
2
1
A − λ0E ,
2
l
B =E−
где E – единичная матрица. Поскольку матрица B – матрица вида
(2.58), то, учитывая (2.59) и (2.60), имеем для любого i = 1, 2 ,…, n
λi ( B ) = 1 −
1
2
l
(λ (A) − λ )
0
i
2
≥ 0.
Отсюда
λ(B) = 1 −
(
1 *
λ ( A ) − λ0
l2
)
2
= 1−
ρ20
,
l2
где λ*(A) – ближайшее к λ0 собственное значение матрицы A. Следо
вательно, расстояние от произвольной точки λ0 спектра матрицы A
до ближайшего собственного значения матрицы
ρ0 = l 1 − λ ( B ).
Замечание. Если λ0 = 0, то ρ0 – это расстояние от нуля до ближай
шего собственного значения матрицы A, т. е. ρ0 = |λn|, где λn – мини
мальное по модулю собственное значение матрицы A.
46
3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ВЫСШИХ ПОРЯДКОВ
И ТРАНСЦЕНДЕНТНЫХ УРАВНЕНИЙ
В этом разделе мы рассмотрим некоторые методы численного ре
шения нелинейного скалярного уравнения вида
f(x) = 0,
(3.1)
где f(x) – заданная функция действительного аргумента x, в частно
сти, f(x) может быть многочленом степени n от x.
Будем считать f(x) непрерывной функцией на отрезке [ a, b ] ⊆ D (a,
b могут быть ±∞). Всякое число x* из D такое, что
f(x*) = 0,
называется корнем уравнения (3.1) или нулем функции f(x). Пусть
f(x) – гладкая функция (т. е. функция, имеющая достаточное число
производных). Число x* называется корнем kкратности, если при
x = x* вместе с функцией f(x) в нуль обращаются все ее производные
до (k – 1)й включительно:
f (x* ) = f ′(x* ) = 1 = f (k−1) (x* ) = 0, f (k) (x* ) ≠ 0.
Однократный корень называется простым.
Алгебраическим уравнением с одним неизвестным называется
уравнение вида
xn + an−1xn−1 + 1 + a1x + a0 = 0,
(3.2)
где n – целое неотрицательное число, называемое степенью уравне
ния; a0, a1, …, an–1 – коэффициенты уравнения; x – неизвестное, кото
рое следует определить.
Известно, что алгебраическое уравнение (3.2) степени n имеет ров
но n корней, если любой kкратный корень считать столько раз, ка
кова его кратность.
Для нелинейного уравнения справедливы утверждения:
Пусть f (x) ∈ C[a,b]. Если f(a)⋅f(b) < 0, то между a и b имеется нечет
ное число корней. Если f(a)⋅f(b) > 0, то между a и b существует четное
число корней или их нет вовсе.
1
Если f (x) ∈ C [a,b], f(a)⋅f(b) < 0 и f ′(x) не меняет знака на [a, b], то
между a и b содержится ровно один корень уравнения.
47
Если уравнение f(x) = 0 не является алгебраическим, то уравнение
(3.1) называется трансцендентным. Например:
f (x) ≡ tg(x) −
1
= 0.
x
Трансцендентное уравнение иногда можно свести к алгебраичес
кому нелинейному уравнению.
При отыскании корня уравнения (3.1) решаем две задачи.
Задача 1. Отделение корней, т. е. отыскание достаточно малых
областей, в каждой из которых заключен один и только один корень
уравнения.
Задача 2. Вычисление корней с заданной точностью.
3.1. Устойчивые решения,
число обусловленности функции одной переменной
Решение прикладных задач, как правило, начинается с создания
приемлемых физических и математических моделей. Пусть есть точ
но поставленная задача: найти на промежутке [a, b] корень уравне
ния
f (x) = 0, a ≤ x ≤ b.
(3.3)
И пусть ищутся корни приближенного уравнения
f (x) = 0, a ≤ x ≤ b
(3.4)
с указанием погрешности в исходных данных:
f (x) − f (x) < ε, ε > 0, a ≤ x ≤ b.
(3.5)
Рассмотрим некоторые трудности, возникающие при решении за
дач (3.3)–(3.5).
1. Пусть имеется уравнение
ex = ε
с решением x = lnε, где ε – малое число. Ошибка в исходных данных
может привести к уравнению ex = 0, не имеющему конечного реше
ния, так как не существует конечного значения x, для которого вы
полнено равенство ex = 0.
48
2. Нетрудно убедиться, что уравнение
x2 – 2x + 1 = 0
имеет корни x1* = 1, x2* = 1, а уравнение
x2 – 2x + 1,01 = 0
обладает корнями x1* = 1 + 0,1i, x2* = 1 − 0,1i. В этом случае погрешность
(ε = 0,01) в исходных данных вывела корни уравнения x2 – 2x + 1 = 0
из области действительных чисел и перевела их в область комплекс
ных чисел.
3. Уравнение
1
x
3
− 2x
1
6
+1 = 0
имеет действительный корень x* = 1, а уравнение
1
x
3
− 2x
1
6
+ 0,99 = 0,
отличающееся от предыдущего лишь погрешностью в свободном члене
ε = 10–2, имеет действительные корни x1* = 1,771561, x2* = 0,531441.
Решение уравнения (3.2) называется устойчивым, если малым из
менениям в исходных данных соответствуют малые изменения в ре
шении. Иными словами, если f (x) − f (x) < ε при всех x ∈[a, b], то и
x * − x* ≤ δ ( ε ),
( )
где x * – точное решение уравнения f x * = 0; х* – точное решение
уравнения f x* = 0 ; ε, δ(ε) – сколь угодно малые числа.
Рассмотрим условия, при которых за решение исходной матема
тической модели [уравнение (3.3)] можно принять решение задачи
(3.4), (3.5). Для этого оценим разность уравнений через погрешность
в исходных данных и функцию (3.4). Справедливо равенство
( )
( ) ( )
или, так как f ( x ) ≡ 0 :
(
)
f x * = f x* − f ′ ( θ ) x* − x * , где θ∈ [x * , x* ] ⊂ [a, b]
*
x* − x * = −
( ) при f ′( θ ) ≠ 0.
f x*
f ′( θ)
49
Обозначим через δf разность в точке x * функций f (x * ) − f (x * ) = δf,
т. е. δf – вариация функций точной и приближенной в точке x *. По
лучим сразу, учтя то обстоятельство, что f (x * ) ≡ 0 :
x* − x * =
δf
.
f ′( θ )
Тогда
(
x* − x * ≤ δf ⋅ min f ′ ( x )
a≤x ≤b
)
−1
.
Назовем величину
M=
(
1
= min f ′ ( x )
a≤x ≤b
min f ′ ( x )
a≤x≤b
)
−1
(3.6)
локальным числом обусловленности функции f(x). При этом вели
чина, которая будет определять характер обусловленности задачи
(3.4), (3.5), есть следующая величина:
(
m = δf ⋅ min f ′ ( x )
a≤x≤b
)
−1
= M ⋅ δf .
(3.7)
Чем больше (меньше) m на отрезке [a, b], где находится корень
уравнения, тем хуже (лучше) обусловленность задачи. Так, в нели
нейном уравнении (п. 3)
1
x
3
− 2x
1
6
+ 0,99 = 0, δf = 0,01;
−5 ⎞
1 ⎛ −2
f ′ ( x ) = ⎜ x 3 − x 6 ⎟;
3⎝
⎠
1
1
minf ′(x) = 0 при x = 1. Следовательно, эта функция f (x) ≡ x 3 −2x 6 + 0,99
плохо обусловлена для функции f (x) ≡ x
1
1
3 − 2x 6 +1, а оценка
x* − x* ≤∞.
Хотя эта оценка носит мажорантный (оценочный) характер, она дол
жна настораживать.
Каждый корень обладает своим числом обусловленности, зная
которое можно оценить наследственную погрешность решения.
Займемся первой из поставленных задач при нахождении корней
уравнения f(x) = 0.
50
3.2. Отделение корней
Существуют несколько различных способов отделения корней
уравнения: теорема Штурма, теорема о перемене знака функции f(x)
и ее производной на промежутке изменения f(x), графический способ
и т. п. Часто для отыскания грубых значений корней можно постро
ить график функции y = f(x) и найти абсциссы точек пересечения гра
фика с осью x. Иногда удобней представить уравнение f(x) = 0 в виде
ϕ(x) = ψ(x), а затем построить графики функций y = ϕ(x), y = ψ(x)
и найти абсциссы их точек пересечения, которые приближенно рав
ны значениям корней. Опишем один из способов отделения корней
(графический) на примере.
Пример. Найти корни уравнения
x sin x = 1 или f (x) ≡ x sin x − 1 = 0.
Решение. Представим функцию f(x) в виде sinx = x −1 и применим
графический способ нахождения корней этого уравнения (рис. 3.1).
Корни уравнения симметричны относительно оси 0y, поэтому мы
можем рассмотреть только положительные корни. Значения x0* , x1* ,1
могут быть определены довольно точно, но xn* , n → ∞, мы определить
не сможем. Тем не менее, по графику ясно, что xn* при n >> 1 близки
к nπ. Все эти полученные приближения к корням уравнения мы мо
( ) ( ) ( )
0
0
0
жем взять за начальные значения x0* , x1* , x2* ,1 и приступить
к задаче 2 уточнения корней.
\
[
\
VLQ [
S
[ [ S
[ [
w
Рис. 3.1. Отделение корней
51
3.3. Методы уточнения корней уравнения
После того, как корни отделены или для них получены хорошие
( )
0
начальные приближения xn* , появляется возможность вычислить
их с заданной точностью. Методы, используемые для уточнения кор
( )
ней уравнения f(x) = 0 при заданном xn*
0
на промежутке [a, b],
в котором находится один единственный корень, разнообразны. Это
может быть метод половинного деления; метод секущих; метод Нью
тона и другие итерационные методы; метод наискорейшего спуска
и т. д. Рассмотрим некоторые из перечисленных выше методов.
Метод простых итераций. Для нахождения изолированного кор
ня x* уравнения f(x) = 0, a ≤ x ≤ b заменим исходное уравнение равно
сильным
х = ϕ(x)
(3.8)
и организуем итерационный процесс вида
xn+1 = ϕ(xn), n = 0, 1, 2, …,
(3.9)
где x0 – заданное начальное приближение.
Можно показать, что, если:
1) |ϕ′(x)| ≤ α < 1 при |x – x0| ≤ r;
2) |ϕ(x0) – x0| ≤ r(1 – α),
то начатый с x0 метод простых итераций сходится к единственному
на отрезке [x0 – r, x0 + r] корню уравнения (3.8). Причем справедливы
оценки погрешности ( ∀n = 1, 2, 1,) :
x* − x n ≤
α
xn − xn−1 ,
1− α
x* − x n ≤
αn
x1 − x0 .
1− α
Выбор функции ϕ(x) в итерационном процессе имеет большое зна
чение и может быть осуществлен не единственным образом.
Предположим, что итерационным методом следует найти корни
уравнения x2 – 1 = 0 на отрезке [0,1; 1,1].
Если итерационную схему записать в виде xn+1 = xn + xn2 − 1, то ϕ(x) =
= x2 – 1 + x. При этом ϕ′(x) = 2x + 1. Очевидно, что ϕ′(x) ≥ ϕ′(0,1) = 1,2
на промежутке [0,1; 1,1]. Следовательно, α > 1 и итерационный про
цесс с такой функцией ϕ(x) не сходится.
52
Если вести итерационный процесс по формуле
1
xn+1 = xn − ⎡⎣xn2 − 1⎤⎦ ,
3
14
1
2
< 1,
т. е. с ϕ ( x ) = x − ⎡ x2 − 1⎤ , то ϕ′ ( x ) = 1 − x и α = max ϕ′ ( x ) =
⎦
0,1≤x ≤1,1
15
3⎣
3
следовательно, процесс сходится.
Итерационный метод Ньютона.
Пусть нужно найти изолированный корень x* ∈ [ a, b] уравнения
f(x) = 0. Предположим, что f ( x ) ∈ C[1a,b] . Тогда функцию f(x) в любой
точке отрезка [a, b], в том числе и в точке x = x*, можно представить
через формулу Тейлора
( )
(
((
)
f x* = f ( x0 ) + x* − x0 f ′ ( x0 ) + O x* − x0
) ),
2
где x0 ∈ [ a, b] – заданное начальное приближение. Но f(x*) = 0, и ко
рень x* удовлетворяет уравнению
(
)
((
f ( x0 ) + x* − x0 f ′ ( x0 ) + O x* − x0
) ) = 0,
2
решить которое точно нельзя, так как неизвестно O((x* – x0)2). От
бросив в предыдущем равенстве последнее слагаемое, приходим
к линейному уравнению, решив которое, получим
x* ≈ x0 −
f ( x0 )
.
f ′ ( x0 )
Это приближенное значение корня обозначим за x1 и найдем сле
дующее приближение x2 по формуле x2 = x1 −
f ( x1 )
и т. д. Таким об
f ′ ( x1 )
разом, метод Ньютона дает следующий итерационный процесс:
xn+1 = xn −
f ( xn )
, n = 0, 1, 2, 1
f ′ ( xn )
(3.10)
На рис. 3.2 показана геометрическая интерпретация метода Нью
тона. Уравнение касательной к графику функции у = f(x), проходящей
через точку с координатами (xn, f(xn)), имеет вид у = f′(xn)(x – xn) +
+ f(xn). Преобразованием формулы (3.10) легко убедиться, что xn+1
53
11
14
13
12
11
1
Рис. 3.2. Приближения к корню нелинейного уравнения методом каса
тельных
является абсциссой точки пересечения этой касательной с осью 0x.
Отсюда другое название метода Ньютона – метод касательных.
Пусть выполнены следующие условия:
1) f(a)⋅f(b) < 0;
2) производные f ′(x) и f′′(x) сохраняют знак на [a, b].
Тогда, если в качестве x0 выбран один из концов отрезка [a, b] так,
что f(x0)⋅f′′(x0) > 0, метод Ньютона сходится.
Пример. Применим метод Ньютона для приближенного вычисле
ния квадратного корня из числа a > 0. В этом случае речь идет о при
ближенном решении уравнения x2 – a = 0, т. е. формулу (3.10) следу
ет применить к функции f(x) = x2 – a.
Поскольку f′(x) = 2x, то для последовательных приближенных
значений xn корня x имеем рекуррентную формулу
x2 − a
xn+1 = xn − n
2xn
или
1⎛
a ⎞
xn+1 = ⎜ xn +
⎟.
2⎝
xn ⎠
В качестве начального приближения можно взять положительное
число x0, такое, что x02 > a. При этом выполняется условие f(x0)⋅f′′(x0) >
54
> 0, так как f′′(x) = 2, и обеспечивается знакопостоянство f′(x), т. е.
метод Ньютона будет сходиться. Например, для a = 4 и начальной
точки x0 = 4 в вычислительных экспериментах получены значения x1 =
= 2,5, x2 ≈ 2,05, x3 ≈ 2,0006 и x4 = 2 с погрешностью менее 10–15. Отме
тим, что в различных языках программирования именно эта проце
дура используется в качестве стандартной функции извлечения квад
ратного корня.
Замечание. Кроме сходимости итерационного процесса необходи
ма и эффективность. Чаще всего в качестве критерия эффективности
итерационного процесса выбирают скорость сходимости.
Так, начиная с некоторого n0, скорость сходимости метода Нью
тона становится квадратичной, т. е. для всех n > n0 выполняется
неравенство
2
x* − xn+1 ≤ c ⋅ x* − xn ,
(3.11)
где с > 0 – некоторая постоянная. (На порядок скорости сходимости
указывает показатель степени в правой части неравенства (3.11).
В данном случае его значение равно 2.)
55
4. ЧИСЛЕННЫЕ МЕТОДЫ ПРИБЛИЖЕНИЯ ФУНКЦИЙ
Приближением функции называется процедура замены по опреде
ленному правилу функции f близкой к ней в том или ином смысле
функцией g из некоторого фиксированного множества. В качестве
приближающего множества часто берут подпространство алгебраи
ческих или тригонометрических многочленов (полиномов). Более
гибкий и мощный аппарат приближения получают, рассматривая
обобщенные полиномы
n
g ( x, a0 , a1, 1, an ) = ∑ ai ϕi ( x ),
(4.1)
i =0
где ϕ0, ϕ1,…, ϕn – некоторая система линейно независимых функций,
выбираемая с учетом конкретных условий задачи и требований,
предъявляемых к функции g. В дальнейшем будем считать, что фун
кция g есть полином вида (4.1).
Мерой погрешности приближения является, как правило, рассто
яние (метрика) ρ(f, g) между приближаемой и приближающей функ
циями. Способ оценки расстояния определяется конкретными усло
виями задачи и имеющейся информацией о приближаемой функции
f. Наиболее часто погрешность приближения на интервале [a, b] оце
нивается в равномерной
ρ ( f, g )C = max f ( x ) − g ( x )
(4.2)
[a,b]
и среднеквадратичной
ρ ( f, g ) L
2
⎧⎪ b
⎫⎪
2
= ⎨ ∫ ⎡⎣ f ( x ) − g ( x ) ⎤⎦ dx ⎬
⎩⎪ a
⎭⎪
1
2
(4.3)
метриках, где [a, b] – интервал, на котором строится приближение.
При выборе метода приближения естественно стремление обеспе
чить более высокую точность приближения и одновременно простоту
процедуры построения g. Первое требование связано с вопросом су
ществования полинома наилучшего приближения, т. е. полинома
g*, удовлетворяющего условию
(
)
ρ f, g * ≤ ρ ( f, g ), ∀g
(знак ∀ означает «любой», для «любого»).
56
(4.4)
Для метрик (4.2) и (4.3) сформулированы и доказаны необходи
мые и достаточные условия существования таких полиномов. В дан
ном пособии рассматривается процедура построения полинома наи
лучшего приближения в среднеквадратичной метрике.
Рассмотрим две постановки задачи приближения.
Пусть в отдельных точках x0, x1, …, xn интервала [a, b] заданы
значения функции f(xj); требуется восстановить ее значения при дру
гих x ∈[a, b]. Если параметры a0, a1, …, an определяются из условий
совпадения значений приближаемой и приближающей функций
в точках xj:
g(xj) = f(xj), j = 0, 1, …, n,
(4.5)
то такой способ приближения называют интерполяцией, а точки xj,
j = 0, 1, …, n – узлами интерполяции. Условия (4.5) носят названия
условий интерполяции. (Первоначально задачей интерполяции на
зывали задачу вычисления значения функции f в некоторой внутрен
ней точке x* интервала [x0, xn] с использованием значений функции
f(xj), вычисленных в точках xj, j = 0, 1 ,…, n этого интервала. До сих
пор задачей экстраполяции называют задачу вычисления значения
функции f в точке x*, находящейся вне интервала [x0, xn].)
Предположим, что алгоритм построения приближающей функции
не требует выполнения условия (4.5) и значения функции вычисляют
ся в произвольных точках, тогда говорят об аппроксимации функции.
2. Пусть функция f задана таблично по результатам измерений, т. е.
значения функции f(xj), вычисленные в точках j = 0, 1,…, n, содер
жат ошибки εj. Построение приближающего полинома интерполя
ционным методом с использованием условий совпадения значений
приближаемой и приближающей функций в узлах привело бы к по
вторению имеющихся ошибок.
Практика показала, что полином g*, минимизирующий погреш
ность приближения в среднеквадратичной метрике, значительно луч
ше представляет функцию f. Таким образом, аппроксимирующий
полином g* строится из условия
(
(
ρ f, g * a0* , a1* , 1, an* , x
))
L2
=
min
1
a0 , a1, , an
ρ ( f, g ( a0 , a1, 1, an , x ) ) L .
2
Процедура построения таких приближений носит название мето
да наименьших квадратов.
Прежде чем перейти к изложению численных методов приближе
ния функций, необходимо сказать несколько слов о характере по
грешности, возникающей в ходе приближения. Наибольший вклад
57
в погрешность приближения вносит, как правило, неустранимая
погрешность. Обычно она обусловлена неполной и неточной инфор
мацией о приближаемой функции (неточность значений функции,
недостаточная информация о гладкости, т. е. о поведении между уз
лами сетки и т. д.). Уточнение информации о приближаемой функ
ции служит основой для выбора численного метода, дающего более
точное приближение. Помимо этого анализ неустранимой погрешно
сти часто приводит к существенному упрощению процедуры вычис
лений. Например, если значения функции получаются в результате
измерений, то нецелесообразно выполнять приближения с точнос
тью более высокой, чем относительная погрешность измерений.
4.1. Интерполяция функций
Интерполяционный многочлен Лагранжа. Пусть на интервале
[a, b] задана система интерполяционных узлов, в которых известны
значения функции f(xj). Задача интерполяции алгебраическими мно
гочленами состоит в построении многочлена
Ln ( x ) = a0 + a1x + a2x2 +1 an xn
степени n, значения которого в заданных узлах интерполяции со
впадают со значениями функции f в этих узлах:
Ln ( xj ) = f ( xj ), j = 0, 1, 1, n.
(4.6)
Покажем, что поставленная задача имеет единственное решение.
Для вычисления коэффициентов ai, i = 0, 1, …, n, составим систему
линейных уравнений
a0 + a1xj + a2xj2 + 1an xjn = f ( xj ), j = 0,1, 1, n.
(4.7)
Определитель системы (4.7) имеет вид
1
x0
x02 1 x0n
1
x1
x12 1 x1n
1 1
1 1 1
1
xn2 1 xnn
xn
и называется определителем Вандермонда. Такой определитель от
личен от нуля, если xi ≠ xj при i ≠ j. В этом случае система (4.7) имеет
единственное решение.
58
Интерполяционный многочлен естественно искать в виде линей
n
ной комбинации Ln ( x ) = ∑ Pi ( x ) f ( xi ) многочленов nй степени Pi(x)
i =0
(i = 0, 1, …, n). Удовлетворяя условиям интерполяции (4.6), получа
ем соотношения
n
∑ Pi ( xj ) f ( xi ) = f ( xj ), j = 0, 1, 1, n,
i =0
которые выполняются, если на функции Pi(x) наложить ограничения
⎧0, i ≠ j;
Pi ( xj ) = ⎨
⎩1, i = j, i, j = 0,1, 1, n.
Последние условия означают, что каждый из многочленов Pi(x),
i = 0, 1, …, n, имеет не менее n нулей на интервале [a, b]. Поэтому Pi(x)
можно записать в виде
Pi ( x ) = Ai ( x − x0 )( x − x1 )1( x − xi−1 )( x − xi+1 )1( x − xn ),
где коэффициент Ai находится из условия Pi(xi) = 1, т. е.
1 = Ai ( xi − x0 )( xi − x1 )1( xi − xi−1 )( xi − xi +1 )1( xi − xn ).
Следовательно, искомые функции Pi(x) вычисляются по формулам
Pi ( x ) =
∏ ( x − xj )
j ≠i
∏ ( xi − xj )
.
j ≠i
Отсюда интерполяционный многочлен имеет вид
n
∏ ( x − xj )
i =0
∏ ( xi − xj )
Ln ( x ) = ∑
j ≠i
f ( xi ).
(4.8)
j ≠i
Форма представления интерполяционного многочлена в виде (4.8)
называется формой Лагранжа. Если ввести обозначение
ωn ( x ) = ( x − x0 )( x − x1 )1( x − xn ),
то многочлен Лагранжа можно записать в виде
59
n
Ln ( x ) = ∑ f ( xi )
i =0
ωn ( x )
.
( x − xi ) ω′n ( x )
Схема Эйткена. Для вычисления значений Ln(x) в заданной точ
ке x* обычно используют алгоритм Эйткена. Обозначим f(xj) = yj.
Построим многочлен первой степени L01(x) по формуле
L01 ( x ) =
( x1 − x ) y0 − ( x0 − x ) y1 .
x1 − x0
Очевидно, что L01(x) является интерполяционным многочленом
для узлов x0 и x1: L01(x0) = y0; L01(x1) = y1. Аналогично вычисляются
многочлены первой степени Li(i+1)(x), i = 1, 2, …, n – 1. Далее постро
им интерполяционный многочлен L012(x) второй степени по правилу
L012 ( x ) =
( x2 − x ) L01 ( x ) − ( x0 − x ) L12 ( x ) .
x2 − x0
Нетрудно проверить, что L012(x) удовлетворяет условиям интер
поляции в узлах x0, x1, x2. Аналогично вычисляются остальные
многочлены второй степени Li(i+1)(i+2)(x), i = 1, 2, …, n – 2 и т. д.
Для всей системы узлов имеем интерполяционный многочлен
L011n ( x ) =
( xn − x ) L011n−1 ( x ) − ( x0 − x ) L121n ( x ) .
xn − x0
В общем случае для интерполяционного многочлена на узлах i,
i + 1, …, i + m имеет место рекуррентное соотношение
Li(i+1)1( i+m ) ( x ) =
( xi+m − x ) Li(i+1)1(i+m−1) ( x ) − ( xi − x ) L(i+1)1(i+m) ( x )
xi+m − xi
.
В целом для вычисления значений Ln(x) в точке x* можно выпи
сать следующую вычислительную схему:
xi
yi
xi – x*
x0
x1
x2
x3
x4
…
y0
y1
y2
y3
y4
…
x0
x1
x2
x3
x4
60
– x*
– x*
– x*
– x*
– x*
…
L(i–1)i
L(i–2)(i–1)i
L(i–3)(i–2)(i–1)i
L(i–4)(i–3)(i–2)(i–1)i
…
L01(x*)
L12(x*)
L23(x*)
L34(x*)
……..
L012(x*)
L123(x*)
L234(x*)
……….
L0123(x*)
L1234(x*)
………..
L01234(x*)
……………
…
Форма записи (4.8) интерполяционного многочлена Лагранжа
удобна тем, что явно выражена зависимость от каждого значения
приближаемой функции в узлах интерполяции. Однако с вычисли
тельной точки зрения она менее эффективна по сравнению с другими
формами. Число операций при расчете Ln(x) составляет 4n – 1, в то
время как для интерполяционного многочлена в форме Ньютона,
который будет обсуждаться ниже, оно равно 3n – 2. Важно также,
что при изменении числа интерполяционных узлов xi все коэффици
енты Pi(x) многочлена Лагранжа должны пересчитываться заново.
Прежде чем перейти к выводу интерполяционной формулы Нью
тона, введем понятие конечных и разделенных разностей.
Конечные разности. Рассмотрим равномерную сетку узлов xk =
= x0 + kh, h > 0, k = 0, 1, …, n, в которых вычисляются значения
функции f: f(xk) = yk. Величина
Δyk = yk+1 − yk
называется конечной разностью первого порядка функции f в точке
xk (с шагом h). Далее
Δ2yk = Δ ( Δyk ) = Δyk+1 − Δyk = yk+2 − 2yk+1 + yk
есть конечная разность второго порядка в точке xk. В общем случае
конечная разность nго порядка функции f в точке xk определяется
по рекуррентной формуле
Δ n yk = Δn−1 ( Δyk ) = Δn−1yk+1 − Δn−1yk .
Свойство конечных разностей. Предположим, что функция
f ∈ C[nxk , xk+n ] , тогда существует такая точка ξ∈ ( xk , xk+n ), что
n
Δ n yk = hn f ( ) ( ξ ).
(4.9)
Доказательство свойства проводится с использованием формулы
конечных приращениий Лагранжа [1]. Имеют место два важных на
практике следствия.
1. Если f ∈ C[nx+01, xn+1 ], то при достаточно малом шаге имеет место
приближенное равенство
Δn+1y0
, где ξ∈ [x0 , xn+1 ].
hn+1
Эта конечноразностная оценка используется на практике, если
выражение производной f(n+1)(x) сложно для вычислений или (n + 1)
раз дифференцируемая функция задана таблично.
f(
n+1)
(ξ) ≈
61
2. Конечные разности nго порядка от многочлена степени n по
стоянны и не зависят от k, а конечные разности более высоких по
рядков равны нулю.
Это обстоятельство используется при выборе степени интерполя
ционного многочлена Ньютона для равномерной сетки узлов, по
скольку приблизительное равенство значений конечных разностей
порядка n некоторой константе указывает соответственно на бли
зость приближаемой функции на заданном участке к многочлену этой
степени.
Разделенные разности. Пусть x0, x1, …, xn – произвольная систе
ма попарно несовпадающих узлов. Значения f(xj) называются разде
ленными разностями нулевого порядка. Величины
f ( xi ; xk ) =
f ( xk ) − f ( xi )
xk − xi
называются разделенными разностями первого порядка функции f(x).
Очевидно:
f ( xi ; xk ) = f ( xk ; xi ) =
f ( xi )
xi − xk
+
f ( xk )
xk − xi
.
Разделенная разность nго порядка определяется через разделен
ные разности (n – 1)го порядка по рекуррентной формуле
f ( x0 ; x1; 1; xn ) =
f ( x1; x2 ; 1; xn ) − f ( x0 ; x2; 1; xn−1 )
xn − x0
.
Свойства разделенных разностей:
1. Разделенная разность nго порядка выражается через узловые
значения функции по формуле
f ( xi )
n
f ( x0 ; x1; 1; xn ) = ∑
i =0
( xi − x0 )( xi − x1 )1( xi − xi−1 )( xi − xi+1 )1( xi − xn )
.
Доказательство проводится по индукции.
2. На равномерной сетке узлов xk = x0 + kh, h > 0, k = 0, 1, …, n, для
разделенной разности nго порядка и конечной разности nго поряд
ка имеет место соотношение
f ( x0 ; x1; 1; xn ) =
62
Δ n y0
n ! hn
.
(4.10)
Пусть [a, b] – минимальный интервал, содержащий узлы x0,
x1, …, xn, f ∈ C[na, b] . Тогда существует такая точка ξ∈[a, b], что
f ( x0 ; x1; 1; xn ) =
n
f( ) (ξ)
n!
.
Это следует из соотношений (4.9) и (4.10).
Интерполяционный многочлен в форме Ньютона. Пусть узлы
x0, x1, …, xn произвольны и попарно не совпадают. Тогда алгебраи
ческий многочлен
ln ( x ) = f ( x0 ) + ( x − x0 ) f ( x0 ; x1 ) + ( x − x0 )( x − x1 ) f ( x0 ; x1 ; x2 ) +
1 + ( x − x0 )( x − x1 )1( x − xn −1 ) f ( x0 ; x1 ; 1; xn )
(4.11)
является интерполяционным, т. е.
ln ( xj ) = f ( xj ), j = 0,1, 1, n.
(4.12)
Покажем, что равенства (4.12) выполняются для n = 1:
ln ( x0 ) = f ( x0 );
ln ( x1 ) = f ( x0 ) + ( x1 − x0 )
f ( x1 ) − f ( x0 )
x1 − x0
= f ( x1 ).
Для n > 1 доказательство проводится по индукции.
Многочлен (4.11) называется интерполяционным многочленом
в форме Ньютона. Он тождественно совпадает с интерполяционным
многочленом в форме Лагранжа, т. е.
ln(x) ≡ Ln(x).
Случай равноотстоящих узлов. Пусть xk = x0 + kh, h > 0, k =
= 0, 1, …, n; f(xk) = yk. Тогда, учитывая соотношение (4.10) для ко
нечной и разделенной разности и вводя безразмерную переменную
x − x0
t=
, перепишем многочлен (4.11) в виде
h
ln ( x ) = ln ( x0 + th ) = y0 + t
Δy0
Δ2y0
+ t ( t − 1)
+
1!
2!
1 + t ( t − 1)1( t − n + 1)
Δ n y0
.
n!
(4.13)
63
Этот многочлен называется интерполяционным многочленом Нью
тона интерполирования «вперед». Здесь начальное значение t = 0
расположено в крайнем слева узле x0. Многочлен (4.13) использует
ся при вычислениях в начале таблиц и для экстраполяции значений
f(x) в точках, лежащих левее x0 (при t < 0).
Интерполяционный многочлен с узлами x–n, x–(n–1), …, x–1, x0, где
x–k = x0 – kh, имеет вид
ln ( x ) = ln ( x0 + th ) = y0 + t
Δy−1
Δ2y−2
+ t ( t + 1)
+
1!
2!
1 + t ⋅ ( t + 1)1( t + n − 1)
Δn y−n
n!
(4.14)
и называется интерполяционным многочленом Ньютона для интер
поляции «назад». Для этого многочлена начало отсчета t = 0 распо
ложено в крайнем правом узле x0. Он используется при интерполяции
в конце таблицы (t < 0) и для экстраполяции правее точки x0 (t > 0).
Чаще всего интерполяционные формулы (4.13) и (4.14) для равно
мерной сетки узлов применяют при t ∈ [−1,1].
Отметим, что для одной и той же системы интерполяционных уз
лов формулы (4.11), (4.13), (4.14) описывают один и тот же интер
поляционный многочлен, ранее записанный в форме Лагранжа.
Оценка погрешности интерполяции. Погрешность приближения
функции интерполяционным многочленом Лагранжа (4.8) описыва
ется разностью
f(x) – Ln(x).
Оценим погрешность интерполяции в точке x*. Будем предпола
гать, что f ∈ C[na+,b1] . Рассмотрим вспомогательную функцию
ϕ ( x ) = f ( x ) − Ln ( x ) − Kωn ( x ),
(4.15)
где ωn(x) = (x – x0)…(x – xn); K – постоянная. Очевидно, что ϕ(xk) = 0,
k = 0, 1, …, n. Выберем K так, чтобы ϕ(x*) = 0, x* ≠ xk, k = 0, 1, …, n.
Имеем
K=
( )
( ).
(x )
f x* − Ln x*
ωn
*
Тогда функция ϕ(x) обращается в нуль по крайней мере в (n + 2)
точках интервала [a, b]: x*, x0, x1, …, xn. Согласно теореме Ролля [1],
64
ϕ′(x) имеет не менее (n + 1) нулей на [a, b], ϕ′′(x) – не менее n нулей
и т. д. Для функции ϕ(n+1)(x) существует по крайней мере одна точка
ξ∈[a, b], в которой ϕ(n+1)(ξ) = 0. Дифференцируя равенство (4.15)
(n + 1) раз, получаем
ϕ(n+1) ( x ) = f (n+1) ( x ) − K ( n + 1) !.
Полагая x = ξ, имеем K = −
f (n+1) ( ξ )
( n + 1) !
. Из (4.15) следует, что по
грешность приближения в точке x* можно представить в виде
( )
( )
f x∗ − Ln x∗ = −
f (n+1) ( ξ )
( n + 1) !
( )
ωn x∗ .
Отсюда получаем выражение для оценки погрешности приближе
ния в точке x*
( )
∗
( ) ≤ max
f x − Ln x
∗
[a,b]
f (n+1) ( x )
( n + 1) !
( )
ωn x∗
и оценки погрешности интерполяции на интервале
max f ( x ) − Ln ( x ) ≤ max
[a,b]
f (n+1) ( x )
[a,b]
( n + 1) !
⋅ max ωn ( x ) .
[a,b]
(4.16)
Поскольку многочлены Лагранжа и Ньютона отличаются только
формой записи, оценка погрешности (4.16) верна для всех формул.
Величину |ωn(x)|, входящую в оценку погрешности (4.16), можно
минимизировать за счет выбора узлов интерполяции. Эта задача ре
шается с помощью многочленов Чебышева: в качестве узлов интер
поляции на интервале [a, b] выбираются корни многочлена Чебыше
ва степени (n + 1) [2], т. е. узлы
xk =
(2k + 1) π, k = 0,1, 1, n.
a+b b−a
+
cos
2
2
2 ( n + 1)
При таком выборе узлов интерполяции оценка (4.16) принимает вид
max f ( x ) − Ln ( x ) ≤ max
[a,b]
[a,b ]
f (n+1) ( x ) ( b − a )n+1
( n + 1) !
22n+1
.
65
О сходимости интерполяционного процесса. На практике часто
возникает вопрос: можно ли, выбрав достаточное количество узлов
интерполяции, приблизить интерполируемую функцию с заданной
степенью точности? Ответ в общем случае отрицателен.
Рассмотрим на интервале [a, b] последовательность систем узлов
интерполяции
{x }, {x
(0)
0
(1)
0 ,
} {
}
n
x1(1) , 1, x0(n), x1(n), 1, xn( ) ,1
(4.17)
Сопоставим последовательности (4.17) последовательность ин
терполяционных многочленов функции f(x)
L0(x), L1(x), …, Ln(x)
{
такую, что многочлен Li(x) соответствует системе узлов x0(i), x1(i), 1,
xi(i)
}. Интерполяционный процесс для функции f(x) называется схо
дящимся, если
Ln ( x ) − f ( x ) → 0.
n→∞
Процесс сходится равномерно, если
max Ln ( x ) − f ( x ) → 0.
[a,b]
n→∞
Сходимость и расходимость интерполяционного процесса зави
сит как от выбора последовательности систем узлов (4.17), так и от
гладкости функции f(x). Например, для функции f(x) = |x|, рассмат
риваемой на интервале [–1, 1], процесс интерполяции не сходится.
На рис. 4.1 представлен график многочлена L9(x) при x ∈[0,1], пост
роенного на интервале [–1, 1] по равноотстоящим узлам.
Интерполяция сплайнами. С увеличением числа узлов интерпо
ляции и соответственно вычисленных в них значениях функции f(x)
естественно желание использовать дополнительную информацию
о поведении приближаемой функции и получить более качественное
приближение. Если число узлов превышает 5–7, то решение постав
ленной задачи посредством интерполяционных многочленов, как
правило, не приводит к успеху. Это объясняется как сильным накоп
лением ошибок при вычислении интерполяционных многочленов
высоких степеней, так и отсутствием сходимости интерполяционно
го процесса в целом. В этом случае целесообразно применять кусоч
66
23145 31314
1
2
1
1
Рис. 4.1. Интерполяция функции |x| многочленом 9й степени
номногочленную интерполяцию, используя многочлены невысоких
степеней (первой, второй, третьей). Наиболее часто применяют ин
терполяцию сплайнами.
Пусть на интервале [a, b] задана непрерывная функция f(x). Рас
смотрим набор узлов a < x0 < x1 < … < xN = b и обозначим f(xi) = fi, i = 0,
1, …, N. Сплайном Sn,ν(x) порядка n дефекта ν, соответствующим
данной функции f(x) и узлам x0, x1, …, xN, называется функция, удов
летворяющая следующим условиям:
а) на произвольном подынтервале [xi–1, xi], i = 1,…, N, функция
Sn,ν(x) является многочленом степени не выше n;
б) функция Sn,ν(x) имеет m непрерывных производных на [a, b]
(ν = n – m);
в) Sn,ν(xi) = f(xi), i = 0, 1, …, N.
Последнее условие является условием интерполирования и опре
деляет название сплайна – интерполяционный сплайн (рис. 4.2).
Примером простейшего сплайна является ломаная. Степень такого
сплайна равна единице, дефект также равен единице. Наиболее ши
роко применяются сплайны третьей степени, имеющие на интервале
[a, b] одну или две непрерывные производные (дефекта 2 или 1 соот
ветственно). Сплайн дефекта 1 называют кубическим.
Рассмотрим процедуру построения интерполяционного кубичес
кого сплайна, удовлетворяющего условиям а–в, и докажем его суще
ствование и единственность.
На каждом из отрезков [xi–1, xi], i = 1, …, N будем искать сплайн
S3,1(x) = Si(x) в виде многочлена третьей степени
67
I[L6QN [
[ D
[
[
[
[ 1w
[1 E
[
Рис. 4.2. Интерполяционный сплайн Sn,ν(x)
Si(x) = ai + bi(x – xi–1) + ci(x – xi–1)2 + di(x – xi–1)3,
(4.18)
где ai, bi, ci, di – коэффициенты, подлежащие определению (i = 1, …, N).
Построим систему алгебраических уравнений для вычисления ис
комых коэффициентов. Из выражения (4.18) для iго звена сплайна
получаем
ai = Si ( xi−1 ), bi = Si′ ( xi −1 ),
2ci = Si′′( xi−1 ), 6di = Si′′′( xi−1 ).
Учитывая условия интерполирования в, имеем
ai = fi–1, i = 0, 1, …, N;
(4.19)
ai + bi(xi – xi–1) + ci(xi – xi–1)2 + di(xi – xi–1)3 = fi.
(4.20)
Далее требования непрерывности первой и второй производных
сплайна S3,1(x) приводят к соотношениям
Si′ ( xi ) = Si′+1 ( xi );
Si′′( xi ) = Si′′+1 ( xi ), i = 1, 2,..., N − 1.
Подставляя выражение (4.18) для сплайна и обозначая xi – xi–1 =
= hi, получаем
2hi ci + 3hi2di = bi+1 − bi , i = 1, 2, 1, N − 1;
68
(4.21)
3h i di = ci+1 − ci , i = 1, 2, 1, N − 1.
(4.22)
Соотношения (4.19)–(4.22) составляют систему 4N – 2 алгебраиче
ских уравнений относительно 4N неизвестных ai, bi, ci, di, i = 1, …, N.
Два недостающих уравнения обычно задаются в виде ограничений на
значения сплайна и его производных на концах интервала и называ
ются граничными условиями. Наиболее употребительными являют
ся следующие три условия.
1. Если известны значения f′ в точках x0 и xN, то естественно по
ложить
′ ( xN ) = fN′ .
S1′ ( x0 ) = f0′, SN
Подставляя выражение для сплайна, получаем два уравнения
b1 = f0′,
2
bN + 2cN hN + 3dN hN
= fN′ .
(4.23)
Полученный в результате решения системы уравнений (4.19)–
(4.23) сплайн называют фундаментальным интерполяционным ку
бическим сплайном.
2. Если в граничных точках интервала известно значение f′′, то
можно задать граничные условия вида
′′ ( xN ) = fN′′
S1′′ ( x0 ) = f0′′, SN
или
2c1 = f0′′;
2cN + 6dN hN = fN′′ .
3. Так называемая интерполяция естественными сплайнами по
лучается из нулевых граничных условий
f0′′ = fN′′ = 0,
т. е.
2c1 = 0;
(4.24)
2cN + 6dN hN = 0.
(4.25)
Такие граничные условия приводят к простейшей системе уравне
ний для искомых коэффициентов. Однако следует заметить, что ис
69
кусственное наложение этих условий на значения сплайна в гранич
ных точках дает погрешность приближения порядка h2 и замедляет
сходимость процесса интерполяции сплайнами на концах интервала.
Покажем, что получаемые системы уравнений имеют единствен
ное решение. Выберем для простоты граничные условия (4.24) и (4.25).
Из соотношений (4.19) определяются все коэффициенты ai, i = 1, … ,
N. Исключим последовательно из уравнений (4.21), (4.22) и (4.25)
переменные di, bi, i = 1, …, N – 1 и получим систему N уравнений
относительно коэффициентов ci, i = 1, … , N.
Из (4.22) и (4.25) имеем
di =
( ci+1 − ci ) , i = 1, 2, 1, N − 1;
3hi
dN = −
cN
.
3hN
Подставляя значения ai и di в соотношения (4.20), можем выра
зить коэффициенты bi:
bi =
fi − fi−1 hi
− ( ci+1 + 2ci ), i = 1, 2, 1, N − 1;
3
hi
bN =
fN − fN −1 2hN cN
.
−
3
hN
Подстановка полученных выражений для di и bi в (4.21) дает ис
комую систему уравнений:
hi −1ci −1 + 2 ( hi −1 + hi ) ci + hi ci +1 =
⎛f −f
f −f ⎞
= 3 ⎜ i i −1 − i −1 i −2 ⎟, i = 2, 3, 1, N − 1;
hi −1 ⎠
⎝ hi
(4.26)
c1 = 0;
(4.27)
⎛ f −f
⎞
f −f
hN −1cN −1 + 2 ( hN −1 + hN ) cN = 3 ⎜ N N −1 − N −1 N −2 ⎟.
hN −1 ⎠
⎝ hN
(4.28)
Соотношения (4.26)–(4.28) определяют трехдиагональную сис
тему N линейных алгебраических уравнений для N неизвестных.
70
Матрица системы имеет диагональное преобладание, т. е. в каждой
строке модуль диагонального элемента превосходит сумму модулей
остальных элементов. Из теории матриц известно, что матрица та
кого вида невырождена, т.е. система (4.26)–(4.28) имеет единствен
ное решение.
Это доказывает единственность интерполяционного сплайна.
О сходимости процесса интерполяции кубическими сплайнами.
Оценки погрешности интерполяции
f(x) – Sn,ν(x)
зависят от сетки интерполяционных узлов. Для простоты ограни
чимся рассмотрением случая равномерной сетки xi = x0 + ih, h > 0, i =
= 0, 1, …, N. Обозначим через S3,1(x, N) сплайн для f(x), построенный
на этой сетке узлов. Приведем без доказательства следующую теоре
му [3].
4
Теорема. Пусть функция f ∈ C[a, b] и в качестве граничных усло
вий выбраны условия 1 или 2. Тогда имеют место оценки
k
(k)
max f ( ) ( x ) − S3,1
( x, N ) ≤ h4−k max f ( IV ) ( x ) ,
[a,b]
(4.29)
[a,b]
где k = 0, 1, 2.
Из оценок (4.29) следует, что при h → 0 (т. е. при N → ∞) соответ
(k)
ствующие последовательности сплайнов и их производных S3,1 ( x, N )
сходятся соответственно к функциям f(k)(x), k = 0, 1, 2.
4.2. Аппроксимация функций в среднеквадратичной метрике
При построении среднеквадратичных приближений обычно рас
сматривают три пространства: пространство En+1 функций, задан
ных своими табличными значениями на множестве точек x0, x1, …,
xn; пространство C[a,b] непрерывных функций, заданных на интерва
ле [a, b]; и пространство L2(a, b) интегрируемых с квадратом функ
ций на интервале (a, b), т.е. функций, для которых выполняется не
равенство (1.12). Способы введения соответствующих скалярных
произведений, норм и метрик описаны в подразд. 1.2, 1.3, 1.5.
Линейно независимые системы в пространстве со скалярным
произведением. Пусть F – пространство со скалярным произведени
ем. Линейная независимость системы элементов f1, f2, …, fn простран
71
ства F устанавливается при помощи определителя Грама, который
имеет вид
〈 f1, f1 〉 〈 f1, f2 〉 1 〈 f1, fn 〉
G ( f1, f2, 1, fn ) =
〈 f2 , f1 〉 〈 f2 , f2 〉 1 〈 f2, fn 〉
11111111111
.
(4.30)
〈 fn , f1 〉 〈 fn , f2 〉 1 〈 fn , fn 〉
Имеет место следующее утверждение.
Теорема. Система элементов f1, f2, …, fn из пространства F линей
но зависима тогда и только тогда, когда определитель Грама этой
системы равен нулю.
Доказательство. Н е о б х о д и м о с т ь. Пусть система элемен
тов f1, f2, …, fn линейно зависима. Покажем, что определитель G(f1,
f2, …, fn) равен нулю. Из линейной зависимости элементов f1, f2, …, fn
следует, что существуют такие числа α1, α2, …, αn, среди которых есть
отличные от нуля, что
α1f1 + α2f2 +…+ αn fn=0.
(4.31)
Умножая равенство (4.31) справа последовательно на fk, k = 1,
2, …, n, получаем
α1〈f1, fk〉 + α2 〈f2, fk〉 +…+ αn 〈fn, fk〉 = 0, k = 1, 2, …, n.
(4.32)
Рассмотрим соотношение (4.32) как систему линейных уравне
ний относительно неизвестных αi, i = 1, 2, …, n. Учитывая предполо
жение относительно αi, заключаем, что система уравнений (4.32)
имеет нетривиальное (ненулевое) решение. Следовательно, опреде
литель системы (4.32), являющийся определителем Грама: G(f1,
f2, …, fn), равен нулю.
Д о с т а т о ч н о с т ь. Пусть определитель G(f1, f2, …, fn) = 0.
Покажем, что система элементов f1, f2, …, fn линейно зависима. Рас
смотрим систему линейных уравнений
〈f1, fk〉β1 + 〈f2, fk〉β2 +…+ 〈fn, fk〉βn = 0, k = 1, 2, …, n,
(4.33)
относительно β1, β2, …, βn. Поскольку определитель этой системы
равен нулю, то она имеет ненулевое решение, т. е. среди чисел βi, i =
= 1, 2, …, n, имеются отличные от нуля. Умножим каждое kе урав
72
нение системы (4.33) на βk и сложим полученные соотношения.
Имеем
n
n
i =1
k=1
∑ βifi , ∑ βkfk
или
n
n
k=1
k=1
∑ βkfk, ∑ βkfk
=0
= 0.
n
В силу аксиом скалярного произведения
∑ βkfk = 0. Учитывая, что
k=1
среди чисел βk, i = 1, 2, …, k, имеются ненулевые, приходим к заклю
чению о линейной зависимости системы элементов f1, f2, …, fn.
Ортогональные и ортонормированные системы элементов. Два
элемента f и g пространства со скалярным произведением называют
ся ортогональными, если их скалярное произведение равно нулю,
т. е. 〈f, g〉 = 0. Элемент f называется нормированным, если ||f|| = 1.
Конечную или бесконечную систему элементов называют ортонорми
рованной, если любые два ее элемента ортогональны и нормированы.
Ортогональная система элементов линейно независима, так как ее
определитель Грама отличен от нуля.
Из системы линейно независимых элементов всегда можно полу
чить ортогональную систему такую, что ее элементы будут представ
ляться в виде линейной комбинации элементов исходной системы
и наоборот. Опишем метод ортогонализации Грама–Шмидта, кото
рый часто используют для получения ортогональных систем.
Пусть f1, f2, …, fn – система линейно независимых элементов. По
строим систему ортогональных элементов g1, g2, …, gn. В силу линей
ной независимости fi имеем 〈fi, fi〉 ≠ 0, i = 1, 2, …, n. Положим g1 = f1.
Определим элемент g2 = f2 + α1(2) g1; подберем коэффициент α1(2) так,
чтобы элемент g2 был ортогонален элементу g1:
〈 g2, g1 〉 = 〈f2 + α1(2) g1, g1 〉 = 〈f2, g1 〉 + α1(2) 〈 g1, g1 〉 = 0.
〈 f2 , g1 〉
. Следующий элемент g3 представим
〈 g1, g1 〉
также в виде линейной комбинации f3, g1 и g2: g3 = f3 + α1(3) g1 + α2(3) g2 .
Из условий ортогональности 〈g3, g1〉 = 0 и 〈g3, g2〉 = 0 следует, что
Следовательно, α1(2) = −
α1(3) = −
〈 f3 , g1 〉
〈f , g 〉
, α2(3) = − 3 2 .
〈 g1, g1 〉
〈 g2, g2 〉
73
Повторяя эту процедуру, можно получить всю ортогональную си
стему. На nм шаге имеют место формулы
n−1
〈f , g 〉
gn = fn + ∑ αk(n) gk , αk(n) = − n k .
〈
gk , gk 〉
k=1
Замечание. Используя процесс ортогонализации Грама–Шмидта,
можно получить ортонормированную систему элементов h1, h2, …,
hn, если сначала положить h1 =
g1
f
= 1 и далее, на каждом iм шаге,
g1
f1
после вычисления коэффициентов αk(i) , k = 1, 2, 1, i − 1, нормировать
gi, т. е. выбирать hi =
gi
.
gi
Ряд Фурье. Система ортогональных элементов называется пол
ной, если ее нельзя пополнить никаким другим ненулевым элемен
том, ортогональным к элементам данной системы.
Пусть ϕ1, ϕ2, …, ϕn, … – бесконечная ортогональная система эле
ментов функционального пространства со скалярным произведени
ем. Числа αi =
〈 f, ϕi 〉
называют коэффициентами Фурье функции f
〈ϕi , ϕi 〉
по ортогональной системе ϕ1, ϕ2, …, ϕn, …. Функции f можно поста
вить в соответствие ряд
f ∼ α1ϕ1 + α2ϕ2 + …+ αnϕn + …
(4.35)
Этот ряд называется рядом Фурье по ортогональной системе фун
кций ϕ1, ϕ2, …, ϕn, …. Можно показать, что если система {ϕi} полна,
то конечные суммы ряда (4.35)
n
Φ n ( x ) = ∑ αkϕk ( x )
k=1
сходятся к функции f в среднем, т. е.
f − Φn → 0
или
n→∞
ρ ( f, Φn ) L → 0.
2 n →∞
Построение многочлена наилучшего приближения в среднеквад
ратичной метрике. Пусть функция задана на некотором множестве
74
A, дискретном или непрерывном. Выберем систему линейно незави
симых функций ϕ0, ϕ1, …, ϕm, заданных на множестве A. Рассмотрим
задачу приближения функции f на множестве A многочленом Pm вида
m
Pm ( x ) = ∑ ak ϕk ( x ).
k=0
При фиксированной системе функций ϕ0, ϕ1, …, ϕm многочлен Pm
однозначно определяется коэффициентами a0, a1, …, am, т. е. Pm(x) =
= Pm(x, a0, …, am). Многочлен Pm∗ называется многочленом наилуч
шего среднеквадратичного приближения (МНСП) функции f на мно
жестве A, если выполняется соотношение
(
ρ2 f, Pm*
)
L2
= min ρ2 ( f, Pm )L .
1
(4.36)
2
a0 , , am
Задача построения МНСП Pm* сводится к задаче вычисления та
ких коэффициентов многочлена Pm, которые минимизируют квадрат
погрешности приближения (4.36). Обозначим эти коэффициенты
∗
и приведем процедуру их вычисления.
через a0∗ , a1∗, 1, am
Квадрат погрешности приближения ρ2 есть функция коэффициен
тов a0, a1, …, am многочлена Pm. Эта погрешность представляет собой
квадратичную форму
m
m
ρ2 ( f, Pm ) L = 〈 f − Pm , f − Pm 〉 = 〈 f − ∑ ak ϕk , f − ∑ ak ϕk 〉 =
2
k =0
k =0
m
m
m
k=0
k =0
j =0
= 〈f, f 〉 − 2〈f, ∑ ak ϕk 〉 + 〈 ∑ ak ϕk , ∑ aj ϕ j 〉 =
m
m m
k =0
k =0 j =0
= 〈 f, f 〉 − 2 ∑ ak 〈 f, ϕk 〉 + ∑ ∑ ak aj 〈ϕk , ϕ j 〉.
(4.37)
В силу аксиом расстояния формула (4.37) является положительно
определенной. Следовательно, она имеет единственный экстремум –
∗
минимум. Это означает, что искомые коэффициенты a0∗ , a1∗ , 1, am
су
ществуют и единственны. Приравняем к нулю частные производные
ρ2 ( f, Pm ) L по переменным a0, a1, …, am:
2
m
m m
⎤
∂ 2
∂ ⎡
ρ ( f, Pm ) L =
⎢ 〈 f, f 〉 − 2 ∑ ak 〈f, ϕk 〉 + ∑ ∑ ak aj 〈ϕk , ϕ j 〉 ⎥ =
2
∂ai
∂ai ⎣⎢
k =0
k = 0 j =0
⎦⎥
m
= −2〈 f, ϕi 〉 + 2 ∑ ak 〈ϕk , ϕi 〉 = 0, i = 0,1, 1, m.
k =0
75
Получаем систему линейных алгебраических уравнений относи
тельно искомых коэффициентов:
a0 〈ϕ0 , ϕ0 〉 + a1 〈ϕ0 , ϕ1 〉 + 1 + am 〈ϕ0 , ϕm 〉 = 〈 f, ϕ0 〉,
a0 〈ϕ1, ϕ0 〉 + a1 〈ϕ1, ϕ1 〉 + 1 + am 〈ϕ1, ϕm 〉 = 〈 f, ϕ1 〉,
11111111111111111111
a0 〈ϕm , ϕ0 〉 + a1 〈ϕm , ϕ1 〉 + 1 + am 〈ϕm , ϕm 〉 = 〈f, ϕm 〉.
(4.38)
Эта система называется системой нормальных уравнений. Ее оп
ределитель является определителем Грама G(f1, f2, …, fn). В силу ли
нейной независимости системы функций {ϕi }0 он отличен от нуля.
m
Поэтому система (4.38) имеет единственное решение. Из теории фун
кций нескольких переменных следует, что это решение дает мини
мум функции ρ2, т. е. решение системы (4.38) – искомые коэффици
∗
енты a0∗ , a1∗, 1, am
.
Таким образом, если система функций ϕ0, ϕ1, …, ϕm линейно не
зависима, то коэффициенты МНСП функции f могут быть вычисле
ны как решение системы нормальных уравнений (4.38). Описан
ный метод построения МНСП называют методом наименьших квад
ратов.
Построение МНСП для функции f упрощается, если система фун
кций ϕ0, ϕ1, …, ϕm ортогональна. В этом случае матрица системы урав
нений (4.38) становится диагональной, и коэффициенты МНСП вы
числяются по формулам
ai =
〈 f, ϕi 〉
, i = 0,1, 1, m.
〈ϕi , ϕi 〉
(4.39)
Коэффициенты ai совпадают с коэффициентами Фурье для функ
ции f по системе функций ϕ0, ϕ1, …, ϕm, а МНСП совпадает с отрезком
ряда Фурье Φm(x) (4.35):
Pm(x) ≡ Φm(x) = a0ϕ0(x) + a1ϕ1(x) +…+ amϕm(x).
(4.40)
Квадрат погрешности приближения описывается формулой
ρ2 ( f, Pm )L = 〈f, f 〉 − a0 〈f, ϕ0 〉 − 1 − am 〈f, ϕm 〉 =
2
= f
76
2
L2
m
− ∑ ai2 ϕi
i =0
L2
.
Из последней формулы явно видно, что с расширением ортогональ
ной системы функций (ростом m) погрешность, вообще говоря, не
увеличивается.
Часто бывает удобнее при построении МНСП использовать орто
нормированную систему функций h0, h1, …, hm. В этом случае, соот
ветственно, коэффициенты имеют вид
ai = 〈 f, hi 〉, i = 0,1, 1, m
и квадрат погрешности вычисляется по формуле
ρ2 ( f, Pm ) L = f
2
2
L2
m
− ∑ ai2 .
i =0
Замечание. Следует иметь в виду, что близость двух непрерывных
функций в смысле среднеквадратичного отклонения не гарантирует
малости их отклонения на всем интервале приближения. Например,
на рис. 4.3 приведены графики функций g(x) ≡ 2 и функции f(x), от
личающейся от g(x) одним пиком высоты 2 + n и длиной основания,
равной 1/n3.
В этом случае расстояние
ρ ( f, g ) L
2
2
⎧⎪ b
⎫⎪
2
= ⎨ ∫ ⎡⎣ f ( x ) − g ( x ) ⎤⎦ dx ⎬
⎩⎪ a
⎭⎪
1
2
1
=
3n
.
5637826637
12324
1
1
4541
2
3
Рис. 4.3. Близкие в среднеквадратичной метрике на интервале [a, b]
функции
77
Очевидно, ρ2 ( f, g ) L → 0, т. е. функции могут быть сколь угодно
2 n →∞
близки в смысле среднеквадратичной метрики, однако сколь угодно
далеки в смысле равномерной метрики: max f ( x ) − g ( x ) = n → ∞.
[a,b ]
n→∞
Приближение функций алгебраическими многочленами. Часто
для среднеквадратичных приближений используют системы алгеб
раических многочленов. Система степенных функций является про
стейшей и имеет вид
(4.41)
1, x, x2,…, xm.
Такая система линейно независима в пространстве C[a,b] непрерыв
ных на интервале [a, b] функций. В пространстве En+1 функций, таб
n
лично заданных на множестве точек xj , система (4.41) линейно
0
независима при m ≤ n + 1 и зависима в противном случае.
Пример. Построить на интервале [0, 1] линейный МНСП для фун
кции f(x) = ex.
Решение. Приближающий многочлен имеет вид P1(x) = a0 + a1x,
т. е. ϕ0(x) = 1, ϕ1(x) = x. Вычислим коэффициенты системы нормаль
ных уравнений по формулам для скалярных произведений (1.6)
{ }
ϕ0 , ϕ0 = 1; ϕ0 , ϕ1 = ϕ1, ϕ0 = 1 2; ϕ1, ϕ1 = 1 3 ;
f, ϕ0 = e − 1; f, ϕ1 = 1.
Запишем нормальную систему
1
⎧
⎪⎪a0 + 2 a1 = e − 1;
⎨
⎪ 1 a + 1 a = 1.
⎪⎩ 2 0 3 1
Решая систему, получаем коэффициенты
α0 = 4e – 10 ≈ 0,87; α1 = 18 – 6e ≈ 1,69.
Отсюда e ≈ 0,87 + 1,69x на интервале [0, 1].
Замечание. В общем случае при вычислении коэффициентов сис
темы нормальных уравнений можно использовать формулу
x
b
〈ϕk , ϕi 〉 = ∫ xk +i dx =
a
78
bk+i +1 − ak+i +1
.
k + i +1
Эта формула достаточно проста. Однако при вычислении правых
частей системы (4.36) интегралы
b
〈f, ϕi 〉 = ∫ f ( x ) xi dx
a
могут оказаться очень сложными. Поэтому часто переходят к диск
ретному варианту метода наименьших квадратов. А именно:
строят
n
МНСП в пространстве En+1 на множестве узлов xj на интервале
0
[a, b], полагая n + 1 равным от 1,5m до 2m. Построенный многочлен
принимают за аппроксимирующий на всем интервале [a, b]. В этом
случае коэффициенты системы нормальных уравнений вычисляют
по формулам (1.7) и погрешность приближения подсчитывают в сред
неквадратичной метрике (1.9).
Основным достоинством системы степенных функций является ее
простота. Однако решение системы нормальных уравнений представ
ляет самостоятельную задачу линейной алгебры. Проблема услож
няется тем, что при больших n (n > 50) матрица системы (4.38) ста
новится плохо обусловленной. Чтобы избежать указанных трудно
стей, часто используют ортогональные системы многочленов. Рас
смотрим некоторые из них.
Многочлены Лежандра описываются выражением
{ }
Ln ( x ) =
1
dn
2n n ! dxn
(x
2
)
n
−1 .
(4.42)
Можно показать [5], что степень многочлена Ln равна n.
Для многочленов Лежандра имеет место рекуррентное соотношение
(n + 1)Ln+1(x) – (2n + 1) x Ln(x) + n Ln–1(x) = 0.
Учитывая, что L0(x) = 1, L1(x) = x, можно вычислить многочлен
любой степени. Например, второй и третий, соответственно, имеют вид
L2 ( x ) = 3 2 x2 − 1 2; L3 ( x ) = 5 2 x3 − 3 2 x.
Покажем, что многочлены Лежандра образуют ортогональную сис
тему на стандартном отрезке [–1, 1], т. е. выполняется соотношение
1
〈 Ln , Lm 〉 =
∫ Ln ( x ) ⋅ Lm ( x ) dx = 0, n ≠ m.
−1
79
Прежде всего, покажем, что многочлены Ln(x) ортогональны од
ночленам xm при m < n. Записывая условие ортогональности и интег
рируя по частям, получаем
〈 Ln ( x ), xm 〉 =
1
dn
∫ (x
2n n ! dxn
1
n
− 1 xmdx =
−1
n−1
n
m
m −1 d
x
x2 − 1 dx.
∫
n
n−1
2 n ! −1
dx
1
=−
)
2
(
)
Выполняя (m – 1) раз интегрирование по частям и учитывая, что
dk
при k < n выражение
(x
2
−1
Ln
2
L2
)
n
обращается в ноль на концах про
dx k
межутка [–1, 1], получаем 〈Ln(x), xm〉 = 0. Отсюда следует, что много
член Ln(x) ортогонален каждому одночлену xk из Lm(x) (k ≤ m) и, сле
довательно, всему многочлену Lm(x).
Замечание. Многочлены Лежандра можно получить из системы
степенных функций, если применить к ним процесс ортогонализа
ции Грама–Шмидта.
Покажем, что система многочленов Лежандра не является орто
нормированной на интервале [–1, 1]. Для этого вычислим норму мно
гочлена Ln:
= 〈 Ln , Ln 〉.
Интегрируя n раз по частям, имеем
Ln
2
L2
=
( −1)n
1
(2 n ! )
n
2
∫ (x
−1
2
−1
)
n
d2n
dx
2n
(x
2
)
n
− 1 dx.
Вычисляя производную порядка 2n от (x2 – 1)n, получаем
Ln
2
L2
=
( −1)n (2n !) 1
(2 n !)
n
2
∫ (x
−1
2
)
n
− 1 dx =
2
2n + 1
или
Ln
80
L2
=
2
.
2n + 1
(4.43)
Можно показать, что многочлены Лежандра образуют полную си
стему функций. Поэтому для произвольной функции f последователь
ность ее МНСП P0, P1, …, Pm,…, построенных для соответствующих
подсистем многочленов {L0}, {L0, L1}, …, {L0, …, Lm}, сходится в сред
нем при m → ∞ к функции f, т. е. выполняются соотношения
f − Pm
→ 0.
L2 m →∞
Пример. Построить на интервале [–π, π] МНСП для функции f(x) =
= sin(x), используя первые три многочлена Лежандра.
Решение. Сделаем замену переменных x = πt, чтобы перевести от
резок [–π, π] в отрезок [–1, 1]. Тогда f ( t ) = sin πt, t ∈ [ −1, 1]. Обозна
чим через P3 ( t ) МНСП для f ( t ) на [–1, 1]. Вычислим коэффициенты
ai многочлена P3 ( t ) по формулам ai =
1
a0 =
〈f, Li 〉
с учетом нормы (4.43):
〈 Li , Li 〉
1
1
3
3
sinπt dt = 0; a1 = ∫ t sinπt dt = ;
π
2 −∫1
2 −1
1
a2 =
1
a3 =
(
)
5 1 2
3t − 1 sinπt dt = 0;
2 −∫1 2
(
)
7 1 3
7 105
5t − 3t sinπt dt = − 3 .
∫
π π
2 −1 2
Следовательно:
P3 ( t ) =
3
15 ⎛ 21 ⎞ 35 ⎛ 15 ⎞
⎛ 7 105 ⎞
L1 ( t ) + ⎜ − 3 ⎟ L3 ( t ) = ⎜ 2 − 1 ⎟ t + ⎜ 1 − 2 ⎟ t3 .
2π ⎝ π
π
π
π ⎠
π ⎠
⎝
⎠ 2π ⎝
Переходя к отрезку [–π, π] с помощью обратной замены перемен
x
ных t = , получаем
π
P3 ( x ) =
15 ⎛ 21 ⎞
35 ⎛ 15 ⎞
− 1 ⎟ x + 4 ⎜ 1 − 2 ⎟ x3 ≈ 0,8x − 0,09x3 .
2⎜ 2
2π ⎝ π
2π ⎝
π ⎠
⎠
Многочлены Чебышева определяются соотношением
Tn ( x ) = cos[n ⋅ arccos x], n = 0,1, 2, 1.
(4.44)
81
Введя обозначение
t = arccos x, ( cos t = x )
(4.45)
и используя метод математической индукции, можно показать, что
Tn(x) – многочлен степени n, а его старший коэффициент равен
2n–1 (n ≥ 1).
Из тождества
cos ( n + 1) t = 2cos t cos nt − cos ( n − 1) t,
где t определяется заменой (4.45), в соответствии с формулой (4.44)
получаем рекуррентное соотношение для многочленов Чебышева:
Tn+1(x) = 2xTn(x) – Tn–1(x), n = 1, 2, ….
Первые четыре многочлена Чебышева имеют вид
T0(x) = 1; T1(x) = x; T2(x) = 2x2 – 1; T3(x) = 4x3 – 3x.
Покажем, что многочлены Чебышева образуют ортогональную
систему на интервале [–1, 1] со скалярным произведением
1
〈Tn , Tm 〉 = ∫ Tn ( x ) ⋅ Tm ( x )
−1
где
ние
1
1 − x2
1
1 − x2
dx,
– весовая функция. Вычислим это скалярное произведе
1
〈Tn , Tm 〉 =
∫ cos ( n arccoc x ) ⋅ cos (m arccoc x )
−1′
1
1 − x2
dx =
u = arccos x
π
dx = ∫ cos nu ⋅ cos mudu =
=
du =
0
1 − x2
=
82
π
π
1
1
cos ( n + m ) udu + ∫ cos ( n − m ) udu = 0, n ≠ m.
∫
20
20
Покажем, что система многочленов Чебышева не является нор
мированной на интервале [–1, 1]:
2
Tn L
2
1
= 〈Tn , Tn 〉 =
∫ cos ( n arccos x )
2
−1′
π
= ∫ cos2 nudu =
0
u = arccos x
1
1 − x2
π
( n ≥ 1), T0
2
dx =
2
L2
du =
=
dx
1−x
2
= π.
Пример. Для функции f (x) = arccos x построить на интервале
[–1, 1] МНСП, используя первые четыре многочлена Чебышева.
Решение. МНСП ищем в виде
P3 ( x ) = a0T0 ( x ) + a1T1 ( x ) + a2T2 ( x ) + a3T3 ( x ).
Вычисляем коэффициенты ak, k = 0, 1, 2, 3, учитывая выражение
для нормы Tk(x):
1
a0 =
ak =
=
1 arccos x
π
dx = ;
π −∫1 1 − x2
2
u = arccos x
1
2 arccos x ⋅ cos ( k arccos x )
dx =
x
=
d
du =
π −∫1
1 − x2
1 − x2
2
2
2 ( −1) − 1
u cos kudu = − ∫ sin kudu = ⋅
, k = 1, 2, 3.
∫
π0
πk 0
π
k2
π
π
k
4
4
Откуда получаем a1 = − ; a2 = 0; a3 = − .
9π
π
Таким образом:
P3 ( x ) =
π 4
4
π 4⎛2
4 ⎞
− T1 ( x ) − T3 ( x ) = − ⎜ x + x3 ⎟.
2 π
9π
2 π⎝3
9 ⎠
Приближение функций тригонометрическими многочленами.
В ряде случаев для среднеквадратичных приближений используется
система функций вида 1, cos x, sin x, cos2x, sin2x, 1. Легко прове
рить, что такая система тригонометрических функций ортогональна
83
на интервале [–π, π]. Соответствующая ортонормированная система
1
1
:
получается из исходной умножением на множители
и
π
2π
1 cos x sin x cos2x sin2x
,
,
,
,
, 1.
2π
π
π
π
π
Из теории рядов Фурье [9] известно, что обе системы являются
полными тригонометрическими системами в пространстве L2(–π, π).
Следовательно, для заданной функции f из L2(–π, π) аппроксимиру
ющий тригонометрический многочлен
Pn ( x ) =
a0 n
+ ∑ ( ak cos kx + bk sin kx ),
2 k =1
(4.46)
коэффициенты которого вычисляются по формулам
ak =
π
1
∫ f ( x ) cos kxdx, k = 0, 1, 2, 1;
π −π
bk =
π
1
∫ f ( x ) sin kxdx, k = 1, 2, 1,
π −π
(4.47)
является МНСП, т. е. погрешность
ρ2 ( f, Pn ) L =
2
π
∫ ⎡⎣f ( x ) − Pn ( x )⎤⎦
2
dx
−π
минимальна. Многочлен (4.46) представляет собой часть известного
тригонометрического ряда Фурье, формулы (4.47) находятся в пол
ном соответствии с формулами (4.39).
Пример. Получить для функции y = x + 1 на интервале [0, 2]
МНСП с использованием системы функций 1, cosx, sinx, cos2x, sin2x.
Решение. Аппроксимирующий многочлен будем искать в виде
P2 ( x ) =
a0 2
+ ∑ ( ak cos kx + bk sin kx ).
2 k =1
Посредством замены переменных
Y = (y – 2)π, X = (x – 1)π
приведем функцию y = x + 1 к функции Y = X, заданной на интервале
[–π, π]. Обозначим аппроксимирующий многочлен для функции Y(X)
84
24156731415
3
2
1
1
1
2
Рис. 4.4. Приближение функции y = x + 1 многочленом P2(x)
через P2 ( X ). Вычислим по формулам (4.47) коэффициенты ak, bk, k =
= 0, 1, 2:
a0 = a1 = a2 = 0;
b1 =
b2 =
π
π
1
2
∫ x sin xdx = π ∫ x sin xdx =2;
π −π
0
π
π
1
2
x sin2xdx = ∫ x sin2xdx = −1.
∫
π −π
π0
Следовательно, P2 ( X ) = 2sin X − sin2X. Выполняя обратную за
мену переменных, получим искомое приближение (рис. 4.4)
2⎛
1
⎞
P2 ( x ) = 2 − ⎜ sin πx + sin2πx ⎟.
2
π⎝
⎠
85
Библиографический список
1. Березин И. С., Жидков Н. П. Методы вычислений. М.: Физмат
гиз, 1962. 465 с.
2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные мето
ды. М.: Лаборатория Базовых Знаний, 2001. 632 с.
3. Волков Е. А. Численные методы. М.: Наука, 2006. 256 с.
4. Фарафонов В. Г. Численные методы решения систем линейных
и нелинейных уравнений: Текст лекций / ЛИАП. Л., 1992. 34 с.
5. Завьялов Ю. С., Квасов Б. И., Мирошниченко В. Л. Методы
сплайнфункций. М.: Наука, 1980. 350 с.
6. Фаддеев Д. К., Фаддеева В. И. Вычислительные методы линей
ной алгебры. М.: Физматгиз, 1963. 502 с.
7. Дьякова Г. Н., Решетов Л. А., Стрепетов А. В. Вычислительная
математика: Программа и методические указания к выполнению кон
трольных работ / ГУАП. СПб., 2004. 26 с.
8. Линник А. Е. Метод наименьших квадратов. М.: Наука, 1980.
201 с.
9. Демидович Б. П., Марон И. А. Основы вычислительной матема
тики. М.: Наука, 2006. 672 с.
86
Учебное издание
Бутенина Дина Викторовна
Стрепетов Алексей Владимирович
ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА
Учебное пособие
Редактор А. Г. Ларионова
Верстальщик С. В. Барашкова
Сдано в набор 07.06.07. Подписано в печать 26.06.07. Формат 60 × 84 1/16.
Бумага офсетная. Печать офсетная. Усл. печ. л. 5,06. Уч.изд. л. 5,53.
Тираж 200 экз. Заказ №
Редакционноиздательский центр ГУАП
190000, СанктПетербург, Б. Морская ул., 67
87
Документ
Категория
Без категории
Просмотров
1
Размер файла
419 Кб
Теги
mat, vychisl, butenina, strepetov
1/--страниц
Пожаловаться на содержимое документа