В.Н. Матвеев
Методы
вычислений
Министерство образования и науки Российской Федерации
Федеральное агентство по образованию
Ярославский государственный университет им. П.Г. Демидова
В.Н. Матвеев
Методы вычислений
Учебное пособие
Рекомендовано
Научно-методическим советом университета
для студентов, обучающихся по специальности
и направлению Математика
Ярославль 2007
УДК 519.6
ББК В 192.1я73
МЗЗ
Рекомендовано
Редакционно-издательским советом университета
в качестве учебного издания. План 2007 года
Рецензенты:
канд. физ.-мат. наук П.А. Корнилов; кафедра прикладной математики
и вычислительной техники Ярославского государственного технического университе-
та.
Матвеев, В.Н. Методы вычислений: учеб, пособие / В.Н.Матвеев; Яросл. гос.
МЗЗ ун-т. — Ярославль: ЯрГУ 2007. — 156 с.
ISBN 978-5-8397-0523-4
В учебном пособии представлены основные темы курса “Методы вычислений", ко-
торый автор течение ряда лет читает в Ярославском государственном университете.
Большинство тем сопровождаются лабораторными работами, в которых предлагается
реализовать изучаемые в лекционном курсе вычислительные методы и алгоритмы.
Предназначено для студентов, обучающихся по специальности 010100 Математика
(дисциплина "Методы вычислений", блок ОПД) и направлению 510(100 Математика (дис-
циплина "Методы вычислений", блок ОПД), очной формы обучения.
Библиогр.: 13 назв.
УДК 519.6
ББК В 192.1я73
ISBN 978-5-8397-0523-4
© Ярославский государственный
университет. 2007
© В.Н. Матвеев, 2007
Оглавление
Основные понятия .............................5
§ 1.	Точные методы решения динейных систем.
Понятие обусловленности ......................9
§ 2.	Частичная проблема собственных значений.
Итерационные методы ..........................13
§ 3.	Итерационные методы решения линейных систем.18
§ 4.	Итерационные методы решения
нелинейных уравнений и систем ...............28
§ 5.	Интерполяционный многочлен Лагранжа ........35
§ 6.	Минимизация остаточного члена.
Многочлены Чебышева .........................38
§ 7.	Разделенные разности.
Интерполяционный многочлен в форме Ньютона .... 43
§ 8.	Интерполирование с кратными узлами..........47
§ Э.	Ортогональные многочлены и их свойства .....51
§ 10.	Численное дифференцирование ...............58
§11.	Численное интегрирование.
Квадратурные формулы Ньюнона — Котеса .......60
§ 12.	Квадратурные формулы Гаусса ...............67
§ 13.	Наилучшее приближение
в гильбертовом пространстве .................72
§ 14.	Наилучшее равномерное приближение..........77
§ 15.	Численные методы решения задачи Коши
для обыкновенных дифференциальных уравнений
и систем.....................................82
§ 16.	Разностные методы решения краевых задач
для обыкновенных дифференциальных уравнений
2-го порядка. Понятие аппроксимации, устойчивости
и сходимости разностных схем ....................89
§ 17.	Методы решения краевых задач 2-го порядка.
Методы прогонки, метод стрельбы ............92
3
§ 18.	Принцип максимума для разностных уравнений 2-го
порядка. Оценка решений первой краевой задачи ... 97
§ 19.	Разностные уравнения .......................101
§ 20.	Собственные значения простейшего разностного
оператора........................................107
§21.	Разностные методы решения краевых задач для
уравнения теплопроводности ......................111
§ 22.	Аппроксимация, устойчивость и сходимость
разностной схемы с весами для уравнения
теплопроводности.................................116
§ 23.	Разностные схемы для уравнения колебаний струны 123
§ 24.	Разностные схемы для уравнения Пуассона ....125
Лабораторная работа	№1 ..........................129
Лабораторная работа	№2 ..........................131
Лабораторная работа	№3 ..........................135
Лабораторная работа	№4..........................138
Лабораторная работа	№5 ..........................141
Лабораторная работа	№6 ..........................144
Лабораторная работа	>7..........................147
Лабораторная работа	№8 ..........................151
Литература........................................155
4
Основные понятия
Понятие о погрешности. При решении тех или иных прикладных задач
важно не только построить алгоритм решения поставленной задачи, но и
оценить погрешность полученного решения. Погрешность принято разде-
лять на неустранимую, вычислительную и погрешность метода.
Неустранимая погрешность связана с тем, что при составлении математи-
ческой модели, которая описывает исследуемый процесс, приходится обычно
прибегать к некоторым упрощающим предположениям, пренебрегая вели-
чинами, влияние которых предполагается несущественным. Далее, при ре-
альном изучении поставленной задачи входные параметры задачи являются
результатом измерений или оценок, выполненных с определенной погреш-
ностью. Поэтому результат решения поставленной задачи будет получен с
некоторой погрешностью,которую и называют неустранимой.
Вычислительная погрешность связана с тем, что вычисления выполня-
ются с конечным числом разрядов. При арифметических операциях на ком-
пьютере число записывается в форме с плавающей точкой, т.е. в виде АЕВ,
где А - мантисса, а В - порядок. Например, число 103,54 записывается в ви-
де 0,10354Е+03. Для записи мантиссы отводится конечное число разрядов,
определяемое типом компьютера и операционной системы. Но при умно-
жении двух чисел, имеющих мантиссу длины п, в результате получается
число с мантиссой длины больше п. Если для записи мантиссы числа ис-
пользуются к разрядов и длина мантиссы больше к, то число округляется в
соответствии с правилами округления: если значение к+1 разряда больше 5,
то значение к-го разряда увеличивается на 1, в противном случае значения
больших разрядов отбрасываются. Таким образом, в арифметических опе-
рациях значения младшего разряда есть результат округления. Но нетрудно
заметить, что погрешность в k-м разряде сомножителей приводит к тому,
что в произведении с погрешностью может определятся уже к-1 разряд.
Следовательно, при большом числе вычислений погрешность, связанная с
округлением в младшем разряде, может постепенно переходить в старшие
разряды, искажая тем самым результат вычислений. Таким образом, вычис-
лительная погрешность возникает как следствие арифметических операций
с конечной разрядной сеткой.
Наконец, погрешность метода возникает в результате применения алго-
ритма приближенного решения поставленной задачи. Изучению таких алго-
ритмов и оценке погрешностей и посвящен этот курс.
5
Напомним основные понятия, используемые в дальнейшем изложении.
Пространство R будем называть линейным, если в нем определена опе-
рация сложения, ассоциативная и коммутативная, операция умножения на
число и существует нулевой элемент 0, такой, что т + 0 = 0 + ж = ж. Эле-
менты я?1, а?2, ...хп называются линейно независимыми, если равенство
C1Z1 + С1Т2 + ... + спхп = 0
возможно тогда и только тогда, когда щ = С2 = ... = сс = 0.
Пространство называется нормированным, если каждому элементу х ста-
вится в соответствие действительное число ||т||, такое, что
1)	|®|| > 0, причем ||а?|| = 0 тогда и только тогда, когда х = 0,
2)	|&ж|| = |&|||®|| для любого комплексного к,
з)	||®+ 2/11 < Ikll + 1Ы1-
Последнее свойство называется неравенством треугольника. Простран-
ство называется строго нормированным, если неравенство треугольника пре-
вращается в равенство тогда и только тогда, когда у=кх, к > 0.
Пространство называется метрическом, если для любых х и у, принад-
лежащих пространству, определено число, называемое расстоянием p(z,?/),
такое, что р(х,у) = р(у,х), причем если = 0, то х=у. Кроме того,
имеет место неравенство треугольника
р(х,у) < p(x,z) + p(z,y).
В нормированном пространстве р(х,у) = ||х — г/||.
Говорят, что в пространстве Н определено скалярное произведение, если
для любых х и у из Н поставлено в соответствие комплексное число (х,у),
такое, что
1)	(х,у)=(®,2/),
2)	(с1ж + с22/,г) = С1(ж,?/) + с2(®,г),
3)	(®,®) > 0. Если (х,х)=0, то х =fl.
Пространство, в котором определено скалярное произведение, называют
гильбертовым. Если (х,у)=(у,х), то его называют вещественным. В гиль-
бертовом пространстве норма определяется естественным образом: ||т|| =
Норма вектора. Норма матрицы
В Rn для вектора х = (®i,®2, ••,хп)т используются следующие нормы.
6
Первая векторная норма ||ж||1 = max|^|-
i
Первую векторную норму часто обозначают ||х||с.
п
Вторая векторная норма ll^lh = 52 lxd-
г=1
/ п	_____
Третья, или евклидова, норма ||т||з = л/52 xi = у/(ж>ж)-
у г=1
Согласно определению, норма оператора А, действующего из простран-
ства Hi в пространство Н2, есть наименьшая из констант С, для которой
||Ае||2 < С||т||1, где Ц.Ц1 — норма в пространстве Hi, а Ц.Ц2 — норма в
пространстве Н2. Из определения очевидно следует неравенство
IIM2 < ШНН^Н,-
Таким образом,
II А II СНГ» II
илн = sup lull •
I ll^lll
Поэтому норма матрицы определяется в зависимости от того, какая норма
выбрана в пространстве Нп. Покажем, что в случае первой векторной нормы
п
Mill = maxV|ay|.
г {
J=1
В самом деле, пусть первая норма вектора х = (xi, Х2, ---ХпУт достигается
на компоненте хр : ||t||i = max |Tj| = |з?р|. Тогда для вектора у=Ах имеем:
п	п
Mi < max ||^| < (maxV |ay|)kP|.
i	i *	*
j=l	j=l
Если введена вторая векторная норма, то
п
шах
з
Действительно, для вектора у — Ах — (yi,yz, ...,уп') имеем ||г/||2 = |з/1| +
...|г/п|. Следовательно,
п	п	п п	п п
IMI2 = У? 152 aiixj\ — 52 52 i°vikji ~ 52^52 la^Dl:rJl•
i=l j—l	г=1 j=l	j=l i=l
7
n	n
Пусть Q = max laO'l- Тогда ||t/||2 < Q £ W-
J t=l	j=l
В случае третьей векторной нормы
||А||з = утах|А;(А*А)| .
Здесь АДА*А) — собственные значения матрицы А*А.
В самом деле,
||АХ||| = (Ах, Ах) = (А*Ах, х).
Собственные вектора матрицы А* А ортогональны. Пусть они еще и орто-
нормированны: (е», ej) = 6ij. Здесь — символы Кронекера. Разложим
вектор х по собственным векторам е* матрицы А*А.
Ж — Щв} Н- С2С2 “1“ ... “Ь
Тогда
|М2 = с2 + с2 + ...+с2.
Но А*Ах = Aiciei + А2с2е2 + ... + ХпСпеп.
Поэтому
п	п	п
||Аж||з = (Ах, Ах) = (А* Ах, х) = (^А^в;, '^aei) = ^XiC* < max |AJ||;r||3
i=l	i=l	г=1
8
§ 1. Точные методы решения линейных систем.
Понятие обусловленности
Метод Гаусса. Рассмотрим линейную систему из п уравнений относитель-
но п неизвестных
аих1 4" <212^2 4" а13хз 4" ••• 4" о>\пхп = ^1
а21ж1 4" 022^2 4" ^23^3 4- ... 4“ 0>2пхп = ^2
«31^1 4“ <232^2 4- <233^3 4“ ••• 4“ а3пхп — Ьз	(1-1)
аП1х1 4- й>п2х2 4” ^пЗхз 4” •••4_ ^ппхп ~ Ьп
Пусть ац отлично от нуля. В противном случае перенумеруем переменные
так, чтобы ац не равнялось нулю. Умножим 1-е уравнение на и вычтем
из 2-го уравнения. Продолжим этот процесс, умножая 1-е уравнение на
и вычитая из 3-го уравнения и т.д.
В результате получим систему, в которой, кроме 1-го уравнения, отсут-
ствует переменная х\
ацХ1 4- ai2x2 4- ахз^з + ... + ainrrn = bi
«22^2 + «23^3 4-... 4- а\пхп = Ь\
^32^2 4- а33х3 + ... + а\пхп = bl
ап2х2 4- а„3х3 4- ... 4- а^пхп = b\.
Продолжим этот процесс, умножая 2-е уравнение на и вычитая из i-ro
°22
уравнения. В случае, если a^2 = 0, перенумеруем переменные Xi начиная с
х-2- Это дает возможность понизить порядок системы еще на единицу и т.д.
В результате придем к треугольной системе:
9
<ЗцТ} + <312^2 + <313^3 4" ••• + dlnxn — ^1
«22^2 + «23^3 + ... + а^пХп = bl
4^3 + ••• + alnxn = bl
a™~i„ A-i + а2п%хп = b™~2
n— In— 1 n 1 1 nn n n
о2~^Хп = 6”-1.
Определяя из последнего уравнения хп, затем из предпоследнего xn„i
и т.д., обратным ходом определим все тг-. Наш процесс может остановить-
ся, если на к-м шаге не удастся найти ненулевой элемент в к строке даже
в результате перенумерации переменных от хь до хп. Это будет означать,
что все коэффициенты этой строки равны нулю. Если bk не равно нулю, то
система не имеет решения, если же и bk = 0, то решение не единственно.
Этот алгоритм требует порядка п2 операций и в силу накопления вычисли-
тельной погрешности неприменим к системам большой размерности. Суще-
ствуют модельные примеры систем размерности >20 с заранее известным
решением, такие, что применение этого алгоритма приводило к результату,
далекому от решения. Особенно этот эффект проявляется в случае плохо
обусловленных систем.
Понятие обусловленности систем. Мера обусловленности
Рассмотрим простой пример
4х + 2у = 6
11.7т+ 6.1?/ = 17.8.
Решение этой системы х=у=1. Изменим правую часть 2-го уравнения
4т + 2у = 6
11.7т+ 6.1?/ = 18
Решение системы х=0.6, у=1.8. Замечаем, что небольшое изменение правой
части системы (примерно на процент) приводит к существенному изменению
решения. Продолжим изучение примера и при той же правой части немного
изменим в исходном примере коэффициенты системы.
10
4х + 2у = 6
12ж + 6.1у = 17.8
Здесь решение х=2.5, у=-2. Обычно считается, что небольшое измене-
ние коэффициентов приводит к небольшому изменению решения. Здесь это
не так. Ситуация крайне неприятная, поскольку в практических расчетах
неизвестно, с какой точностью определены коэффициенты системы. Изучим
этот пример подробнее и найдем обратную матрицу для исходной и изме-
ненной матрицы. Для матриц А и В
4	2 \
11.7 6.1 )
6.1
вычислим обратные матрицы
/ 6.1	—2 \
-11.7 4 )
/ 15.25 -5 \
\ -29.25 10 ) ‘
Таким образом, замечаем, что небольшое изменение коэффициентов ис-
ходной матрицы приводит к непропорционально сильному изменению об-
ратной матрицы. Если вычислить нормы, то получим ||А“ 1||i = 15.7, а
В-11|1 = 39.25. Исторически матрицы, у которых небольшое изменение
коэффициентов приводит к сильному изменению обратной матрицы, назы-
вают плохо обусловленными матрицами, а системы с такими матрицами —
плохо обусловленными системами. Для числовой характеристики этого эф-
фекта вводится понятие меры обусловленности.
Пусть заданы системы
Ах = b и Ау — Ь + г].
Определение. Пусть г=у-х. Мерой обусловленности системы назовем ве-
личину
Поскольку Аг = ту, т. е. г = A то, согласно определению нормы опе-
ратора,
11
Так как Ь=Ах, то
ЦЬ|| ИМ
IMI IMI
< им
Поэтому
Величину

7 = И11|а-1ц
называют мерой обусловленности матрицы. Эта величина характеризует
свойства системы только через свойства матрицы. В практике для харак-
теристики меры обусловленности матрицы используют величину
шах|Л(А)|
min |Л(А) | ’
Легко видеть, что 7 < 7, а в случае симметричных матриц и использова-
ния 3-й векторной нормы 7 — 7. В самом деле, если матрица А симметрич-
ная, то ее 3-я норма совпадает с ее максимальным по модулю собственным
значением, а шах|Л(А)| = |minA(A-1)|. В остальных случаях неравенство
7 < 7 следует из того, что любая норма больше (либо равна) максимального
по модулю собственного значения. Последнее утверждение с очевидностью
следует из неравенства ||А|| ||z|| > ||Az||, если вектор х есть собственный век-
тор, соответствующий максимальному по модулю собственному значению.
В разобранном выше примере Ai = 10, Л2 = 0.1 и 7 = 100.
12
§ 2. Частичная проблема собственных значений
Итерационные методы
Под частичной проблемой собственных значений понимается задача отыс-
кания максимального по модулю собственного значения. В силу ее важности
в приложениях, исследования в этой области имеют большую историю как
отечественных, так и зарубежных математиков. Здесь следует упомянуть в
первую очередь А. Н. Крылова, Д. К. Фадеева, Ф. М. Данилевского.
Напомним, что собственными значениями матрицы А называются корни
ее характеристического многочлена, т. е. корни уравнения
	ап — А	012	&1п	
|А - АЕ| =	^21	022 ~ А ...	^2п	= 0.
	^nl	Оп2	• • •	&пп	
Ниже будет рассмотрен алгоритм приближенного нахождения наибольшего
по модулю собственного значения матрицы А. Далее мы будем предполагать,
что собственные значения матрицы А пронумерованы в порядке убывания
их модулей:
|Ai| > |Л2| > ••• > |АП|.
Рассмотрим итерационный алгоритм
X(fc+1> = АХ{к\	,	(2.1)
выбрав в качестве вектора Х^ любой ненулевой вектор. Пусть разложение
вектора Х^ по собственным векторам ek матрицы А имеет вид
Х^ = ciei + С2^2 + ... + спеп.	(2-2)
Тогда
= AjCiei + А^ег + ... + А*спеп.	(2.3)
Для любой фиксированной компоненты вектора Х^°) имеем разложе-
ние
•Е* )	01^1* Ч" С2&2* Т ••• “Ь
Здесь е,* — соответствующая фиксированная компонента собственного век-
тора ег-. Пусть dj = Ci-e;*, т. е.
#*0) = d\ + d2 + ... + dn.	(2.2')
13
Для этой же компоненты х* вектора Х^
— AjCiei* + А2с2е2* + ••• + А*спеп*
Или
= A, di + Л^2 +   • + A‘d„.	(2-3')
Предположим
1.	Максимальное по модулю собственное значение вещественно
и единственно.
Пусть
(t) = (A^>,.YW)
1	(Х<‘>,XW) 	’
Покажем, что > Ai со скоростью порядка O((fyk).
В самом деле,
X(t+D = Л»+1С1е1 + у Аувд = л*+1[С1е1 + О((^)*+* 1)]
п
Х(‘> = Afcjei + 52 A‘Cie,
i=2
Тогда
(Х»+1\Х^) = A?‘+1[c1e1 + O((^)';+1)][C1e1 + O((^)‘)]
(W‘>) ~ A2‘[ciei + O((^)‘)]2
Из последнего выражения видно, что А^ стремится к А,.
Максимальное по модулю собственное значение можно вычислить и по
более простой формуле
т(*+1)
А, «	(2.4')
Напомним, что х^ — любая фиксированная компонента вектора В
самом деле, из (2.3')
г№+1) = д{+1(</1 + O((^)fc+1)
Ai
и
= A?+1(di + O((^)t+1)
Af(di + O((^)
14
При практических вычислениях отмечено, что итерационный процесс не
даст возможности определить Ai, если в качестве начального вектора
выбран вектор, ортогональный собственному вектору, соответствующему Ар
Кроме того, процесс может остановиться, если Ai и А2 близки в пределах вы-
числительной погрешности.
2.	Два собстственных значения равны по модулю и различны по знаку.
Тогда с помощью итерационного процесса (2.1) А2 приближенно можно
вычислить с помощью выражения
~ (X(2t+2)iX(2t))
1 ~ ' 1
или выражения
(2fc+2)
А? «=	(2.5')
Д'*
Напомним, что хк — любая фиксированная компонента вектора Х^к\
Докажем формулу (2.5х).
С учетом разложения (2.3/) при А2 = —Ai
x(2k+2) = A4fc+2[</i +	+ О((^)2/=+2)]№ + d2 + О((£)2*)]
"	A«[d1 + d2 + O((^)2‘)]2
Отсюда замечаем, что при увеличении к это выражение стремится к А2.
Формула (2.5) доказывается аналогично.
3.	Наибольшие по модулю собственные значения образуют простую
комплексную пару. В этом случае в разложении (2.3х) два первых слагае-
мых Ajdi + А^г комплексно сопряжены. Если Ai = геш, di = Re^, то
Ai = ге~га и di = Re~l@. Поэтому
= 2Rrcos(ka 4- (3) + АзС3е3* + Хкс^е^ +х...
Последнее выражение показывает, что наличие наибольших по модулю
комплексных корней приводит к тому, что величина х^ при итерациях бу-
дет сильно колебаться по величине, меняя знак.
Пусть Ai и Ai — корни квадратного уравнения
А2 + рХ + Q = 0.
Тогда
Ai + Аг = —р и А1А2 = д.
15
Коэффициенты р и q можно приближенно определить из следующих со-
ображений. Пусть к настолько велико, что я* «	+ d2Xk. Тогда
:c(fc+1) р хк + q « d1Aifc"1^[Ai + р Ai + q] + ^2[^2 + Р Аг + <?] = О-
Следовательно,
+ р xk + q ж^-1) « 0.	(2.6)
Такое же приближенное равенство будет справедливо для любой другой ком-
поненты z* вектора Х^к\
4fc+1) + Р zk + q 4fc“1} « 0.	(2.6')
Решая (2.6) и (2.6'), получим
X[k-P)z{k+\) _ z(k-i)x{k+i)
Р	_ z(k-l)x(k)
и
x(k)z(k+l) _ z(k)x(k+l)
q ™ x(k-i)z(k) _ z(k-i)x(ky	(2*8)
4.	Наибольшее по модулю собственное значение вещественно и находится
в жордановой клетке 2-го порядка.
В этом случае в качестве базиса возьмем базис, составленный из собствен-
ных векторов и присоединенного вектора, который соответствует собствен-
ному значению Ai. Пусть ei — присоединенный вектор, а в2 — собственный
вектор, соответствующий собственному значению Ар Тогда
Aei = Авх + е2, Ае2 = Ае2, . . . , Аеп = Хеп.
Поэтому
Afcei = Afei + fcA^-1)e2
Ак е2 = А2в2
71 вп — Апвп
Пусть — начальный вектор. Будем предполагать, что проекция этого
вектора на корневое подпространство, которое соответствует собственному
значению Ai, не равна нулю и не является собственным вектором. Пусть это
вектор h. Примем его за первый вектор базиса. Тогда
Хт = АкХт = А{Л + кХ(к-1)в1 + Х3с3е3 + ... +
16
Для любой фиксированной компоненты вектора получим
= drf + d2fcAj-x + ... + d„.A‘.
т- !*+!
Из последнего разложения видно, что отношение по-прежнему стре-
х*
мится к Л1, но медленнее, чем любая геометрическая прогрессия из-за на-
личия множителя к во втором слагаемом. Именно
^Г = А1(1 + О(-)).
Обычно в этом случае сильно влияет вычислительная погрешностей в
такой ситуации следует поступить так же, как и в случае комплексной па-
ры собственных значений, а именно вычислить коэффициенты р и q
квадратного уравнения по формулам (2.6) и (2.7). Если — р « 4g2, то это
свидетельствует о правильности гипотезы о кратном собственном значении
с жордановой клеткой размерности 2.
В заключение заметим, что в случае 1 вектор Х^ стремится к собствен-
ному вектору, соответствующему Лр который следует нормировать на еди-
ницу. В случае разных по знаку, максимальных по модулю собственных
значений вектор Х^ стремится к вектору, который принадлежит линей-
ной комбинации их собственных векторов. В последнем разобранном случае
вектор Х^ стремится к вектору, принадлежащему корневому подпростран-
ству собственного значения Ар В случае комплексной пары максимальных
по модулю собственных значений можно показать, что вектора
- A2X(fc)
И
%(fc+D _ хгХ^
при увеличении к стремятся к собственным векторам, соответствующим соб-
ственным значениям Ai и А2.
Разумеется, описанные алгоритмы не исчерпывают все возможные слу-
чаи в задаче отыскания максимального по модулю собственного значения,
а лишь позволяют приближенно найти максимальное по модулю собствен-
ное значение в наиболее часто встречающихся ситуациях. Простой пример:
Al = 1, Л2)з = ±г показывает, что ни один из рассмотренных итерационных
процессов не позволяет найти наибольшее по модулю собственное значение.
17
§ 3. Итерационные методы решения линейных систем
Итерационные методы решения линейных систем представляют собой
некоторые алгоритмы построения вектора Х^к\ который является прибли-
женным решением системы п линейных уравнений относительно п неизвест-
ных.
АХ = Ь.	(3.1)
В отличие от точных методов, итерационные методы не являются универ-
сальными. Матрица А системы должна удовлетворять некоторым условиям,
чтобы последовательность векторов Х^ сходилась к точному решению си-
стемы. Однако для систем большой размерности применение итерационных
методов предпочтительнее, так как они более устойчивы к вычислительной
погрешности.
Метод простой итерации.
Систему (3.1) преобразуем к виду
X = ВХ 4- д.	(3.2)
Это можно сделать, переписав систему (3.1) в виде -АХ+Ь=0, и затем до-
бавить и отнять вектор X. Иногда систему (3.1) умножают на некоторую
неособенную матрицу Н, а затем добавляют и отнимают вектор X.
Зададим произвольный ненулевой начальный вектор и построим ите-
рационный процесс по формуле:
X(fc+1) = вх(к) + д	(3.3)
Теорема 1. Пусть ||В|| < q < 1. Тогда существует единственное решение
системы (3.2) и итерационный процесс сходится к решению системы с любого
начального приближения.
Доказательство. Докажем единственность. Предположим противное: су-
ществуют два решения Xi и Х%. Тогда их разность R = Х\ —Х2 есть решение
однородной системы
R = BR.	(3.4)
Отсюда ||Я|| =	< ||В||||7?|. Следовательно, ||В||(1 — ||В||) < 0, т. е.
||Я|| < 0. Но ||.R||| > 0. Поэтому |Я|| = 0 и по свойству нормы R=0.
Достаточность. Пусть R^k+^ = Х^к+^ — Х^. Вектор удовлетворяет
однородной системе R^k+^ = BR^k\ Поэтому
18
и
Пусть ||Я(1)Ц = а. Покажем, что последовательность Х^ фундаментальна.
В силу неравенства треугольника
Поэтому - Х^ || <
< (qk+p 1 + qk+p 2 + ... + qk)a = qka(l + q + ...qp T) < qka^ _ .
Поскольку q< 1, последнее выражение может быть меньше любого числа при
к > ко. Последовательность Х^ фундаментальна, и, следовательно, суще-
ствует предел X. Пусть теперь R^ = Х^ — X. Последовательность R^
удовлетворяет однородному уравнению R^k+1^ = BR^. Поэтому Н-Й^Ц <
^!|Д(°)||. Это означает, что последовательность Х^ сходится к X со скоро-
стью qk. Теорема доказана.
Достаточные признаки сходимости метода простой итерации.
Лемма. Пусть собственные значения матрицы В удовлетворяют условию
 AJ < Чч причем Xi = q соответствует жорданова клетка размерности 1. То-
гда существует матрица Н и неособенная матрица D такая, что Н = D~rBD,
причем ||#||1 < q-
Доказательство. Пусть р — q — max |АД. Рассмотрим матрицу В± =
|А;|<д
р~1В. Собственными значениями матрицы Bi будут p-1Aj. Преобразуем
матрицу Bi к нормальной жордановой форме
Hi = D-'BiD = D^p-'BD.
Но
/ Р TAi а12 О
Я1 = I О Р-1А2 «23
Здесь = 0, если Aj - простое собственное значение, или 1, если Xj
соответствует клетка Жордана размерности > 1. Пусть Н = рН^. Тогда
/ Ai р cvi2 О
Н = pD^p^BD = I О Х2 р «23 • •
19
Вычислим первую норму матрицы Н.
Если |Лг| = q, то ckfi,i+i = 0 по условию леммы. Если же |Аг| < Q, то возможно
= 1. В этом случае |Xi| + р = |Aj| + q — max |Л; < q.
|A<|<g
Поэтому ||#||i < q. Лемма доказана.
Теорема 2. Пусть собственные значения матрицы В по абсолютной вели-
чине меньше единицы. Тогда система (3.2) имеет единственное решение и
итерационный процесс (3.3) сходится к решению системы с любого началь-
ного приближения.
Доказательство. Докажем сначала единственность. Пусть X=DY, где D
— некоторая невырожденная матрица. Тогда систему (3.2) можно записать
в виде
Y = D^BDY + D^g.	(3.5)
Пусть матрица D такова, что матрица Н = D~1BD есть нормальная
жорданова форма матрицы В. Если |А<| < q < 1	||Я||1 < 1, то, со-
гласно теореме 1, система (3.5) имеет единственное решение, и это можно
утверждать и для исходной системы в силу невырожденности матрицы D.
Достаточность. Пусть Н — нормальная жорданова форма матрицы В: В =
DHD~\ Но
Вк = DHD-\..DHD~1 = DHkD-\
Если X — решение системы, то для Rk = Х(к^ — X имеем
||й№+1)||1<||О||1||Я‘||1||Д-1)||1<а/,
где а = ||Z>||||Z)_1||. Таким образом, доказано, что Xfc+1) —> X со скоростью
qk.
Сделаем еще несколько замечаний по поводу теоремы 2. Пусть |AJ < 1
i=2,3,...,n, а |Ai| > 1. Тогда, если выбрать начальное приближение так, что
R^ = ei, где ei — собственный вектор, соответствующий Ai, то очевидно,
что итерационный процесс будет расходиться. С другой стороны, если ре-
шение системы принадлежит линейной комбинации векторов 62, ез, ..,еп, а
вектор таков, что = 6262 4- 6363 + ...Спеп, то итерационный процесс
будет сходиться. Поэтому в некоторых учебниках теорема 2'формулируется
следующим образом
Теорема 2'. Для того чтобы итерационный процесс (3.3) сходился с любо-
го начального приближения, необходимо и достаточно, чтобы все собствен-
ные значения матрицы В были по модулю меньше 1. Как мы заметили,
фраза “с любого начального приближения*4 является существенной.
20
Теорема 3. Пусть матрица А такова, что
У |Т - 9 < 1	« = 1,	(3.6)
А. Ы
ИЛИ
р4<?<1 J = l,..n.	(3.7)
|«w I
Тогда метод простой итерации сходится.
Доказательство. Запишем систему (3.1) в виде
«12	«13	. а1п	bi
Xi = 	Х2 +	—Хз + .	• • Н	4	
«11	«11	ац	«11
«21	«23	. ^2п	&2
Х2 =	Z1 +	—Хз+ .	• • Ч	4	
«22	«22	^22	«22
ап1 . «п2	,	, ®пп—1	,
хп —-----Xi 4----®з 4- ... 4-------хп_1 4-----
апп «пп	«пп	«пп
Условие (3.6) означает, что ||B||i < q < 1, а условие (3.7) означает, что
ВЦг < <7 < 1 и в силу теоремы 1 метод простой итерации сходится.
Метод Зейделя. Пусть матрица А такова, что ац 0. Для нахождения
решения системы (3.1)
«п^1 + «12^2 4- 0,13X3 4- ... 4- a-[nxn = bi
«21*^1 + 6X22^2 4" «23^3 4" ... 4" О,2пХп = &2
0>п1Х1 + ап2Х2 4“ О13Х3 4~ ... 4“ «пп*^п — Ьп
о рганизуем итерационный процесс по формулам
21
au^ifc+1) + ^12^2^ + а13хз^ + ... + ainxW = b\
a2ix^+1^ 4- a22^2fc+1) 4-	+ ... + a2nx^ = 62	(3.8)
n ,T(fc+1) -I- a JW) 4- rtl,T(*+1) 4-	4- a r^k+1>> - b
^nl^i * ^n2^2	’ ^13^3	। • • • 1 ^nn^n	— ^п*
Этот итерационный процесс называется методом Зейделя.
Этот же процесс, записанный в виде
р-1	п	1
дД+1) _	^P3_x[k+1) у^ ^рУд.^) °Р
Р	( ft 3	< л. 3	ft 5
J=1 арр	j=p+lUpp	pp
в отечественной литературе иногда называют методом Некрасова.
.Введем в рассмотрение матрицы
( d\i 0	0. . . .	0
_	021 «22	0	. . О
\ О>п1 ®п2	®пЗ	• •	О'пп /
/ 0 012 а13	. . .	ain \
„ _ О	О О2з	. .	d2n
с~ ..................... ..	О’
\ 0	0	0	. .	О /
С помощью матриц В и С итерационный процесс (3.8) можно записать в
виде:
X(fc+1) = ^B~^CXW 4- В-16.
Поэтому метод Зейделя есть метод простой итерации с матрицей, равной
--В-1С, и правой частью В-1Ь . Как показано выше, для сходимости
этого метода с любого начального приближения необходимо и достаточно,
чтобы собственные значения матрицы —В~ХС были по модулю меньше
единицы. Но
det{-B-xC - ХЕ) = det(-B~x)det{C 4- АВ),
и мы приходим к теореме
22
Теорема 4. Для сходимости метода Зейделя с любого начального при-
ближения необходимо и достаточно, чтобы корни уравнения
ацЛ ai2 ... ain
6Z21A 0,22^ ... в2п	_ Q
• • • О’пп^
были по модулю меньше 1.
Докажем еще несколько достаточных признаков.
Теорема 5. Пусть матрица А -- матрица с диагональным преобладанием,
т. е.
П I I
£ И < 9 < 1 i = 1, -п	(3.9)
j=iw. 1 й|
Тогда метод Зейделя сходится.
Доказательство. Пусть X — точное решение ВХ + СХ = Ь, а Х^
— приближение по методу Зейделя: ВХ^к+^ + СХ^ = Ь. Для =
y(fc) _ х имеем:
BR(k+V + CRM = 0.	(3.10)
Покажем, что ||^fc+1^||i	0.
Пусть первая норма R(k+^ = (г[к+1\ г^+1\.ггп+Г))т достигается на ком-
поненте Грк+1\ т. е. ||7?(fc+1)||i = kp*+1Y
Запишем р-е уравнение системы (3.10) в виде
р-1	п
г(к+1) _ _ ^ipr(fe+l) _ QtPr(fc)
i=l аРР	г^р+1 аРР
Отсюда
|г<‘+1)| <а||Л<‘+1>||1 + /3||Я<‘>||1,
Следовательно, так как |rpfe+1^| =	по предположению, получаем
||Л(,;+1,||1 <а||Л('г+1,||1+ЯЛ(Ч111.
23
Отсюда
1 — а
Но (3 < q — а. Поэтому
Таким образом,
||Л(‘+1)||1 < ?1||Л(‘>||1 < - < <7?||Л(0>||1.
Теорема доказана.
Следующий достаточный признак мы докажем при изучении метода ско-
рейшего спуска.
Теорема 6. Пусть матрица А симметричная и положительно определен-
ная. Тогда метод Зейделя сходится к решению системы (3.1).
Напомним, что матрица А называется положительно определенной, если
для любого ненулевого вектора имеет место неравенство
(АХ,Х)>72(Х,Х).
Метод скорейшего спуска. Для вещественной, симметричной и положи-
тельно определенной матрицы А рассмотрим функционал
F(X) = (АХ,Х) -2(6,Х).	(3.11)
Теорема 7. Пусть X — решение системы
АХ = Ь.	(3.12)
Тогда вектор X сообщает функционалу F(X) минимум. Обратно, вектор X,
на котором функционал достигает своего минимума, есть решение системы
(3.12).
Доказательство. Пусть X — произвольный вещественный вектор, X —
решение системы (3.12) и X — X = R. Поскольку b = АХ и в силу сим-
метричности матрицы А
F(X) = (АХ, X) - 2(Ь, X) = (А(Х + R), (X + 2?)) - 2(6, X + R) =
= (АХ, X) + (АХ, R) + (AR, X) + (AR, R) - 2[(АХ, X) + (АХ, R)] =
= (АХ, X) - 2(6, X) + (AR, R) = F(X) + (AR, R).
24
В силу того, что матрица А положительно определенная, F(X) > F(X),
что и доказывает, что X сообщает функционалу F(X) минимум. Докажем
обратное утверждение, что элемент X. на котором F(X) достигает своего
минимума, есть решение системы (3.12). Функционал F(X) есть функция
координат вектора X - переменных a?i, Xi,..., Xi и в точке минимума grad
F(X)=0, т. е. = 0	» = 1,2,...,п. Но
п п	п
z=l j=l	2—1
И
гх т—!	ТЬ	П
dF \	у—Л
= / apjxj + / J aipXj — 2bp.
p	j=i	i=i
Матрица А симметричная: aij = aji, поэтому
OF	J1 >
= 2(^2 apjxj ~ M-
p	j=i
Следовательно,
gradF^X) = 2(AX - b)	(3.13)
и равенство нулю градиента F(X) при некотором X = X означает, что ком-
поненты этого вектора есть решение системы (3.12). Теорема доказана.
Таким образом, задача отыскания решения системы (3.12) сведена к зада-
че отыскания минимума функционала F(X). Рассматриваемый ниже метод
скорейшего спуска есть один из алгоритмов отыскания минимума функци-
онала F(X).
Построим следующий итерационный процесс
X(fc+i) = Х)к) _ skgradF(X^).	(3.14)
Согласно этому алгоритму, для минимизации F(X) от вектора Х^ мы дви-
жемся в направлении — gradF(X^) — направлении наибольшего убывания
функционала — на некоторое расстояние.которое будет определено из усло-
вия минимума функции
Ф(<5) = F(Xlk} - 6kgradF(X^)).
Согласно (3.13) gradF(X) = 2(АХ-Ь). Пусть R^ = b-AXk\ а 25к =
Хк. Процесс (3.14) можно записать в виде
25
X(fc+i) = х^>+ xkR(k).
Определим Xk из условия минимума функции
Ф(Л) = F(X^ + XR™),
Ф(А) = (А(Х^ + XRW), (Х)к) + XRW)) - 2(6, X(fc) + XRW).
Вычислим Ф'(А).
Ф'(А) = (AR(k\ X(fc) + XRW) + (А(Х^ + XR<M), R^) - 2(6, Rw) =
= 2(AX<k\ RW) + X(AR(k\ R^)) - 2(6,7?(fc)) =
= 2((AXW - b, Rw) + X(AR{k\ R^Y).
Ho AX^ — 6 = — R(k\ и из условия Ф'(А) = 0 получаем
х _ (R^,R^)
к~ (ARW,RWy
Соберем воедино формулы алгоритма метода скорейшего спуска для линей-
ных систем:
х(*+1) =X^ + XkR^k\	(3.15)
где
Д<*>)
А‘ (ARW,RWy	(3'16)
В заключение заметим, что метод Зейделя есть метод покоординатной
минимизации функционала F(X). В самом деле, пусть заданы начальные
приближения а;®, ^з, вектора X. Приближение к компоненте xi - ве-
личину выберем так, чтобы F(X) = F(xi,X2, •••,£„) было минимально. В
этом случае
дху
При доказательстве теоремы 7 мы заметили, что
dF(X)
dxj
есть удвоенное j-e уравнение системы (3.12). Следовательно, определя-
ется из уравнения
auz{ + а12х^ + ... + ain^n = ^1-
26
Пусть теперь заданы начальные приближения ж}, ж®,х„. Минимизация
функционала по переменной ж2 приводит к тому, что х\ определяется из
уравнения
а21^} + а22Х2 + ••• + <*2пХп = Ъ2.
Продолжая этот процесс и считая, что компоненты ж}, х£,,, х\_2 найдены,
определим х\ из условия минимума F(X) по переменоой хп, что приводит к
уравнению
«nl^l “Ь ®п2^2	••• ^ппХп = Ьп.
Таким образом, замечаем, что метод Зейделя есть метод покоординат-
ной минимизации функционала F(X), что и доказывает признак сходимости,
сформулированный в теореме 6.
27
§ 4. Итерационные методы решения нелинейных уравнений и систем
В начале этой темы необходимо отметить, что ситуация с методами ре-
шения нелинейных систем существенно отличается от той, которая имеет
место при изучении итерационных методов для линейных систем. В преды-
дущем параграфе мы доказывали различные теоремы о сходимости тех или
иных методов с любого начального приближения. При исследовании итера-
ционных процессов для нелинейных уравнений и систем выбор начального
приближения играет существенную роль, поскольку неудачный выбор на-
чального приближения к решению может привести к неверному результату
или расходимости исследуемого итерационного процесса. При практических
вычислениях начальное приближение к решению определяется стандартны-
ми методами исследования функций или из иных соображений с учетом
специфики задачи.
При изучении этой темы функции, входящие в уравнения и системы, бу-
дут предполагаться непрерывными и нужное число раз дифференцируемы-
ми.
Метод простой итерации. Пусть задана система
/1(^1, ж2,	= О
/2(^1, ^2,-,жп) = О
/п(®ъ ж2, •••> %п) — О,
или в векторной форме
F(X) = 0.	(4.1)
Преобразуем систему к виду
X = G(X)	(4.2)
и рассмотрим итерационный процесс
X(fc+i) = G(X)k)y	(4.3)
Определение. Пусть в полном метрическом пространстве определено отоб-
ражение G(X). Отображение G(X) называется сжимающим, если для
любых Ai и Аг из области определения G(X) имеет место
p^X^X^qp^X^ д<1.
28
Теорема 1 (принцип сжатых отображений). Пусть G(X) — сжимающее
отображение. Тогда уравнение (4.2) имеет единстенное решение и итераци-
онный процесс (4.3) сходится к решению системы (4.2).
Доказательство. Единственность решения доказывается элементарно. Пусть
есть два решения Х\ и Х2 системы (4.2). Имеем
р(Х1,Х2) = р(С(Хх),С(Х2)) < др(Х1,Х2),
т. е. р(ХьХ2)(1 — q) < 0. Следовательно, р(Хг,Х2) = 0.
Достаточность. Докажем, что последовательность Х^ фундаментальна.
Поскольку
p(X{k+1\Xw) < qp(X{k\X{k-1}) < qka,
р(Х^\Х^\
и силу неравенства треугольника
Р(Х^,ЛЯ) < р(Х^к+р\Х{к+р-1}) + р(Х^к+р-1\Х^к+р-2^ +. . .
+p(X(fc+1), Х^) < aqk(l + q + ... + qp~^ = aqk\^.
1 - g
Следовательно, р(Х1к+1\ X^) < € при к > fco- Согласно критерию Коши,
последовательность Х^ имеет предел, и при р -> оо из последнего нера-
венства имеем
р(Х,Х^) < адк—^—.
1 - q
Теорема доказана.
Метод секущих. Этот достаточно простой метод применяется при на-
хождении корней уравнения
/W = 0.
Пусть на концах некоторого отрезка [#i, z2] функция /(ж) принимает значе-
ния разных знаков. Это означает, что внутри отрезка имеется, по крайней
мере, один корень уравнения (4.1). Для его отыскания проведем прямую
через точки Mi(xi, /(^i)) и М2(т2, /(ж2)). Уравнение этой прямой
У ~ /Ы = —^-(Ж) - /01))-
Х2 —
Точка Х3 пересечения этой прямой с осью х равна
хз =
29
Вычислим значения функции f(x) в точке х% и в качестве следующего от-
резка для построения секущей возьмем тот из отрезков [a?i, Ж2] и [373,^2],
на концах которого функция /(а?) принимает значения разных знаков. В
результате приходим к алгоритму
хк+1 = хк- /(ж&)	** _	(4.4)
J\xk) f\xk-l)
Сравним знаки функции /(ж) в точках xk-i и xk+i и в точках хк и ж^-щ.
В качестве точек хк и хк-\ следующего шага выбирается та пара точек,
в которой значение функции f(x) противоположно по знаку. Сходимость
этого алгоритма для непрерывной функции очевидна.
Метод Ньютона. Приведем эвристические соображения, связанные с по-
строением этого метода. Пусть х — решение уравнения
= О
(4.1')
и xq — приближение к ж, принадежащее некоторой области Q : {|ж — ж| < а},
и пусть
X = h + XQ.
Разложим /(ж) = /(жо + ti) по степеням h, ограничиваясь лишь линейными
членами.
/(ж) = /(ж0 + h) = /(ж0) + hf'(xo) 4- O(h2) = 0.
Отсюда
Поэтому
- _ Ж) , /Q/l2\
№>)+О(Л)’
Следовательно, величина
будет отличаться от решения ж уже на величину <2(/г2). Эти соображения
приводят к построению итерационного процесса уточнения начального при-
ближения	।
/(fffe)
+ l — %k	р,/ \
(4.5)
30
Метод Ньютона иногда называют методом касательных, поскольку следую-
щая точка определяется как точка пересечения касательной к f(x) в точке
Xk- В самом деле, уравнение касательной к f(x) в этой точке имеет вид
У ~ fM = f'(xk)(x - хк}
и точка Xk+i есть точка пересечения касательной к f(x) с осью ох.
f&k)
Хк+1~Хк f(Xty
Аналогичные соображения приводят к построению метода Ньютона для
систем.
Пусть задана система уравнений
/1(Я1,Х2,-,хп) = О
f2(x1,x2,...,xn) = 0	(4.6)
fn(x\, 3?2, •••; Xfi) 0.
Или в векторной форме
F(X) = 0-	(4.6')
Пусть вектор X = (xi, х2,..., хп)т есть решение системы и X = Xq4-/i, то
??ТЬ
X = (ж? 4- hi, ж® + h2,xQn 4- hn)T.
Разложим {(х^ 4- h^x® 4- h2,...,x^ 4- hn) по степеням hi, ограничиваясь
линейными членами. Так же, как и для одного уравнения, получим
dxi
C/(Z?2	<^71
-- г° т° т0х. %М-^п),	%(^1-^п), , df^xl-.x*.
' ^^-0+——h>+—
Г?. ?ть определитель матрицы Якоби
hn+O(h2) = 0.
С/Ж72
31
(dfi	ЭЛ	\
dxi	’	'	'	дхп	1
............................. (47)
9fn	dfn	I
dxi	’	‘	*	dxn	/
отличен от нуля. Найдем из этой системы	вектор	h
h = W~1(x0)F(x0) +
Поэтому вектор
X, = Хо - W-1(X0)F(X0)
будет отличаться от Хо на величину порядка O(h2). Эти соображения при-
водят к построению итерационного процесса
Хк+1 =Хк-VV-1(Xt)F(Xt).	(4.8)
В случае скалярного уравнения легко доказывается теорема о сходимости
метода Ньютона.
Теорема 2. Пусть х — решение уравнения (4.1') и пусть в некоторой об-
ласти Qo : |ж —	< а выполнены условия
1. Первая производная отделена от 0:	|//(z)|-1 < Ь.
2. Вторая производная ограничена:	|/"(ж) | < d.
Тогда метод скорейшего спуска сходится в области Qo при
1
а < —.
bd
Доказательство. Рассмотрим выражение
/(ж) - f(xk) - f(xk)(x - хк).
Поскольку f(x) — 0 и так как в силу (4.5) f'(xk)(xk+i — Хк) = —f(xk)
|/(я) - fM - f'(xk)(x - ®fe)| = |/'(rrfc)| |zfc+i - ж|.	(4.9)
С другой стороны, по теореме Лагранжа
\f(x)-f(xk)-fXxk)(x-Xk)\ = |(/'(0-/'(^))(^-^)1 = KfM^-X^^-Xk)
Здесь х < £ < х^ <rj < хк. Следовательно,
|/(ж) - f(xk) - f'(xk)(x - xfc)| < |(/"(^)| |я - Zfc|2. (4.10)
32
Сравнивая (4.9) и (4.10), получаем
|/'Ы1- г| < |i - ^|2.	(4.11)
Пусть с = bd. Неравенство (4.11) означает, что
кл+1 - я| < с\х - xfc|2.	(4.12)
Пусть pk = c|:rfc — ж|. Умножая (4.12) на с, получим pk+i < и
.9^4	2к
Рк+1 < Рк < Рк-1 — •" — Ро •
Следовательно,
к ~ я*+1| < с-1(с к - £о| )2*-
Для сходимости надо, чтобы с к — #о| < 1, т.е. область Qo должна быть
такой, что а < Теорема доказана.
Метод наискорейшего градиентного спуска. Известно достаточо боль-
шое число различных градиентных методов. Мы рассмотрим один из про-
стейших вариантов метода наискорейшего градиентного спуска.
Системе (4.6) поставим в соответствие функцию
Ф(я) = /1(*ь	®п) + ••• + /пкь •••’ *п).	(4.13)
Решение системы (4.6) сообщает этой функции нулевой минимум, и, обратно,
значения вектора X = (zi, #2, • • •, хп)т, в котором функция Ф(я) достигает
значения нулевого минимума, является решением системы (4.6).
Построим следующий итерационный процесс
X(fc+i) = %(Л) _AfcV(fc)?
где = gradV(X№.
Таким образом, поскольку Х^к+1^ — Х^ = —Xkgrad^{X^\ движение от
точки Х^ к точке происходит в направлении наибольшего убывания
функции Ф(Х) на некоторый шаг Хк- Величину Хк следует определить из
условия минимума функции Ф(Х^ — XV^), что опять приводит к решению
нелинейной задачи отыскания корня производной функции
Т'Л) = Ф(Х^ — XV^). В простейшем варианте метода для отыскания Л
приближенно отыскивается корень уравнения
Ф(Л) = 0.
33
Если разложить Ф(А) по степеням Л, пренебрегая квадратами и высшими
степенями Л
Ф(Л)«Ф(0) + ЛФ'(0)),
то
Ф(0)
‘	Ф'л(О)'
Поскольку
Ф(Л) = Ф(а?1 - Ащ,..., хп - Avn),
где
£	Ф'(А) =	= ^rad^X^gradi’^)').
1=1	2 = 1
В итоге приходим к следующему итерационному процессу
X(fc+1) = X(k) _ хку(к),	(4.14)
где
И‘> = дгаЯЩХМ), Ф(Х) = Л2 + /22 + ... + Л2, А* = (^^)}. (4.15)
Следует еще раз подчеркнуть, что при неудачном выборе начального при-
ближения к решению системы построенный процесс может расходиться или
сходиться, но не к решению системы, о чем свидетельствуют допущения при
вычислении А.
Отметим без доказательства, что более точный метод отыскания А как ко-
рень уравнения Ф'(А) = 0 приводит к формуле
х(к+1) = х(к) _ xkW\xM)F(x^,	(4.16)
где
(F(x^), W{x^)W\x^)F(x^))
к~ (W(x^)W'(x^F(x^), W(x^)W\x^)F(x^))'	'
Здесь вектор F = (/i, /2, • • • fn)T, W ~ матрица Якоби, определенная в (4.7),
a W' — транспонированная матрица Якоби.
34
§ 5. Интерполяционный многочлен Лагранжа
Интерполяционный многочлен Лагранжа. Пусть на отрезке [а,Ь] задана
функция /(а?) и пусть внутри этого отрезка заданы п различных точек
Х2,хп. Требуется построить многочлен минимальной степени
Рк(х) = аохк 4- aix^1 + ...ak-ix + а* ,
,1‘
такой, что
= f(xj) j = l,2,...,n.	(5.1)
Теорема 1. Существует единственный многочлен степени п — 1, такой,
что Pn(xi} = f(xi) i = 1,п.
Доказательство. Условие (5.1) означает, что
a0Zi-1 + ala?”-2 + ... + an_i = /‘(яд)
ao^2-1 + aix2~2 + ••• + an-i - /(^2)
1 +	2 + ... + an_i = f(xn).
Эта система из n уравнений относительно п неизвестных ад, а1 , •••> ®п-1-
Определитель этой системы
1	п—2
	Х1	5
п—1	
^2 5	^2 >
zpfl 1	„п-2 7
есть определитель Вандермонда. Он отличен от нуля при Х{ Xj, если
I j. Следовательно, система имеет единственное решение относительно
ад. «1,..., ап-1- Теорема доказана.
Многочлены, удовлетворяющие условию (5.1), называют интерполяцион-
ными многочленами. Лагранж предложил простой способ построения таких
многочленов.
Рассмотрим
35
Многочлен щ,п(х) степени п-1 в точках Xj равен 0 при j 7^ i, а в точке I
Xj равен 1. Следовательно, многочлен
П (* -
п	j = l
=	-------- (5.2)
г=1	П	— xj)
j=l
j^i
степени п-1 удовлетворяет условиям (5.1). Многочлен, определяемый фор-
мулой (5.2), называется интерполяционным многочленом Лагранжа. Введем
в рассмотрение многочлен
wnW = (х - ят1)(ат - хг).'.(х - хп).
№
С его помощью интерполяционный многочлен Лагранжа можно записать
следующим образом
п
Ln(x) =^f(Xi)
i=l
Unfa)
(5-3')
(Ж - X^'^Xi)
В самом деле,
п п
^n(*C) = V П(*С —
г=1 j=l
и при х = Xi в сумме все слагаемые равны нулю, кроме одного, где отсут-
ствует множитель х — Xi.
Величину if>(x) = Ln(x) — f(x} называют остаточным членом интерполя-
ционной формулы
Ln{x) « f(x).
Следующая теорема является важным инструментом для оценки вели-
чины tp(x).	>
Теорема 2. Для любой точки z Xi найдется т^кая точка £, что
---------- /
z) - Ln(z) = —
Доказательство. Рассмотрим функцию
Ф(х) = f{x} - Ln(x) - Аып(х),
36
где А — некоторая константа. Эта функция обращается в нуль в точках
Xi. Выберем константу А так, чтобы эта функция была равна нулю еще в
одной точке z. Тогда
Л -	~ Ln^
^п(^)
С другой стороны, функция обращается в нуль в п точках Xi и точке z, т. е.
в п+1 точке. По теореме Ролля между двумя последовательными нулями
функции есть, по крайней мере, одна точка, где производная этой функ-
ции обращается в нуль. Следовательно, для Ф'(я) имеется не менее п точек
si,^2», где Ф'(£) = 0 обращается в нуль. Продолжая эти рассужде-
ния, замечаем, что существует хотя бы одна точка £, в которой = 0.
Но Ln(x) — многочлен степени п-1, а а?п(я) — многочлен степени п со
старшим коэффициентом 1. Поэтому
- Ап'..
Равенство нулю этого выражения в точке £ приводит к тому, что
А =
п!
Сравнение двух полученных выражений для А и доказывает теорему.
37
§ 6. Минимизация остаточного члена. Многочлены Чебышева
Займемся задачей минимизации остаточного члена интерполяционной фор
мулы. Далее, если не оговорено иное, под нормой функции понимается
11/11 = max 1/(ж)|.
а<я<о
Из (5.4) следует
II/-М < ^IKII.	(5.4')
Поэтому задача минимизации остаточного члена для широкого класса функ-
ций связана с нахождением такого многочлена
шп(х) = (х - Xi)(x - Х2)...(х - Хп),
норма которого минимальна. Иными словами, требуется так подобрать узлы
интерполяции з?г-, чтобы минимизировать ||о>п||. Эта задача решается с
помощью многочленов Чебышева.
Многочлены Чебышева. При — 1 < х < 1 рассмотрим функцию
Тп(х) = cosnarccosx.	(6.1)
Заметим, что Tb(z) = 1, 7\(х) = х.
Теорема 1. Тп(х) есть многочлен степени п со старшим коэффициен-
том 2П-1, четная функция при четном п и нечетная при нечетном п и
удовлетворяющий рекурентному соотношению
ТП+1И = 2хТп(х) - Т^х).	(6.2)
Доказательство. Докажем рекурентное соотношение (6.2), из которого сле-
дуют остальные свойства.
Пусть а = arccosz. Тогда
Tn+i(x) 4- Tn-i(a;) = cos(n 4- l)cv 4- cos(n — l)cv = 2 cos na cos a = 2xTn(x).
Остальные утверждения легко доказываются по индукции с использованием
реку рентного соотношения (6.2).
В самом деле, докажем, что старший коэффициент Тп(х) равен 2П-1. Стар-
ший коэффициент ГДж) равен 1. Пусть старший коэффициент Тп равен
2П-1 . Тогда старший коэффициент Тп+1 определяется старшим коэффи-
циентом выражения 2хТп(х), который равен 2п.
38
Так же по индукции доказывается свойство четности — нечетности. За-
метим, что Т0(х} = 1, т. е. четная функция, Ti(z) = х — нечетная функция.
Предположим, например, что Тп_1(т) — нечетная функция, а Тп(х) — чет-
ная функция. Из (6.2) получаем
Тп+1( х) — 2жТп( х) Тп-^х) — 2хТп(х) + Тп_\(х^ — Тп^^х).
Определение. Многочлен
Тп(х) = cos n arccos z,
(6.3)
определенный на отрезке [-1,1], называется многочленом Чебышева.
Следует заметить, что в некоторых учебниках определение многочлена
Чебышева дается с использованием рекурентного соотношения (6.2). Имен-
но пусть То(х) = 1, Ti(z) = а?, и при любом п Тп(х) определим из
(6.2). Эти два определения эквивалентны. В одну сторону этот факт мы
установили. При изучении разностных уравнений будет показано (пример
2 § 19), что решением разностного уравнения (6.2) на отрезке [-1,1] является
функция cos п arccos х.
Свойства многочлена Чебышева.
п 7^ m
п = m 7^ 0 .
п = тп = О
Первое свойство доказывается с использованием формулы cos па = -—~^е
второе свойство — непосредственным интегрированием. В качестве упраж-
нения предлагается доказать эти свойства самостоятельно.
Найдем корни многочлена Чебышева. Из уравнения
Tn(z) = cos п arccos х — О
получаем
7Г ,	(2k + 1)7Г ,	, ч
пarccosх = — + кя 2^ = cos-—-—— к = 0,1,..., п — 1.	(6.4)
2	2п
Поскольку Тп(х) = cosпarccosх, то |Tn(^)| < 1-
39
Найдем точки, где |Tn(z)| = 1. Имеем
777» 7Г	z _ ч
п arccos хт = ттг хт = cos- m = 0,(о.о)
71
Подчеркнем, что в (6.4) получено п точек, т. е. найдены все корни Тп(х). То-
чек же максимума модуля — п + 1, что будет использовано в дальнейшем.
Оптимальный выбор узлов. Многочлены, наименее уклоняющиеся от нуля
Определение. Рассмотрим класс многочленов фиксированной степени п и
с фиксированным старшим коэффициентом. Многочлен Qn(x) называется
многочленом, наименее уклоняющимся от нуля на отрезке [а, 5], если на
этом отрезке он имеет наименьший максимум модуля среди всех многочле-
нов данного класса.
Таким образом, многочлен Qn(x) наименее уклоняется от нуля на [а,Ь], ес-
ли
max Юп(ж)| < max |Pn(z)l,
a<x<b	a<x<b
где Pn(x) — всевозможные многочлены той же степени и с тем же старшим
коэффициентом.
Теорема 2. Многочлен Тп(х) = 21~пТп(х) является многочленом, наи-
менее уклоняющимся от нуля на отрезке [—1,1], среди многочленов со
старшим коэффициентом, равным единице.
Доказательство. Предположим противное, что существует многочлен
Рп(х), такой, что
max |Рп(ж)| < max ITn(rr)l.
а<я<6	а<х<Ъ
Рассмотрим многочлен
Qn-i(P) = Тп(#) Рп(ж).
Поскольку старшие коэффициенты у Рп(х) и Тп(х) одинаковы, степень
многочлена Qn_i(x) — п — 1. Рассмотрим значение Qn-i(x) в точках хт,
т. е. в точках, определяемых равенством (6.5). В этих точках Тп(х) имеет
максимум модуля и знак многочлена Qn_i(a?) будет определяться знаком
многочлена Тп(х). Но в этих точках Тп(хк) = ±21-п. Следовательно, мно-
гочлен Qn-ifx) имеет п + 1 перемену знака и по теореме Больцано —
40
Коши имеет не менее п корней. Но это невозможно, если Qn_i(x) не есть
тождественый нуль. Теорема доказана.
Таким образом, для отрезка [—1,1] можно утверждать, что если в каче-
стве узлов интерполяции выбрать корни многочлена Чебышева, определя-
емые формулой (6.4), то норма остаточного члена интерполяционной фор-
мулы будет минимальна и имеет место оценка
|| г || < \\Г\\
||/(ж) Ln\\ < 2п-1пг
Займемся задачей минимизации остаточного члена для произвольного от-
резка [а, Ь]. Рассмотрим на этом промежутке многочлен Рп(х). Замена
b — а
х=^~
переводит отрезок [—1,1] в отрезок
_ 2х
b — а
2
[а, Ь]. С помощью обратной замены
b + а
b — а
получим многочлен
„ . .	„ . 2х Ь + а.
— Рп(т 7	)>
о — а о — а
аргумент которого меняется на [—1,1]. Но у многочлена
- . 2х Ь + а
при такой замене изменится старший коэффициент. Для того чтобы он был
п	(Ь—а)п г>
равен 1, многочлен надо умножить на величину v 2/ . В итоге получаем
многочлен
р М - (д~а)"т (	- b + a't
22п-1 1п'Ь_а Ь-а’’
старший коэффициент которого равен 1 и который наименее уклоняется от
нуля на [а, Ь]. Если в качестве корней многочлена выбрать
b-a (2& + 1)тг b + a
= —— cos---------—+ —— к = 0,1,..., п — 1,
2	2п 2
получим многочлен Тп^а^(х)

(Ь — а)п 2х Ь + а
22п-1 1п'Ь_а Ь_а)'
41
наименее уклоняющийся от нуля на [а, &].
Итогом этих рассуждений является
Теорема 3. На промежутке [а, 6] многочленом, наименее уклоняющимся
от нуля, является многочлен
т - (Ь~а^т < 2х - Ь + а\
Тп^ь){х)- 22n_1 Tn(b_a b_ah
и при оптимальном выборе узлов интерполяции
42
§ 7. Разделенные разности.
Интерполяционный многочлен в форме Ньютона
Вопросы, изучаемые в этом параграфе, необходимы для решения задачи
построения интерполяционного многочлена, который совпадает с функцией
и ее производными до некоторого порядка в заданных точках. Эта задача
решается с помощью интерполяционного многочлена в форме Ньютона.
Определение. Будем называть разделенной разностью нулевого порядка
значение функции f(xi) в точке жг-. Разделенную разность первого порядка
определим равенством
f(T ч _ fM ~ f(xi)
J \Xiy Xk) —
Xk Xj,
Разделенную разность к-го порядка определим с помощью разностей к — 1-
го порядка следующим образом:
Г( .	.	\	fffc+l) ~ /(#!> •••) хк)
J \х1у х2) ••••! хк+1)
хк+1 Х1
Теорема 1. Для разделенных разностей справедливо равенство
ffxiW,
г у v Ж)
i=1 П (xi ~ xj)
j=l
(7-1)
Доказательство будем проводить по индукции. Для разделенной разности
первого порядка имеем:
Х2) = Ж) -	= JkiL +	.
х2 ~ Х1	Х1 ~ х2 х2 ~ Х1
Предположим, что формула верна для разделенной разности к — 1-го по-
рядка. Тогда для разделенной разности к-го порядка имеем:
/(х2;	•••; ^+i) - /foi; ^2; -хк)
•£fc+i xi
В соответствии с индуктивным предположением
№15 ^25 •••  хк+1) =
43
fc+i
П (Xi ~ X'j)
3=2
3&
A fte)
2^ k
1-1 П (Xi ~ Xj)
J=1
j/i
/fofc+l)
k
П C^fc+l — xj)
3=2
№1)	A f(Xj)________________
к	2_^	fc+i
П ("El — xj)	i-2	П (Xi ~ Xj)
j=2	j=2
/(#*)
i=2 П (xi ~ Xj)
3 = 1
3&
Рассмотрим
1 2 П (xi xj) i 2 П {Xi Xj) i 2 П (Xi Xj)
3=2	j=l	j=2
j/i	j/i	j/i
f(xi)
к
П “ хз)
3=1
j/г

и рассмотрим множители при f(xi)
1	1	\	Х], ЗС\ CCf	ЗС^
fc+1	к	'	к+1	fc+1
I~J (Xi — Xj)	IJ (Xi	— Xj)	PJ	(Xi	— Xj)	PI (a?i	—	Xj)
3=2	j=l	j=l	j=l
j^z	j/i	j/i	j/i
Поэтому
f (•£]., Xr2) ••• • *£fc+l)
=	1________f(xk+i) +	1________f(Xj) + A f(xj)
(Xk+l-Xi)^	(Xi-Xk+1) *	i^2 ГТ1/	\
Pl (^fc+l	Xj)	Il (Ж1	Xj) II (.Xi Xj)
j=2	j=2	j=l
j/i
В итоге получаем
44
. _	. „ \	f(xi)
f \xli x2\ — - xk) / y n
t=i jj (a?j Xj)
j=i
З^г
Теорема доказана.
Из формулы (7.1) вытекает ряд следствий.
Следствие 1. Разделенная разность не меняется при перестановке своих
аргументов.
/(жи ...Xi- ...xj- ..хк) = /(^i; ...Xj-, ...хк-, ..хк).
Следствие 2. Разделенная разность есть линейный оператор от функции
/W-
(ai/i + a2/2)(^i; .-хп) =	-хп) +	-хп).
С помощью аппарата разделенных разностей интерполяционный многочлен
Лагранжа можно записать иначе, в форме Ньютона, что, в свою очередь,
позволит решить задачу интерполяции с кратными узлами.
Теорема 2. Интерполяционный многочлен представим в форме Ньютона
Ln(x) =	...\xn')(x-xi)(x-xi)..\x-xn-i').
(7-2)
Доказательство. Заметим, что
f(x) - Ln(x) = f(x; хг,...; хп)шп(х).	(7.3)
В самом деле,
П (х “ хз)
п	j — i
fix') - Ln(x) = /(г-) - £	----------=
г=1	jj ((Ej — Xj)
3=1
з/i
45
С учетом того, что шп{х) = (ж — х\}(х — х{)...(х — хп), равенство (7.3)
доказано.
С другой стороны,
Ln(x) = (Ln(a?) - Ln-iW) + (Ln-iM) 4-... + (Т2(*) - ТДж)) + М(ж) (7.4)
Но разность Lk(jc) — Lk-i(x) есть многочлен степени А; —1, обращающийся
в нуль в точках х^ i = 1, 2,..., k — 1. Поэтому
-^fc(^)	^к—1 (*^) = Ак—1^к— 1(*е)>	(7.5)
где' Шк-1(х) = (х- Я1)(я - х2)...(х - Xk-1).
Рассмотрим Lk(x) — Lk_i(x) в точке хк.
f(%k)	-^fc—l(^fc) Ак_]Шк_ 1 (х&).	(7.5 )
Отсюда
л ffak)	A'fc—i(^fc)
A-k-l = -------7—x------•
Wk-l(Xk)
В силу (7.3)
Ak-i = f(xk-,xx;...-,xk-iY	(7.6)
Поэтому (7.4) с учетом (7.5) и (7.6) и того, что I/i(x) = /(a?i), можно запи-
сать следующим образом
Ln(x) = f(x^+f(xi; ж2)^1 W+/(a?i; х2; x3)u>2(x)+...+f(xr, х2\...; xn)wn_i(x).
Теорема доказана.
Следствие. Для разделенной разности f(x\\x2 : ...;#&) справедливо пред-
ставление
ft . ™	\ _ f(Xk) ~ Lk-l(Xk)	/„ „X
/(Ж1,Ж2 . ...уХк)	.	.	0’7)
Мк-\\хк)
Практическое построение интерполяционного многочлена в форме Ньюто-
на можно проводить с помощью формулы (7.1), вычислив соответствующие
разделенные разности. Существует схема Эйткина для построения разде-
ленных разностей и интерполяционного многочлена в форме Ньютона, из-
ложенная, например, в [1].
46
§ 8. Интерполирование с кратными узлами
В этом параграфе будет рассмотрена задача построения интерполяцион-
ного многочлена, совпадающего с заданной функцией и ее производными до
некоторого порядка в точках Xi.
Итак, требуется построить многочлен Р3(х) степени s — 1, удовлетво-
ряющий условиям
Л(®1) = Щ1),
Л(®2) =
р'(^) = f(X1),
Л'Ы = Г(Х2),

р.ы=/ы, р;(хт) = гы, . .	=
Здесь s = Zi + l2 -I-... + 1т.
Рассмотрим пример. Пусть требуется построить многочлен, который сов-
падает с заданной функцией в точке #1, а в точке х2 совпадает с функцией
и ее производными до второго порядка.
РзЫ = /(Х1),	=	Р!Ы = Г{гг).
В ^-окрестность точки х2 добавим еще 2 точки х% и х% и построим по
точкам ЯГ1, х2, ^2, х%интерполяционный многочлен Рз(х)
= f(xi) + /(xi;x2)(x - Ж1)+
+№15 х2; х%)(х - xi)(x - х2) +	х2; ж2; х%)(х - х^(х - х2)(х - х|).
Устремим теперь е к нулю. При этом х% —> х2, х% —> х2.
^з(ж) = /(^1) 4- /(®i;x2)(x - Xi)+
4-/(a?i; ж2; x2)(z - xi)(x - х2) + /(хи х2; я2; х2)(ж - хх)(х - х2)2.
Выясним смысл выражений
/(xi;x2;x2),	/(^х;^2;ж2;х2).
По определению
г/ . ч _ /(^2! х2} - f(X1; х2)
Д/2	1
47
Но, сравнивая (7.7) и (5.4), замечаем, что
/(zijtfi) =
Более того, если точка перечисляется т раз, то
,,	у	/ЯП
m раз '	'
Поэтому
J (^1> ^2) ^2)
Х2 ~ Xi
Что касается f(xi; #2; <г2; ж2), то
/(xi; х2] х2; х2) = —-—(f(x2; х2‘, х2) - /(®г, х2; я2)).
Х2 — XI
Первое слагаемое определено формулой (8.1), a /'(xi; Х2; я2) только что
определено.
В общем случае процесс построения многочлена Ps(x) степени s — 1 за-
ключается в следующем. В с-окрестность кратной точки Xj помещаются
lj-i вспомогательных точек xj, xj,	xl-~l. При г —> 0 получим много-
член вида
Ps(z) = Ajf(xi) + А}(а: - a?i) + . . . + А^1-1^ - xi)Z1-14-
—zi^+A.^z ——z2)+. . . + А2-1(я —xi)Z1-1(z —• • •
+^m U (x - xi)li + Am П (X ~ Xi)li(X ~ X™) + Am И (Ж “ х^\х ~ xm)2+
i=l	г=1	i=l
• • • + Arn ntx-Xi^x-x^1™-1.	(8.2)
1=1
Здесь A? = 1. Далее Af = f(xy, #1;...; x^) и, согласно (8.1).
k+1 раз
A"= k = 1J1 - L
kl
Коэффициент A® = /(^i; a?i;...; a?i; £2)- Покажем, что это выражение имеет
h раз
смысл. Доказательство будем проводить по индукции. Заметим, что /(a?i; х2)
48
имеет смысл. Это обычная разделенная разность 1-го порядка. Пусть нами
определена разность /(a?i; a?i;...; xij х2) • Но
к раз
f(xi;xi] ...;zi;z2) = —-—(№i;zr, ;xr,x2) - /(ад;
------V	'	X2 — Xi	4-~"v
k+1 раз-------------------------k раз	k раз
Первое слагаемое в скобках определено по индуктивному предположению.
Второе слагаемое определено формулой (8.1). Таким обазом, коэффициент
А® определен. Продолжая аналогичные рассуждения для других коэффи-
циентов формулы (8.2), приходим к заключению, что все коэффициенты в
(8.2) имеют смысл.
Коротко остановимся на оценке остаточного члена интерполяционного
многочлена в форме Ньютона. Если построить по точкам Xi и точкам х[
интерполяционный многочлен, то, согласно (5.4), получим
f{x) - Ра(х) = -сиДж).
S!
При е —> 0 многочлен u)s(x) имеет вид
— _|_ J_(*^ Xi)
i=l
и оценка (5.4) в этом случае имеет вид
=	(8.3)
В заключение этой темы заметим, что задача интерполяции функции
многих переменных значительно сложнее и в настоящее время еще недоста-
точно разработана. Для иллюстрации трудностей при интерполяции функ-
ции многих переменных рассмотрим 2 простых примера.
Пусть требуется построить интерполяционный многочлен, совпадающий
с функцией двух переменных f(x,y) в трех точках. Ясно, что многочлен
1-го порядка — уравнение плоскости, проходящей через эти три точки, — и
будет этим интерполяционным многочленом. Но если эти три точки лежат
на одной прямой, то этот многочлен не единственнен, их бесконечно много.
Пусть теперь требуется построить интерполяционный многочлен, совпа-
дающий с функцией f(x,y) в четырех точках
49
^1(^1,У1), М?(Х1,У2)> Мз(х2,Уг),	Ясно, что это должен быть
многочлен двух переменных, степени выше первой, поскольку для постро-
ения многочлена 1-й степени Pi (ж, у) = ао + № + Ь\у надо определить три
коэффициента aQ,ai,bi. Но функция должна совпадать с многочленом в
четырех точках, т. е. получаем систему из четырех уравнений для опреде-
ления 3-х неизвестных, которая в общем случае неразрешима. Что касается
многочлена 2-й степени Р2(х,у) = ао + ааХ + Ъ±у + Яс-^ху 4- а^х2 + Ь2У2 ,
то для определения шести коэффициентов получаем систему из четырех
уравнений, и ясно, что без дополнительных условий решение этой системы,
вообще говоря, не единственно. Таким образом, уже на начальном этапе изу-
чения мы сталкиваемся с трудностями, связанными с теоремами существова-
ния и единственности интерполяционного многочлена для функции многих
переменных. Изучению проблематики теории приближений посвящен ряд
работ отечественных и зарубежных математиков. Здесь следует назвать, в
первую очередь, работы П.Л. Чебышева, А.Н. Колмогорова, П.С. Урысона.
С отдельными результатами по этой теме можно ознакомиться, например, в
II].
50
§ 9. Ортогональные многочлены и их свойства
Рассмотрим пространство функций таких, что
ь
У p(x)f2(x)dx
а
ОО,
где р(х) > 0. Зададим в этом пространстве скалярное произведение по
формуле
ь
(f,9) = У p(x)f(x)g(x)dx.
а
Как обычно, если в гильбертовом пространстве для пары элементов имеет
место (/,g) = 0, то говорят, что эти элементы ортогональны.
Рассмотрим пространство многочленов степени пив этом простран-
стве рассмотрим линейно независимую систему элементарных одночленов
1	/т*
±,	da ,	,	,	, da .
Теорема 1. Существует линейно независимая система многочленов со стар-
шим коэффициентом, равным единице,
Рк(х) =Рк,к + Рк,к-1% + ••• + +Pk,ixk~1 + хк к = 0,1,..,п,	(9.1)
таких, что
(Pj,Pj) = 0 при г j.	(9.2)
Здесь индекс к введен для того, чтобы идентифицировать принадлежность
коэффициента pi к многочлену заданной степени.
Доказательство. Пусть Ро(т) = 1. Предположим, что система таких
многочленов построена при некотором к и пусть
к
Рк+1(х) = irfe+1 - 57 akiPi(x).	(9.3)
г=1
Коэффициенты выберем из условия (Pk+i,Pi) = 0 i = 1,2
Отсюда
^ki —
(9.4)
Следовательно, построенный многочлен будет ортогонален всем многочле-
нам РДт) i=l,2,...,k. Далее, представление многочлена (9.1) означает, что
51
многочлен Рк(х) может быть представлен в виде скалярного произведе-
ния двух векторов (рк,к,Рк,к-1,-,Рк,1, 1) и (1,ж,a?fe), и переход от
системы элементарных одночленов к системе ортогональных многочленов
осуществляется преобразованием
/ 1 \
х
„2
X
\ Рк,к Рк,к-1 Рк,к—3 ...	1 / \Хк J
с матрицей В
(9-5)
\ Рк,к Рк,к-1 Рк,к-3 •••	1 /
Таким образом, переход от линейно независимой системы элементарных од-
ночленов к системе ортогональных многочленов осуществляется невырож-
денным преобразованием, поэтому и система ортогональных многочленов
линейно независима.
Теорема доказана.
Определение. Систему многочленов, удовлетворяющих условию (9.1),
называют ортогональной системой.
Определение. Матрица В называется матрицей ортогонализации. Об-
ратный переход осуществляется также невырожденной треугольной матри-
цей А
А =	/ 1 «1,1 «2,2	0 1 «2,1	0 0 1	...	0	\ ...	0 ...	0
	\ ^к,к	«fc,A:—1	«fc,fc-3	... 1 /
Следствие 1. Ортогональный многочлен Рк(х) ортогонален любым мно-
гочленам меньшей степени. Доказательство легко провести по индукции, что
предлагается проделать самостоятельно в качестве легкого упражнения.
Следствие 2. Система ортогональных многочленов образует базис в про-
странстве многочленов.
Следующая теорема является основной в теории ортогональных функций.
52
Теорема 2. Для ортогональных многочленов справедливо равенство
Pk+i(z) + dk+i,kPk(x) + ^+i,fc-iPfc-i(x) - жРДт) = 0,	(9.6)
причем
dk+i,k-i > О
и
ь
f xp(x)P£(x)dx
dk+i,k = ~b------------,	(9.7)
/р(ж)Р^2(х)с?ж
a
fp(x)P%(x)dx
dk+i,k-i = у--------------•	(9.8)
f P(x)Pk-i(x)dx
a
Доказательство. Пусть построена система ортогональных многочленов
Po(z), Pi (ж),.... , Pk+i(x). Рассмотрим многочлен Qk+i(x) = хРк(х), где
РДж) — ортогональный многочлен степени к. Qk+i(x) уже не относится
к системе ортогональных многочленов. Разложим Qk+i(x) = хРк(х) по
системе ортогональных многочленов.
^Pfc(^) — dk+lyk+lPk+l Н- dk+ltkPk(*е) “I” dk+ltk—1Рк— 1 (^) Н- • • • *Р Фг+1,О-^*Ь- (9.9)
Поскольку у хРк(х) и Рк+1(х) старшие коэффициенты одинаковы и
равны единице, dk+i,k+i = 1-
Коэффициенты </fc+i)Tn этого разложения определяются по формулам
ь
f xp(x)Pk(x)Pm(x)dx
—ь-----------------	(9.10)
fp(z)P£,(x)dx
а
С другой стороны, согласно следствию 1, ортогональные многочлены орто-
гональны всем многочленам меньшей степени. Поэтому
(яРь Pm) = (Рк, хРт) = 0 при к > т + 1.
Согласно следствию 1, dk+i,m — 0 при к > т + 1 ив разложении (9.9)
равны нулю коэффициенты	Фь+1Л-з, • • • , Фк+1,о- Разложение
53
(9.9) при этом имеет вид (9.6). Докажем теперь, что dk+i,k-i > 0. Согласно
(9.10)
л _ (xPk,Pk-i)
dk+i,k-i- {Pk_bpk_lY
Знаменатель этой дроби > 0. Покажем, что и числитель > 0. Разлагая
xPk-i(x) по системе ортогональных многочленов
•ЕРк—1(®) = -Pfc(^) Н” frfc—1-Pfc—1 (*^) + • • • bOPo,
получим
(хРк,Рк-1) = (Рк,хРк-1) = (Рк,Рк) > 0.
Теорема доказана.
Теорема 3. Все корни ортогонального многочлена Рк(х) различны и
принадлежат открытому интервалу (а, 6).
Доказательство. Предположим противное, на (а, &) находятся т корней
а?1,Т2, . . . ,хт многочлена Рк(х) и т < к. Предположим сначала, что
все корни различны, и рассмотрим многочлен
т
Qmfa') = J~[(*^ xi)‘
i=l
Поскольку р(х) > 0, подынтегральное выражение в
ь
У p(x)Pk(x)Qm(x)dx
а
не меняет знак на (а, &) и поэтому значение интеграла > 0. С другой сто-
роны, ортогональный многочлен ортогонален всем многочленам меньшей
степени (следствие 1) и поэтому значение интеграла = 0. Получено проти-
воречие. В случае кратного корня xq множитель (х — яо) в многочлене
Q(x) можно исключить, если корень четной кратности, и оставить его в
первой степени, если он нечетной кратности.
Теорема доказана.
Теорема 4. Корни многочленов Рк+\{х) и Р^(т) разделяют друг друга.
Иначе, пусть ж* < я* < • • • < хк ~ корни многочлена Рк (т) и
а < х*+1 < х%+1 < ... < х%+1 <	< b — корни многочлена Pfc+i(x).
Тогда
а < a?i+1 < х± < а?2+1 < х% < . . . < х^ < х^^ < Ъ. (9.11)
54
Доказательство. Докажем это утверждение по индукции. Докажем сначала,
что для корней многочленов P2(z) и Pi (ж) имеет место цепочка неравенств
а < Xi < х\ < а?2 < Ь. Рассмотрим рекурентное соотношение (9.6) для
многочленов Pq(z) = 1, Pi(x) и РДх) в точке, где РДж^) = 0.
P2(z}) +d2)oPo(z}) = О-
Поскольку с?2,о > 0, то РДж}) < 0. Но Р2(ж) ~ многочлен второй степени
со старшим коэффициентом 1, и вне корней он принимает положительные
значения. Поэтому должны существовать точка х% : а < х% < х\ такая,
что P2(a?j) = 0, и точка х% :	< х% < Ь, что РДя^) = 0- Следовательно,
для многочленов РДт) и Р2(х) утверждение теоремы имеет место.
Предположим теперь, что утверждение теоремы имеет место для корней
многочленов Pk-i(x) и Pk(x).
п 'T'k 1 <•'	<<	1 <"*	<" т^~ 1	h
Cl \	Д/j \ Д<2	**/2	Д/д._j \ Ц/Д. и.
Докажем, что оно имеет место для корней многочленов РДт) и Р^+Дх).
Заметим,что многочлены Р^-Дя) и Pk+i(x) имеют одинаковый знак
вне корней, так как старший член у этих многочленов равен соответственно
аД-1 и rrfc+1 , т. е. х одновременно или в четной, или в нечетной степени.
Рассмотрим (9.6) в точках РД;г|) = 0 j = 1, 2, . . . , к. В этих точках
значения Pfc+i(a?j) и Р^_Джу) имеют разные знаки. Согласно индуктивно-
му предположению, между двумя последовательными корнями РДх) есть
один корень Pfc-i(r) , и, следовательно, в точках Xj и х*+1 значения
Pjt-i(x) имеют разные знаки. Но тогда и Pk+i(x) будет иметь разные знаки
в точках х1- и 3j+1. Тем самым мы заметили, что между двумя последо-
вательными корнями многочлена Pk(x) есть корень многочлена Pjt+Дх).
Итак, на к — 1 интервале (ж|,ж*+1) j = 1,2, . . . , к есть, по крайней
мере, один корень Р^+1(а;). Рассмотрим теперь интервал [х%, Ъ]. В точке
х% значения Рк-i и Р^+Дж) имеют разные знаки, поэтому Р^+Дх) на
(z'£,6) должен сменить знак, чтобы при х > b у многочленов Р^_Дх) и
Pfc+i(z) были одинаковые знаки. Значит, на этом интервале также должен
быть корень Р^+Дя). Такая же ситуация на интервале (а, гт*). Итак, най-
ден к+1 интервал, где есть корни Р^+Дя). Одновременно замечаем, что
там может находиться только один корень Pfc+Дя) и корни многочленов
Рк(х) и Pfc+Да;) удовлетворяют условию (9.11). Теорема доказана.
55
Теорема 5. Пусть весовая функция р(х) является четной функцией от-
носительно середины отрезка. Тогда многочлены Рк(х) являются четны-
ми функциями относительно середины отрезка при четном к и нечетными
функциями относительно середины отрезка при нечетном к.
Доказательство. Замена
_ 2х — (Ь + а)
Ъ — а
переводит отрезок [а,Ь] в отрезок [-1,1], поэтому достаточно доказать утвер-
ждение теоремы для отрезка [-1.1] . Имеем Ро(х) = 1 и в (9.4)
1
f p(x)xdx
aio = ——5------= о,
поэтому РДж) = х — нечетная функция. Предположим теперь,что Pk-iix)-
нечетная функция, а Рь(ж)- четная функция. Тогда в (9.6) dk+i,k = 0, так
как (9.7) в числителе стоит интеграл от нечетной функции на симметричном
относительно нуля интервале. Поэтому
Pfc+1(-^) — ^'Pfc(*^)	Фг+1,&—1-Pfc—1 (*^)
Pfc+1( ^) =	З') ^fc+l,fc-1-Pfc—1 ( *^) “
Т	1-Pfc—1 (з?)	Pfc+1(^')'
Другой случай рассматривается аналогично.
Ортогональные многочлены являются важным инструментом в теории
интерполяции и смежных задачах. Изучая вопросы, связанные с оптималь-
ным выбором узлов интерполяции, мы ввели в рассмотрение многочлены
Чебышева. Нетрудно заметить, что многочлен Тт(а?) является ортогональ-
ным многочленом на [-1,1] с весовой функцией
₽w = Tib-
При построении квадратурных формул мы убедимся, что наиболее точные
формулы для приближенного вычисления интеграла получаются при ис-
пользовании в качестве узлов квадратур нулей ортогональных многочленов.
Приведем наиболее употребительные в вычислительной практике ортого-
нальные многочлены.
56
Многочлены Якоби. Эти многочлены получаются при построении систе-
мы ортогональных многочленов на [-1,1] с весовой функцией
р(х) = (1 — #)Q(1 + х)&	а, (3 > — 1.
Многочлены Якоби имеют вид
(_ 1\П	лп
2пп!	ахп
Частным случаем при р(х)=1 являются многочлены Лежандра, с которыми
мы встретимся при построении формул приближенного вычисления инте-
гралов.
57
§ 10. Численное дифференцирование
Задача приближенного вычисления производных от функции, заданной
таблично, обычно решается вычислением производных от интерполяционно-
го многочлена. Если известны значения функции в точках Xi, Xi. . . ,хп,
то
к L^(x).	(10.1)
Оценим погрешность формулы (10.1). Из (7.3) имеем равенство
f{x) - Ln(x) = f(xi; х2;...; хп; x)wn(x).
Отсюда
к
f»\x) - L»\x) = £ Cl(f(zv,   ;	х))^-‘>(х).
г=0
Выясним смысл выражения
(/(xi;...;xn;x))w.
В § 8 получена формула
f{xi,Xj^...,Xi) - /	_	’
m раз v '
из которой следует, что
(/(я?1;...; хп; ж))(г) = г!/(хг,...; хпх;...-,х/).
i+1 раз
Поэтому
к
f<k\x) - LW(x) = £ СЦ(х1;...;	Нк-'\х).
1=0	i+1 раз
В практических вычислениях обычно используются более простые сооб-
ражения. Так, если в точках х^ и х2 заданы f(xi) и f(x2) , то
L2(x) = Дх,)^-+ f(x2)^-
Х~[ Х2	Х2 Х~у
58
и
f (х) «
X2 — Xi
При этом
£
Если h = |o?2 — 2^1 h TO
\\f - l'2\\ = ишь.
Более подробно с разностными аналогами производных и соответствующи-
ми погрешностями мы познакомимся при изучении сеточных уравнений и
разностных методов решения краевых задач.
59
§ 11. Численное интегрирование.
Квадратурные формулы Ньютона - Котеса
В этом разделе мы получим формулы для приближенного вычисления
интеграла
ъ
1(f) = Уp(x)f(x)dx.
а
Функция р(х) называется весовой, и предполагается, что р(х) > 0. Квад-
ратурные формулы бывают простые и составные. Методика получения про-
стых квадратурных формул состоит в том, что функция f (х) на заданном
отрезке под интегралом заменяется на ее интерполяционный многочлен, что
позволяет получить достаточно простые формулы приближенного вычисле-
ния интеграла и оценить погрешность. Составные формулы получаются,
когда весь отрезок интегрирования разбивается на отрезки длины h, и на
каждом отрезке используются простые квадратурные формулы. Займемся
выводом простых квадратурных формул. Все построения мы будем прово-
дить на отрезке [0, h]. Поскольку в некоторых случаях мы будем исследо-
вать дополнительные свойства квадратурных формул в случае симметрии
функции р(х) и симметричного расположения узлов, удобно ввести отоб-
ражение отрезка [—1,1] на [0, h].
z = £(t + l).	(11.1)
£
Изложим общую схему построения этих формул. Рассматриваемые ниже
формулы приближенного вычисления интеграла называются квадратурны-
ми формулами Ньютона — Котеса. Все квадратурные формулы имеют вид
/ f(x)dx « y^Cjf(xj).
о	2=1
Зададим некоторые узлы Xi = ^(ti + 1) i — 1,2,... ,п и построим интер-
поляционный многочлен
П (ж - Xj)
П	j = l
ш =	------
г=1	jj (Xi — Xj)
J=1
J#
60
Будем считать, что
h	h
1(f) = I f(x)dx « I Ln(x)dx = Sn(f).	(11.2)
о	0
Здесь
n
sn(f) = £сшо,
i=l
а коэффициенты Cj вычисляются по формуле
c. = / p(z)~--------------dx	(11.3)
О II (xi ~ Xj)
J=1
j/i
или, сделав замену (П-1),
.	п (* - А?)
1 A	3=1	1
Ci = z / Р° W2?-----------dt = 3-Dit
3=1
3&
где
.	П (^ - h)
I	3=1
Di = / P°(t)^-------------dt.	(11.3')
-1	II (ti ~ tj)
3=1
3^i
Оценим погрешность построенной квадратурной формулы. В силу (5.4)
1Яп(/)1 max|/WJW f |р(а.ЯК(а.)|^. (п 4)
7*i.	/
0
После замены (11.1)
61
I*W)I < hn+'™^fJ— [ |p°WIH«|dt. (11-4')
»v!	</
-1
Следующая теорема иллюстрирует свойства квадратурных формул в случае
симметричного расположения узлов.
Теорема 1. Пусть весовая функция р(ж) -- четная относительно середины
отрезка [0,7г], а узлы расположены симметрично относительно середины
отрезка. Тогда коэффициенты Ci, соответствующие симметричным узлам,
равны.
Доказательство достаточно проделать для отрезка [—1.1]. Пусть узлы
на [—1,1] расположены симметрично: ti == £n-z+i- Рассмотрим
.	П(*	fy)
J.	з=1
Di = /	-------dt
-1	П (^1 ~ fy)
3=1
3/г
шп(1) j.
где cun(t) = (t — ti)(t —	— В силу симметричного расположения
узлов функция wn(t) есть четная функция при четном п и нечетная
функция при нечетном п. Далее, производная от четной функции есть
нечетная функция, а от нечетной — четная. Поэтому w'n(ti) = —iv'n(tn-i+i)
при четном п и w'n(ti) = a^(tn_i+i) при нечетном п. Сделаем замену
переменной интегрирования t = —т и, заменив ti на —получим
1 1
[ тРа\___Шп^'>__dt = [ v°(-t)_______________________dr
Iр 1 >(t-Jp( т>(-г+«„_,+1H(-t„_i+1)dT-
1 -1
В обоих случаях получаем
Теорема доказана.
Будем говорить, что квадратурная формула точна для некоторого клас-
са функций, если с ее помощью вычисляется точное значение интеграла.
62
Квадратурные формулыНьютона — Котеса точны для многочленов степени
п — 1. В самом деле, интерполяционный многочлен, построенный по п уз-
лам для функции, которая сама есть многочлен степени т < п, совпадает
с этом многочленом в силу теоремы единственности для интерполяционного
многочлена.
Симметричное расположение узлов позволяет не только уменьшить рабо-
ту по вычислению коэффициентов квадратурных формул, но в ряде случаев
увеличить точность квадратурной формулы. В самом деле, пусть р(х) — чет-
ная функция и число узлов п нечетно. Тогда такая квадратурная формула
точна уже для многочленов степени п. Для доказательства этого утвержде-
ния достаточно заметить, что для элементарного одночлена хп
1
у* p(x)xndx = О
-1
з силу нечетности подынтегральной функции. С другой стороны,
п
Е с.х" = О
г=1
з силу нечетности п и симметричного расположения узлов.
Займемся выводом наиболее употребительных формул при р(х)=1. Сра-
зу заметим, что константа является четной функцией относительно любого
интервала.
Формула прямоугольников. Различают формулу прямоугольников с ле-
восторонним, правосторонним и центральным узлом. Выведем формулу пря-
моугольников с центральным узлом. В этом случае ti = 0, х-^ = | и
1
G = | У* Idt = h
-1
и получаем
h
[
J	"
1
Оценка погрешности согласно (11.4)
Ri = maxf'(x)—.
[О,Л] k J 4
63
Интересно отметить, что здесь мы сталкиваемся с описанной выше ситуаци-
ей: число узлов нечетно и весовая функция четна. Это обстоятельство поз-
воляет предположить, что на самом деле оценка погрешности может быть
лучше. Действительно, построим квадратурную формулу прямоугольников
с кратным центральным узлом. Интерполяционный многочлен в этом слу-
чае имеет вид
Z	Z	Z
Вычисляя коэффициенты С\ и . замечаем, что С\ = h, а С% = 0.
Таким образом, получена та же квадратурная формула
h
J	Z
1
но с оценкой погрешности
1
h3 шах |///(зг)|
24
Rt =
16 J
-1
Составная формула прямоугольников получается, когда весь отрезок ин-
тегрирования [а,Ь] разбивается точками Xj. на отрезки длиной h. с последу-
ющим применением полученной формулы на каждом интервале. Составная
формула прямоугольников имеет вид:
п—1	,
г=0
Здесь
b — а
h =------, Xi = а + ih, i = 0,1, . . .п — 1.
п
Поскольку nh = Ь — а, погрешность составной формулы прямоугольников
24
Формула трапеций. Аппроксимируя функцию f(x) многочленом L^tx)
с узлами Xi = 0 и Х2 — h (ti = — 1 и fg = 1), получаем
1
С =h
1	2 J 2
-1
Z
64
В результате приходим к формуле
h
[ f(x)dx « £(/(0) + /(Л))
о
с оценкой погрешности
Д3 max|/"(a?)| Г 2 _ 11 Л _ h3 тах|/"(ж)|
16 J 1	1	12
-1
Составная формула трапеций имеет вид
5П(У) -	4- f(xj) Xi = a + ih h = -—-	(11.6)
2	t—'	n
z=l
с оценкой погрешности
h2 max | f" (z) |
R~ 12
Формула Симпсона. Так же, как и в случае формулы прямоугольников
с центральным узлом, можно заметить, что интерполяционный многочлен
с тремя узлами Xi = 0, Х2 = а?з = h ( ti = —1, <2 = 0, <з = 1 )
и интерполяционный многочлен с четырьмя узлами = 0, Х2 — х3 =
X4 = h (<i = —1, <2 — <3 = 0, <4 = 1) дает одну и ту же квадратурную
формулу. В самом деле,
L4(z) = L3(x) + /(0; h) х (х - -) (х - h),
где L3(x) — интерполяционный многочлен, построенный по трем точкам
Т1 = 0, Х2 = х3 = h. Но
поскольку под интегралом нечетная относительно середины отрезка функ-
ция. Поэтому квадратурная формула Симпсона с тремя узлами такая же,
как и формула с четырьмя узлами, когда средний узел двукратный.
65
Вычислим коэффициенты квадратурной формулы Симпсона, соответ-
ствующие узлам ti = — 1, <2 = 0) h = 1 :
r r h f ^(< + l) ,, h
C3 — Ci = — /    at = — .
2 J 2	6
-i
1
x-v h f, 9ч ,	4/z
c2 = - / (1 - t2)dt = —.
& J	о
-1
Получена квадратурная формула
h
[ f(x)dx^[f(0) + 4f& + f(h)].
J	о	2
о
Оценим погрешность квадратурнойформулыСимпсона. Согласно (11.4')
/г5тах|/(4)(я)| Г 2 2	П|,. max|/(4)(х)|
Ri - -----2^5!-----Pt{t ~ 1)|d< =--------2880-----'
-1
Составная формула Симпсона имеет вид
SN(f) =
=	+f(b)+2f(Xi) + 4f(x2) + 2f(x3) + 4/(t4) + ... + 2f(xN^)] (11.7)
Здесь N четно, Xi = a + ih, h = Погрешность составной формулы
трапеций
h4 max|/(4)(z)|
2880
66
§ 12. Квадратурные формулы Гаусса
Эффект увеличения точности квадратурных формул при специальном
выборе узлов позволяет предположить, что возможно построение квадра-
турных формул, точных для многочленов степени 27V — 1. Это предположе-
ние обосновывается следующими соображениями. В этом случае квадратур-
ная формула должна давать точное значение при вычислении интегралов
от элементарных одночленов 1, х, х2,... ,x2N-1. Следовательно,
ь
p{x)dx = C*i + С*2 Н- • • • + Cn
а
Ь
у* р(х) xdx = Cixi + С2Х2 + . • • + Cn^n
а
Ъ
у* р(х) x2dx = CiX} + С2Х2 + ... + C^x2N
а
b
у* р(х) X2N~rdx = C\X2N~X + C2X2N-1 + • • . + C^X™-1 .
a
Получена нелинейная система из 2N уравнений относительно 2N неиз-
вестных Ci, С2, • • •, Cni #i, ^2, • • •, %n- Если эта система имеет решение и
узлы Xi принадлежат отрезку [а,Ь], то такая квадратурная формула может
быть построена.
Предположим, что такая квадратурная формула существует.
Теорема 1. Пусть квадратурная формула
SN(f) = Cifixt) + C2f(x2) + ... + CNS(xN)
точна для многочленов до N-1 степени и пусть
= (х — Х1)(т — Х2) .. .(х — Xn)
67
— многочлен, корнями которого являются узлы этой квадратурной форму-
лы. Тогда для любого многочлена Р^-^х) степени не выше N-1 имеет место
ь
р(х}шм(х}Рм-\(х}(1х = 0.
а
Доказательство. Рассмотрим многочлен
Qzn-i(x) = wk(x)Pn-i(x),
где Pn-i(x) — многочлен степени не выше N-1. По предположению теоремы
квадратурная формула точна для многочленов до 2N-1 степени. Следова-
тельно,
ь
У p{x)ujN{x}PN_i(x)dx = Ciujn(x\)PN-i(xi) + ... + CNuN(xN)PN-\(xN).
а
Но (jj{xj) — 0 при j = 1, 2,..., N, что и доказывает теорему.
Эта теорема показывает, что узлы квадратурной формулы, точной для
многочленов до 2N — 1 степени, являются корнями ортогонального мно-
гочлена, которые рассмотрены в § 9.
Теорема 2. Пусть
SN(f) = ОДи) + C2f(x2) + ... + CNf(xN)
— квадратурная формула Ньютона — Котеса, узлами которой являются кор-
ни ортогонального многочлена Ру (ж). Тогда эта квадратурная формула
точна для многочленов степени 2N — 1.
Доказательство. Квадратурная формула Ньютона — Котеса точна для
многочленов до N-1 степени. Пусть Xi i = 1,2,..., N — корни ортогональ-
ного многочлена, которые являются узлами квадратурной формулы, и пусть
шу = (х — Х1)(х — а?2) • • • (я — жу).
Поскольку старший коэффициент у многочлена равен единице, то мно-
гочлен шу и есть ортогональный многочлен Ру (ж). Разделим произволь-
ный многочлен Q2N-i(x) на шу(я)
Шу (ж)
68
Здесь Tn-i(x) и Rn-iIx) — многочлены степени не выше N-1. Поэтому
Q2n-i(x) можно представить в виде
ф2лг-1(я) = wn(x)Tn--[(x) + Rn-[(x).
Покажем, что квадратурная формула точна для многочлена Q2W-i(z)- В
силу того, что и>дг(ж) = Pn(%)
ь
p(x)uj^{x)T^^\(x)dx = О
а
И
Ь	Ь
У p(x)Q2N-i(x)(x)dx = У p(x)Rn-i(x)(x).
а	а
Но квадратурная формула Sx(f) точна для многочленов степени N-1, и
поэтому она точна для многочленов степени до 2N-1. Теорема доказана.^
Таким образом, квадратурные формулы Гаусса — это квадратуры Ньюто-
на — Котеса, узлами которых являются корни ортогонального многочлена.
Из других свойств квадратур Гаусса отметим положительность коэффи-
циентов Ci. В самом деле, рассмотрим многочлен
N
T2N_2(x) =	- Xi)2.
г=1
Этот многочлен степени 2N-2 равен нулю в точках х^ i к. Квадратурная
формула Гаусса точна для таких многочленов. Поэтому
ь
У p(x)T2N_2(x)dx = CkTN_2(xk)
а
и положительность Ск при Ь>а очевидна.
Построим в качестве примера квадратурную формулу Гаусса с тремя уз-
лами для вычисления интеграла
1
У f(x)dx «	+ c2f(x2) + с3/(ж3).
-1
69
Эта квадратура должна вычислять точное значение интеграла в случае,
если
f(x) = 1, /(ж) = т, f(x) = X2, f(x) = х3, f(x) = т4, f(x) = X5.
Подставляя хт, т = 0,1, 2, 3,4,5 в квадратурную формулу, получаем
Cl + С2 -Ь с3 = 2
С1Т1 4- С2Х2 + С3Т3 = О
С1Т? + c2xl 4- С3Т3 = |
С1Т? + с2х% 4- С3Т3 = О
4	2	4
Cl^i + С‘2^4 + с3#3 = 7
5
cizf 4- с2Х2 4- С3Т3 = О
Поскольку р(х)=1 — четная функция, узлы расположены симметрично (тео-
рема 5 § 9), а коэффициенты, соответствующие симметричным узлам, равны
(теорема 1 § 11). Второе, четвертое и шестое уравнения — следствие этого
свойства. По этому же свойству
Х2 = 0 и Xi = —Х2 и получаем
2	2
2ci 4- С2 = 2, 2ciZi = -
2
5
Решив эту простую систему, получим
Х2 = 0 Тз -
Х1 =
5	8	5
С1 = 9 с2 = - с3 = -.
2
Так как отрезок [-1,1] отображением
h
х = —
2
переводится в отрезок [0,h], то для интервала интегрирования [0,h] в квад-
ратурной формуле
~ dif(xi) + d2f(x2) + ^3/(^3)
для узлов и коэффициентов этой квадратуры получаем
h h [з
X1 ~ 2 ~ 2 V 5’
h h /З
I3=2 + 2V5
h
х2 =
70
. 5h 1 4h . 5h
dl=18 d2 = T d3=lS-
Составная формула Гаусса с тремя узлами, которая является аналогом
формулы Симпсона, имеет вид
п-1
/ f(x)dx «	+ a?i) + d2f(yi + x2) + d3f(yi + x3)].	(11.8)
a	2=0
Здесь
n
С оценками погрешности квадратурных формул Гаусса можно ознако-
миться, например, в [2].
71
§ 13. Наилучшее приближение в гильбертовом пространстве
Рассмотрим следующую задачу. Пусть в гильбертовом пространстве Н за-
дана линейно независимая система gi,gi,... ,дп- Требуется заданный эле-
мент f приблизить линейной комбинацией элементов #1, </2, ••• ,<7п так,
чтобы
п	п
А =	= inf II/-22^11-
г=1	г=1
Элемент
п
h=z^2 c^9i
1=1
называется элементом наилучшего приближения, а величина А называется
величиной наилучшего приближения. Далее мы будем предполагать, что f
не принадлежит линейной оболочке <71, <72, • • •, <7П> то есть А > 0. Задача
отыскания такого элемента решается достаточно просто, если в простран-
стве задано скалярное произведение.
Теорема 1. Вещественное гильбертово пространство строго нормирован-
но.
Доказательство. Напомним, что пространство называется строго норми-
рованным, если равенство
Ik + sll = Ikll + Ikll	(13.1)
выполняется тогда и только тогда, когда у = Хх Л > 0. В одну сторону,
если у = Хх, Л > 0 неравенство треугольника очевидно превращается
в равенство. Пусть теперь неравенство треугольника превратилось в равен-
ство. Докажем, что существует такое Л > 0 , что у = Хх. Заметим, что из
||я + г/|| = ||ж|| + ||2/|| следует, что
(*,</) = 1И1 Ы-	(13.2)
В самом деле, возводя (13.1) в квадрат, получаем
lk + 2/Ц2 = (х + у,х + у) = ||z||2 + Ц1/112+ 2(ж,?/) = ||ж||2+ ||i/||2 + 2||z|| ||j/||.
В случае строгого неравенства треугольника
(я,1/) < ||х|| ||2/||.
72
Рассмотрим ||г/ — Лж||:
||Х/ - Ла?||* 2 * * * * 7 = (у - Ат, у - Хх) = А2||ят||2 - 2А(®, у) + ||?/||2.
Последнее выражение есть многочлен второй степени относительно Л.
Найдем его корни:
Л2||я||2 - 2Х(х,у) 4- ||з/||2 = 0.
В силу (13.2) дискриминант этого уравнения равен нулю и
А = ГК > °-
ГН
При выполнении строгого неравенства треугольника дискриминант отрица-
тельный и уравнение вещественных корней не имеет. Таким образом пока-
зано, что при выполнении (13.1) и только в этом случае существует Л > 0,
что \\у — Ат|| = 0. Теорема доказана. Эта теорема справедлива и для ком-
плексного гильбертова пространства. Доказательство проводится по той же
схеме, что предлагается проделать самостоятельно.
Теорема 2. В гильбертовом пространстве элемент наилучшего прибли-
жения единственнен.
Доказательство. Предположим противное: существуют два элемента наи-
лучшего приближения Tii и /г.2, такие, что \\f — /г.1|| = Д и \\f — /12Ц = Д-
Тогда
II/-мп. II/-м
2	2
Но
II/-	> Д
по определению величины наилучшего приближения. Следовательно,
f~h2 ||/-МН ||/-М
2	2	2	2'
Поскольку гильбертово пространство строго нормированно,
/ ~ М _ х / ~ М
2	2 ‘
Если Л = 1 , то /г.1 = /?2- Если же А 7^ 1, то
«__h\ — Л/?2
7 = 1-Л ‘
73
Ho hi и /г.2 есть линейная комбинация элементов gi, дъ,..., дп, и тогда и f
есть линейная комбинация этих элементов, т. е. принадлежит их линейной
оболочке. Получено противоречие. Теорема доказана.
Следующие две теоремы дают практический аппарат построения наилуч-
шего приближения в гильбертовом пространстве.
Теорема 3. Пусть в гильбертовом пространстве Н существует элемент
ho, реализующий наилучшее приближение к элементу f. Тогда разность
f — ho ортогональна всем элементам гильбертова пространства Н.
Доказательство. Предположим противное, что в Н есть элемент х,
такой, что
(/ - ho, х) = а 0.
Без ограничения общности можно считать, что ||ж|| = 1. Пусть у = ho—ах.
Тогда
11/-//Ц2-
(f -ho- ах, f - ho- ax') = \\f - ho ||2 + aa|!z||2 - a(x, f - h0) - a(f - ho, x).
Ho (f — ho,x) = а и ||ж|| — 1. Поэтому
II/ “ ?/||2 = II/ “ ho||2 + aa - aa - aa.
В итоге получаем
\\f-y\\2<\\f-ho\\2.
Последнее неравенство означает, что ho не есть элемент наилучшего при-
ближения. Это противоречие и доказывает теорему.
Теорема 4. Пусть (/ — ho,x) = 0 для любого х 6 Н. Тогда х есть
элемент наилучшего приближения.
Доказательство. Для произвольного х G Н имеем
ll/-z||2 = (f-ho + ho-xj-ho + ho-x) = (f-h0, f-ho) + (f-ho,ho-x)+
+(/i0 - x, f - ho) + (h0 - x,h0- x).
Ho f—ho ортогонально всем элементам из Н, в том числе и ho—х. Поэтому
второе и третье слагаемые обращаются в нуль. В итоге имеем
II/- < = II/ - holl2 + ||Ло - < > ||/ - Ло||2.
Это неравенство и доказывает теорему.
74
Эти две теоремы дают практический алгоритм нахождения наилучшего
приближения для заданного элемента f. В самом деле, пусть
п
Hq — Cj9i'
г=1
Умножая скалярно выражение
п
f -^Ci9i
г=1
на gi г = 1, 2,..., п, получаем систему для определения коэффициентов
Cj.
Cl(^l,gi)+c2(p2,Pl) + --- + Cn(^n,Pl) = (Л Pl)
Ci(pi, рг) + с2(рг, Рг) + •. • + Сп(рп, 92} = (/, Рг)
ci(pi,рп) + с2(р2, дп} + ... + сп(дп, дп) = (J, дп}.
Матрица этой системы есть матрица Грамма. Если система pi, р2,..., дп ли-
нейно независима, то ее определитель отличен от нуля. Поэтому полученная
система разрешима единственным образом.
Рассмотрим пример. Пусть в пространстве 2>2,(о,4) со скалярным произ-
ведением
4
(f,9) = У f(x}g(x}dx
о
для функции f(x} =	— 1| нужно найти наилучшее приближение мно-
гочленом первой степени. Для элемента наилучшего приближения ho(x} —
+ а^х получаем систему для определения коэффициентов во, ар
«о(Р1,Р1) + «1(Р2,Р1) = (к “ !|,Р1)
ao(pi,p2) +ai(p2,P2) = (к ~ 11,Р2)-
Здесь gk = хк, к = 0,1, а
4	4
(Рг,Р>) = / ХгХ^Х (к~ 1,Рг) = / к — l\xZdx.
о	о
75
Вычисляя интегралы, получаем систему для определения коэффициентов
многочлена наилучшего приближения.
	4oq 4* 8ai — 5, 64	2 8«о +	= 13-. О	о
Отсюда	1 11 ao--g, а1-уё’
и многочлен имеет вид	м . 11 1
76
§ 14. Наилучшее равномерное приближение
Рассмотрим пространство С вещественных непрерывных функций, за-
данных на отрезке [а, Ь], с нормой
[|/|| = max |/(z)|
[а,6]
Будем искать наилучшее приближение заданной функции f(x) многочле-
нами степени п. Далее будем предполагать, что f(x) не есть многочлен
степени п .
В пространстве С норма элемента определяется не через скалярное
произведение, что значительно осложняет задачу построения многочлена
наилучшего приближения.
Итак, для заданной функции /(я?) будем искать многочлен наилучшего
приближения
п
^?n(*^)	Qk% >
к=0
такой, что для любых многочленов Рп(х) степени п имеет место
En(f) = \\f-Qn\\<\\f-Pn\\.
Многочлен Qn(^) будем называть многочленом наилучшего равномерного
приближения, а величину En(f) — величиной наилучшего приближения.
Теорема 1 (Валле - Пуссен). Пусть существуют п 4- 2 точки
а < хо < Xi < ... < хп+\ < Ь, такие, что величина /(ж) — Рп(Е) после-
довательно меняет знак при переходе от точки Xj к а;г-+1, где Рп(х) —
некоторый многочлен степени п, и пусть
V = min \f(xi) - Pn(xi)|.
0<г<п+1
Тогда En(f) > р.
Доказательство. Пусть ц > 0. (Случай р = 0 —> Еп > 0 тривиа-
лен.) Предположим противное, т. е. En(f) = \\f — Qn|| < щ и рассмотрим
многочлен
Тп(ж) = Рп(х) — Qn(E)
77
в точках xi i = 0,1,..., п + 1. Здесь Qn(x) — многочлен наилучшего
равномерного приближения. Имеем
Рп{х) - Qn(x) = (Рп(х) - /(®)) - (Qn(z) - /(ж)).
Поскольку Qn(x) — многочлен наилучшего равномерного приближения,
I/W - QnMI < I/W - PnMI
во всех точках отрезка [а, Ь]. Поэтому знак Pn(x) — Qn(x) в точках Xi i =
0,1,... , n + 1 определяется знаком выражения (Рп(х) — f(x)). Но в этих
точках последнее выражение п + 2 раза меняет знак. Мы получили, что
многочлен Тп(х) степени п меняет знак п+2 раза и, следовательно, имеет
п+1 корень, что невозможно. Теорема доказана.
Теорема Чебышева. Для того чтобы многочлен Qn(x) был многочленом
наилучшего равномерного приближения непрерывной функции /(ж), необ-
ходимо и достаточно существование на отрезке [а, 6], по крайней мере,
п + 2 точек xq < Xi < ... < xn+i таких, что
f(xk) - Qn^Xk) = а(-l)n|| f - Qn|| к = 0,1,..., n + 1
a — 1 или a = — 1 одновременно для всех к.
Точки, удовлетворяющие условиям теоремы, называются точками чебы-
шевского альтернанса.
Иными словами, теорема утверждает, что для многочлена наилучшего
равномерного приближения в точках чебышевского альтернанса разность
/(я) — Qn(x) должна по абсолютной величине достигать своего максималь-
ного значения, равного \\f — Qn||, и менять знак при последовательном
переходе от точки хк к точке Xk+i-
Доказательство. Мы докажем только достаточность. С доказательством
необходимости можно ознакомиться в [2]. Пусть Qn(x) — некоторый много-
член степени п, для которого в точках чебышевского альтернанса выполнено
условие теоремы, и пусть L = \\f — <5п||- Согласно теореме Валле — Пус-
сена L = р, < En(f). Но, по определению, En(f) < L. Следовательно,
En(f) = L и многочлен Qn(rr) есть многочлен наилучшего равномерного
приближения.
В случае гильбертова пространства была доказана теорема единственно-
сти. Докажем эту теорему для пространства С.
Теорема. Многочлен наилучшего приближения непрерывной функции еди
ственнен.
78
Доказательство. Предположим противное; существуют два многочлена
наилучшего приближения
Qn(x) ± QnW такие> чт0 II/ “ Qill = II/ ~ Qnll = Еп(П.
Тогда
следовательно, многочлен
также является многочленом наилучшего равномерного приближения. Пусть
то, Xi,..., xn+i — точки чебышевского альтернанса многочлена <?„(#). В
этих точках
ж)_Ш+Ш=ад	fc = 0,l,...,n + l
XU
или
|(/(zfc) - Qn(xk)) + f(xk) -Q2(Tfc)| = 2En(/) fc = 0,1,..., n 4-1.
Но в этих точках
\f(xk)-Q^xk)\<EM	г = 1,2.
Поэтому последнее равенство возможно, лишь если
f(xk) - Q„(xk) = f(xk) - Q2n(xk) к = 0,1,..., п + 1.
Но это означает, что два многочлена совпадают в п + 2 точках, что воз-
можно лишь в случае, если они равны. Теорема доказана.
Рассмотрим примеры наилучшего равномерного приближения.
1. Непрерывная функция /(т) приближается многочленом нулевой сте-
пени.
Пусть М = max f{x) и т = тт/(т). Тогда таким многочленом будет
[а,6]	[а,Ь]
многочлен
М — т
у = —2~'
поскольку для него есть две точки чебышевского альтернанса xi и а?2,
такие, что /(#i) = т и /(^г) = М.
79
2. f(x) = |rr — 1| на отрезке [0,4] приближается многочленом наилуч-
шего равномерного приближения первой степени у = ах 4- Ь.
Пусть /(0) — b = d, /(1) — (а 4-6) = — d, /(4) — (4а4- &) = d. Очевидно,
что точки zi = 0, Х2 = 1 и хз = 4 являются точками чебышевского
альтернанса. Решая систему
1 — b = d, — (а 4- b) — —d, 3 — (4а 4- Ь) = б/,
находим а = 1/2, b — 3/4 и многочлен первой степени наилучшего
равномерного приближения имеет вид
х 1
у 2	4
Сравните с многочленом, полученным в § 12 для этой же функции в гиль-
бертовом пространстве L2 на том же отрезке.
Эта же методика построения многочлена наилучшего равномерного при-
ближения 1-й степени может быть применена в случае, если функция f{x)
есть многочлен 2-й степени.
3. Выпуклая (вогнутая) дифференцируемая функция f(x) на отрезке
[а,Ь] приближается многочленом 1-й степени.
Построим прямую, проходящую через точки (a,f(a)) и (b.f(b)). Уравнение
этой прямой
,, ДО-Да),
У = --Г-----(z - о) + JW-
о — а
Найдем точку xi, касательная в которой к графику функции f(x) парал-
лельна к этой прямой, т. е. решим уравнение
гы-Ы
Ь CL
Уравнение касательной в этой точке
о — а
Проведем теперь третью прямую, параллельную касательной на равном рас-
стоянии от первой прямой и касательной. Уравнение этой прямой и будет ис-
комым многочленом первой степени наилучшего равномерного приблжения.
В самом деле, для этого многочлена есть три точки чебышевского альтер-
нанса a, Xi, Ь, где разность f(x) — РДя) максимальна и меняет знак. Най-
дем уравнение этой прямой. Для этого достаточно найти точки пересечения
80
1	г	о — а
Хл = — а 4- xi-—-----
2	1 f(b )-f(
(f(a) - f(X1 ) .
первой прямой и касательной с осью ОХ и провести прямую, проходящую
через середину полученного отрезка с заданным угловым коэффициентом
к=т-ж>
b — а
Первая прямая пересекается с осью ОХ в точке
_ г/ \ Ь~а
Х2~а
а касательная — в точке
Точка Х4 пересечения искомой прямой с осью Ох есть полусумма точек
Х2 И Хз
«)
Наконец, уравнение прямой, пересекающей ось ОХ в точке Х4 с угловым
коэффициентом к, имеет вид
,. ffiizza,, _
о — а
Изложенная методика позволяет построить многочлен 1-й степени наи-
лучшего равномерного приближения для дифференцируемой функции f(x) и
достигающей максимума (минимума) во внутренней точке отрезка [а,Ь]. По-
строение этого алгоритма предлагается проделать самостоятельно.
81
§ 15. Численные методы решения задачи Коши
для обыкновенных дифференциальных уравнений и систем
Пусть на отрезке [to,?1] требуется найти решение дифференциального
уравнения
x = f(t,x),	(15.1)
удовлетворяющего начальному условию
x(t0) = х0.	(15.2)
Разобьем отрезок [to, Т] на отрезки [tj, ty+i], j = 0,1,..., п — 1, tj+i — tj =
h и рассмотрим алгоритмы приближенного нахождения Xj = x(tj) —
значения решения задачи (15.1)— (15.2). Простейшим методом нахождения
приближенного решения задачи Коши является метод Эйлера, согласно ко-
торому значение решения в точке tj+i находится по формуле
Xj+i = Xj + hf(tj,Xj).	(15.3)
Эту формулу можно получить из следующих соображений. Разложим Xj+i =
x(tj + Л) в ряд по степеням h, ограничиваясь лишь линейными членами
разложения
x(tj + h) = x(tj) + hx(tj) + <9(Д2).
Так как i(tj) = f(tj,Xj), получаем формулу (15.3). При этом замечаем,
что прогрешность на шаге имеет величину порядка h2. Для практиче-
ских расчетов эта формула применяется достаточно редко из-за высокой
погрешности на шаге. Если, например, требуется найти решение на отрез-
ке длины 100 с шагом 0.01, то при расчете по методу Эйлера значений в
конце интервала итоговая погрешность будет порядка единицы, что обычно
неприемлемо. Более того, если искомое решение неустойчиво по Ляпунову,
то расчет решений с помощью такого грубого метода просто теряет смысл.
Эти соображения привели к построению более точных методов численного
интегрирования дифференциальных уравнений.
Определение. Метод интегрирования дифференциальных уравнений бу-
дем называть методом к-го порядка, если погрешность на шаге порядка
O(hk+1).
Рассматриваемое ниже семейство методов носит название методов Рунге
— Кутта.
82
Рассмотрим равенство
h
x(t + h) = x(t) 4- x(t + r)dr.	(15.4)
о
Воспользуемся формулами приближенного вычисления интеграла, напри-
мер формулой трапеций
x(t + Д) = x(t) 4- ^[i(i) 4- x(t 4- h,)] 4- О(Д3)
или
x(t 4- h) = x(t) 4-	x(t)) 4- f(t 4- h, x(t 4- h))] 4- O(7i3),
XU
что приводит к формуле
•£j+i Xj 4-	Xj) 4- f(tj+i, ж?4-1)]
xu
с погрешностью на шаге порядка <9(/i3). Однако полученная расчетная
формула неприменима, поскольку правая часть неявно зависит от Xj+i.
Заменим в правой части Xj+i на некоторую другую известную величину,
такую, чтобы итоговая погрешность на шаге по-прежнему была порядка
О(№). Если в качестве такой величины взять х*, вычисленную по методу
Эйлера
х* Xj 4- hf(tj, Xj),
то итоговая погрешность по-прежнему останется порядка O(h?). В самом
деле, x(t 4-/г) — х* = O(h2) и по теореме Лагранжа о конечных приращениях
+ h,x(t + h)) - f(t + h, X*) =	+ ft) - x").
ox
В результате этих рассуждений приходим к расчетным формулам
Qi — hf (j'ki хк) > Q2 — hf (tk	h, хк 4” Qi)>
. Qi + Q2	Z1K r\
Xk+1 = Xk +---------•	(15.5)
Вычисляя интеграл (15.4) по формуле прямоугольников
^j+i — xj 4"	4" о> x(tj 4" о))
83
и заменяя x(tj 4- ^) на
х — xj
вычисленную по формуле Эйлера, получаем другие расчетные формулы
Qi = hf(tk,xk), Q2 = hf(tk +	+ ^)’
%к+1 = Хк + <?2-	(15.6)
В случае системы дифференциальных уравнений, например
x(t) = fi(t,x,y)	x(t0) = я?о
y(t) = f2(t,x,y)	y(to) = yo,
те же соображения приводят к формулам
Qi = hfi(tk, хк.ук), Q\ = /1/2(4, хк.ук),
Q2 = hfi{tk + h,xk + Qf, ук + Q\),
Q2 = hf2(tk + h,xk + QI, yk + Ql),
SM=yt + ^.	(15.5').
Формула (15.6) в случае системы двух уравнений имеет вид
Qi = hfi(tk, хк.ук), Q\ = hf2(tk, хк.ук),
Q2 — hfi(tk +2’ Хк + Ук +
h	Ох	Оу
Qy2 = hf2(tk + ^, хк + ^, ук + ^),
Zl	Zl
2-/;+1 — Хк 4" Q2i
У к+1 =yk + Q2-	(15.6')
Полученные формулы реализуют методы Рунге — Кутта 2-го порядка. J
Изложим схему получения вычислительных формул для методов более
высокого порядка. Будем последовательно вычислять
Qi(h) = hf(t,x),
84
Q2(/i) = hf(t + a2h,x +fiziQilJi)),
Qp(h) = hf(t + aph, x + /3piQi(h) + ... + Z^p,p-i Qp-iW)
и положим
p
x(t + h) « z(h) = x(t) + diQi(h).
Параметры 09, Д-j, di выберем так, чтобы величина = y(t + h) — г(Д)
была порядка O(hs+1). Величина s называется порядком метода. Условие,
наложенное на <p(h), означает, что 92(0) = <//(0) = ... = 92^(0) = 0, а
^(з+1)(о) ^о.
Пусть р=1. В этом случае z(Jt) = x(t) + d]hf(t, х) и
<p(h} = x(t + ti) — x(t) — dihf(t, x).
Тогда
92(0) = 0,	</>'(0) = x'(t) - di/(t, x) = f(t, z)(l - di).
Для того чтобы равенство 92z(0) = 0 выполнялось при любой функции f(t, ж),
надо, чтобы di = 1, и мы получаем формулу Эйлера
x(t + /г) = x(t) + hf(t, х).
Построим формулы 2-го порядка. Пусть р=2. В этом случае
z(ti) = x(t)+diQi[h)+d2Q2(h) = x(t)+dihf (t, x)+d2hf (t+o^h, x+faihfit, x))
и
9?(/i) = x(t + A) — x(t) — dihf(t, x) — d2hf(t + a2h, x 4- 02ihf(t, x\).
Введем обозначения
t = t + &2h, x = x + folk f(t,x).
Вычислим 9?'(^) и ^(O) :
ip'fji) = x\t + Д) - x) - d2f(t, x) - d2h[a2f'fit, x) + Д21Л(Л x)],
^"(A) = x"(t + h) - 2d2[a2f^t, x) + Д21 f'x(i, x)f(t, ж)] - d2h[ . . . ].
Поскольку мы будем вычислять 9/(0) и «//'(О), выражение в последних квад-
ратных скобках для нас интереса не представляет.
85
При h = 0 получаем
9?(0) = x(t) — z(£) = О,
<р'(°) = x'(t) - dif(t,x) - d2f(t,x) = f(t,x)(l - dr- d2).
Далее, поскольку x"(t) = ft(t, x) 4- fx(t, x)x'(t) = f[(t,x) + fx(t,x)f(t,x)
4>"(h) = fttt’x) + /Ж x) ~ Zd2a2f{(t, x) - 2d2fx(t, x)f(t, x).
Вычислим эти выражения при h 0.
(/>(0) = x(tO — x(x) = 0,
¥?'(0) = (1 - di - d2)f(t, x),
</(0) = (1 - 2d2a2)f[(t, x) + (1 - 2d2/32i)f'x{t, x)f(t, x).
Равенство нулю </(0) и </?"(0) при произвольной функции f(t, х) возможно,
если
1 — di — d2 = 0,
1 — 2d2a2 = 0,
1 — 2d2/32i = 0.
Мы получили три уравнения относительно четырех неизвестных. Задавая
один из параметров, получаем семейство методов Рунге — Кутта 2-го по-
рядка^ В частности, при
di = d2 = -, а2 = (321 = 1
получаем формулу (15.5)
Qi = hf(tk,xk), Q2 = hf(tk + h,xk + Qi), xk+1 = xk +
Формула (15.6) получается при
di = 0, ‘d2 — 1,	= Д21 = x-
ла
h Q1
Ql	hf (j'ky ‘%k) ?	Q"2 — hf (th 4" T? %k 4 n~) ’	^Ar+1 =	4“ Q2*
Мы построили методы второго порядка. Методы более высокого порядка
получаются при р > 2. Построение расчетных формул для методов связано с
86
громоздкими вычислениями. Приведем наиболее употребительные формулы
для методов 3-го и 4-го порядка.
Методы 3-го порядка.
1.
h	Q
Qi = hf(tk,xk), Q2 =	+	= hf(tk + h,xk-Q1 + 2Q2),
Z	z
_ Qi + 4Q2 + Q3.
^k 4”
(15-7)
2.
Qi = hf(tk,Xk), Qi = hf(tk+^, xk+^), Qi = hf(tk+^, xk-Qi + ^—),
00	00
xk+i — xk
Методы 4-го порядка.
1.
Qi = hf(tk,Xk),
h	О
Qz — hf(tk + —,xk + —
z	z
= %k
Qi + 3Q3
4
Q2 = hf(tk + —,xk + —),
z z
Q4 = hf(tk 4- /1, xk -|- Q3),
2Q2 + 2Q3 + Q4
6	’
(15.8)
(15.9)
Напомним, что в случае системы
У№ = h(t,x,y)
87
величины Qi становятся двумерными векторами Qi = (Qf, Q^)T, причем
при вычислении Q^ используется функция fi(t,x,y), а при вычислении
Qf — функция f2(t,x,y) так же, как и при получении формул (15.5х) и
(15.6') для систем из формул (15.5) и (15.6) для скалярного уравнения.
В заключение отметим, что методы Рунге — Кутта относятся к так на-
зываемым одношаговым методам. В одношаговых методах решение в точке
ti = to + h строится с помощью информации о решении, полученной на
предыдущем шаге. В многошаговых методах при вычислении Xk использу-
ется информация о решении, полученная на нескольких предыдущих шагах.
Опишем коротко одну из процедур получения двушагового метода. Пусть
для простоты fi+i = ti + h. Первый шаг сделаем по одному из методов Рунге
— Кутта, вычислив xi в точке ti = to 4- h. Построим интерполяционный
многочлен в форме Ньютона по двукратой точке to и t\ :
L3(t) = xq + xo(t — to) 4—-—-^2-~ *o)2-
FL
Положим <2 — ti 4- h и будем считать, что Х2 = x(ti 4- h) = L$(ti 4- h) .
Получим
X2 = —Зжо 4- 4zi — 2hio
и мы приходим к формуле
Xk+i = -3xk-i 4- 4zfc - 2hxk-i-
Из оценки интерполяционного многочлена следует, что построенный метод
есть метод 2-го порядка.
Для увеличения точности нужно строить интерполяционный многочлен
более высокого порядка, сделав предварительно несколько шагов по методу
Рунге — Кутта. Естественно, что при увеличении числа точек для построе-
ния интерполяционного многочлена формула становится более громоздкой.
Кроме того, какой смысл переходить к новым расчетным формулам, а не
продолжать расчеты по методу Рунге — Кутта? Все эти причины и при-
водят к тому, что многошаговые методы используются в вычислительной
практике значительно реже, чем методы Рунге — Кутта.
88
§ 16. Разностные методы решения краевых задач
для обыкновенных дифференциальных уравнений 2-го порядка
Понятие аппроксимации, устойчивости и сходимости разностных схем
Рассмотрим следующую задачу. Требуется найти решение дифференци-
ального уравнения
у" + р(х)у' + q(x)y - f(x) 0 < х < 1,	(16.1)
удовлетворяющее на концах отрезка краевым условиям вида
«12/(0) + &2/'(0) = 71,
«22/(1) + /?2?/'(1) = 72-	(16.2)
Будем говорить, что уравнение (16.1) и краевые условия (16.2) задают в
пространстве дважды непрерывно дифференцируемых функций дифферен-
циальный оператор
Ly = f-
Разобьем отрезок [0,1] точками х^
Xi = a + ih, г = 1,2,	h = —
п
и аппроксимируем в точках х^ i = 0,2,..., п производную У(жг) выражени-
ем
у'(х{) ыУг+1^ Уг, где yi = y(xi),	(16.3)
или выражением
// \ Уг Vi—1
у	—h—	(16-4)
Введем обозначения
Ayi = yi+1 ~ Уг, ^yi = yi~ yi-i,
yi+\ - yi
l/r i — —— = -----------.
- ^yi - Vi ~ yi~x
y$,i h ~ h
89
Второй производной y"(xi) поставим в соответствие выражение
_	_ АУ?/г _ УДуг _ Уг-l - 2^ + уМ
Ухх,г Ухх,г д2	^2	Jj2
и дифференциальному уравнению (16.1) в точках Xj можно поставить в
соответствие следующее уравнение, которое будем называть разностным
yi—1 tyi “Ь yi+1 . yi+1 Уг .	г -io	1	/1й
------тъ	+Pi	т-+ QiVi = fi г = 1, 2...., n — 1	(16.5)
rr-----------------------------------------------------------а
После приведения подобных членов уравнение (16.5) запишем в виде
kyi+i - Ciyi + diyi-x = -ipi i = 1,2.... ,n - 1,	(16.6)
а краевые условия в виде
Уо = Х1У1 + Pi Уп = Х2Уп-1 + М2-	(16-7)
Здесь
bi = l-\-pih, c-i — 2 + pih - qih2, di = 1, Pi = fih2,
==	_	72^
X1	a\h-(3\	+	a2h + /?2
При аппроксимации y'(xn) в правом граничном условии использовались
формулы (16.4).
Задачу (16.6)— (16.7) будем в дальнейшем называть разностной или се-
точной краевой задачей, а функцию целочисленного аргумента yi — сеточ-
ной функцией. Эта задача представляет собой с учетом граничных усло-
вий (16.7) систему п+1 уравнений относительно п+1 неизвестного yi i —
0,1,..., п. Ниже будут рассмотрены методы решения таких задач, однако,
прежде чем приступить к их рассмотрению, следует выяснить, насколько
решение непрерывной задачи (16.1) - (16.2) отличается от решения сеточной
задачи (16.6) - (16.7).
Пусть щ = u(xi} — значение решения задачи (16.1) - (16.2) в узле xi, а
yi — значение решения сеточной задачи в том же узле.
Определение 1. Будем говорить, что решение задачи (16.6) - (16.7) схо-
дится (сходится со скоростью hk ), если Zi = yi — щ —> 0(^ = O(hk)} при
h —> 0.
Обозначим задачу (16.1) - (16.2) оператором
Lu = f,
90
а задачу (16.6) - (16.7) — оператором
lyi = 4>i-
Пусть по-прежнему Z{ = yi — щ, т. е. Уг = Zi + щ. Тогда для Zj получаем
задачу
lZi — Фг,
Где фг = <Pi~ 1Щ.
Определение 2. Будем говорить, что разностная задача lyi = <pi ап-
проксимирует задачу Lu = f (аппроксимирует с порядком hk), если для
Zi = yi — т в задаче Izi = фг для правой части имеет место оценка
Л = О(Л‘).
Определение 3. Будем говорить, что разностная схема
lyi = <Рг
устойчива, если имеет место оценка
Ы < .4W,
причем постоянная А не зависит от h. Если А зависит от h, схема называется
условно устойчивой.
Из аппроксимации разностной схемы и ее устойчивости очевидно следует
ее сходимость.
Выясним порядок аппроксимации разностной схемой (16.6) - (16.7) исход-
ной задачи (16.1) - (16.2). Пусть u(xj) — решение задачи (16.1) - (16.2) в
узле Х{. Тогда
и
д2	Ь3	Л4
(Xi ± h) = u(xi) ± hu'(xi) + -^u"(xi) ± -Qu'”(xi) + — u^\xi) + O(/i5).
Поэтому
^г+1	//	\ . z-)/r\
----------= и [Xi) + <9(Д),
f L
Щ+1	4- Ui—1	/t. .
h2
и левая часть в (16.5) есть
u"(xi) +Piu'(xi) + qiU^Xi) + 0(h).
Таким образом, задача (16.6) - (16.7) аппроксимирует задачу (16.1) - (16.2)
с порядком O(h), причем если р(х) = 0 и Д = /?2 = 0, то порядок
аппроксимации O(h?).
91
§ 17. Методы решения краевых задач 2-го порядка.
Методы прогонки, метод стрельбы
Как отмечалось, краевая задача
biyi+1 - СгУг + diyi-1 = -(pi i = 1,2.	— 1	(17.1)
Уо = Х1У1 + Mi Уп = Х2Уп-1 + М2	(17.2)
представляет собой линейную систему п -I- 1 уравнения относительно неиз-
вестных yo,yi,... ,уп. Особенностью матрицы этой системы является то,
что эта матрица является трехдиагональной, т. е. у этой матрицы отлич-
ны от нуля только элементы на главной и двух соседних диагоналях. Ниже
будут предложены простые методы решения таких систем.
Метод прямой (правой) прогонки.
Теорема 1. Решение задачи (17.1) - (17.2) может быть представлено в
виде
Уг = Лг+12/i+i + Pi+I i = П -- 1,п - 2,..., 1,0,	(17.3)
где
а «г и Д определяются из рекурентных соотношений
Ьг
«г+1 =--------Т г - 1,2,... ,п - 1,	«1 - XI,	(17.5)
Сг ai dj
A+i = di^_+	г = 1,2,..., п — 1,	/31 = й. (17.6)
Cj di
Доказательство. Левое граничное условие у^ = Х1У1 + Mi имеет вид (17.3),
где ai = xi , a /?i = Ml- Уравнение (17.1) при i=l
Ь1У2 ~ сгуг + с^уо =
после подстановки уо приобретает вид
Ъ1У2 ~ С1Уг + ск(х1У1 + Mi) = -<М1-
Или
_	di + pifii
У1 —------—У 2 Ч-----------
ci - 01X1	- otidi
92
т. е. при i=l уравнение (17.1) представимо в виде (17.3), где и Д имеют
вид (17.5), (17.6). Предположим теперь, что i — 1-е уравнение (17.1) пред-
ставимо в виде (17.3). Подставим yi-i в г-е уравнение (17.1), получим
biyi+i - ауг 4- di^yi 4- Д) =
или
difii 4" <-Pi
Уг = -----V&+1 + --------V-
oiidi	Cj aidi
Следовательно, и i-e уравнение (17.1) представимо в виде (17.3) и коэффици-
енты и Д представимы в виде (17.5), (17.6). Последнее уравнение (17.3)
и правое граничное условие представляют собой систему двух уравнений
отностельно уп-1 и уп
Уп = Х2Уп—1 + Р>2
Уп—1 = ОСпУп “Ь Дг?
решая которую получаем (17.4).
Теорема доказана. Эта теорема дает удобный алгоритм нахождения решения
задачи (17.1) - (17.2). В самом деле, левое граничное условие задает ai
и Д. Далее, в цикле i = 1,2, ...,п — 1 определяются про формулам
(17.5), (17.6) прогоночные коэффициенты cej+i и Д+ь Затем по (17.4)
определяется уп и далее обратным ходом по формулам (17.3) определяются
yi i = п — 1, п — 2,..., 0.
Рассмотренный в теореме 1 метод, при котором определение уг про-
изводится последовательно справа налево, называют методом прямой или
правой прогонки. Аналогично доказывается
Теорема 2 (метод обратной или левой прогонки). Решение задачи (17.1)
(17.2) представимо в виде
Уш = &+1У1 4- yi+i i = 0, l,...,n- 1,	(17.7)
где
__ di
Ci bi^i+l
^i^i+l +
Vi —
Ci *гьг+1
Д1 + Х1У1 Уо=	c , 1 - 6X1		(17.8)
i = n — 1, n — 2,..., 1,	Cn — X2y	(17-9)
г = n — 1, n — 2,..., 1,	Уп = Д2-	(17.10)
93
Иногда используется метод встречной прогонки, который состоит в следу-
ющем. Зафиксируем некоторое i = I и, используя метод прямой прогонки,
определим щ и Д i = 1, 2,..., I + 1. Затем, используя метод обратной
прогонки, определим и i = п — 1, п — 2,..., I + 1. Далее, из системы
У1 = oii+iyi+i + (31+1
Ум = Cz+iS/z + Vm
определим yi и j//+i и обратным ходом определим остальные yj.
Очевидно, что методы прогонки применимы, если в процессе вычислений
знаменатели в формулах (17.4) — (17.6) и аналогичных им в методе обрат-
ной прогонки не обращаются в нуль. Достаточный признак применимости
метода прямой прогонки фомулируется в следующей теореме.
Теорема 3. При выполнении условий
|ci| > \bi\ + \di\ г = 1,2,..., n — 1	(17.11)
|X1| < 1,	|Х2| < 1,
причем одно из этих неравенств строгое, знаменатели в формулах (17.4) -
(17.6) не обращаются в нуль и |ct,| < 1, i — 1,..., п.
Доказательство. Заметим, что <то = Х1 и по условию |ао| < 1- Предпо-
ложим, что |аД < 1, и покажем, что и |с*г-ь1| < 1- Рассмотрим формулу
_ bi
^г4-1	j
Ci
и покажем, что — <4а?г| > |Ьг|. В самом деле, используя (17.11) и индук-
тивное предположение, получаем
|Сг - di&i\ - \bi| > |Cj| - |б?г|lofi) -	> Щ1 - |cv,|) > О,
что и доказывает, что |с^ — di(Xj\ > |6г| и |агг+1| < 1.
Обратимся к формуле (17.4). Заметим, что если |сц>| < 1, то |сег| < 1,
при i > к. Поэтому, если |xi| <1, то и |an| < 1 и знаменатель в (17.4)
не обращается в нуль. Если же |an| < 1, то условие |хг| < 1 гарантирует
отличие от нуля знаменателя в (17.4).
Теорема доказана.
Метод стрельбы. Рассмотрим систему уравнений (17.1)
ЬгУм ~ Ciyi + diy^x =-ipi г = 1, 2...., n — 1.
94
Последовательность этих уравнений можно рассматривать как систему п—1
уравнений относительно п + 1 неизвестного уо, yi,..., yn+i- Рассмотрим
сначала однородную систему
biyi+1 - СгУг + diyi-i = 0 i = 1, 2...., п - 1.	(17.12)
Очевидно, что разность двух решений неоднородной системы есть решение
однородной системы. Поэтому, для того чтобы описать все решения неод-
нородной системы, надо к формуле общего решения однородной системы
прибавить любое частное решение неоднородной системы (17.1). Мы будем
предполагать, что коэффициенты системы (17.12) 6j, Cj, di таковы, что ранг
этой системы в точности равен п — 1. Тогда эта система имеет два линей-
но независимых векторных решения — вектора U = (uq,ui, ... ,un+i)T и
V = (г?о, гц, • • •, vn+i)T . Если U и V есть решения (17.12) и если поло-
жить uq = 1, а щ = 0, то последовательность уравнений (17.12) стано-
вится рекурентной для определения U2, из,..., ип . Аналогично, если поло-
жить Vq = 0, Vi — 1, то также из (17.12) можно рекурентно определить
V2.. г>з,... ,vn. Очевидно, что два этих вектора линейно независимы и их ли-
нейная комбинация
Citf + CaV
есть общее решение однородной системы (17.12). Частное решение неодно-
родной системы (17.1) — вектор W = (wq,wi, ..., wn)T также рекурентно
можно определить из (17.1), положив Wq = 0, гщ — 0. Таким образом,
общее решение (17.1) можно представить в виде
yi = ciUi 4- c2Vi 4- Wi.	(17.13)
Определим теперь ci и сг так, чтобы решение удовлетворяло краевым
условиям (17.2)
Уо = Х1У1 + Mb Уп = Х2Уп-1 + М2-
*
Из левого краевого условия получаем
О = Х1С2 + Д1,
из правого краевого условия —
ci«n + c2vn + wn = X2(ciun_i 4- c2fn_i 4- wn_i) 4- ц2-
Таким образом, получена система
95
Cl = X1C2 4- //1
С1(«Ц- X2^n-1) + C2(vn - X2*>n-1) = ^2 + X2Wn-i - Wn
Найдем из этой системы ci и с2 и подставим в (17.13). Получим решение
уравнения (17.1), удовлетворяющее краевым условиям (17.2).
Мы изучили методы решения краевых задач для разностных уравнений.
Но для того, чтобы решение разностной схемы сходилось к решению исход-
ной задачи, необходимо, чтобы эта разностная схема была устойчивой. Для
исследования на устойчивость нужно получить оценку решения разностной
задачи через правую часть. Для этого исследуем некоторые общие свойства
разностных уравнений.
96
§ 18. Принцип максимума для разностных уравнений 2-го порядка.
Оценка решений первой краевой задачи
Рассмотрим следующую краевую задачу для разностного уравнения
2-го порядка:
lyi = biyi+i - CiVi + diyi-i = -(Pi г = 1,2,..., n-1,	(18.1)
Уо = М, Уп = Р2,	(18.2)
причем
bi > 0, di > 0, Ci > bi + di.	(18.3)
Теорема 1 (Принцип максимума). Пусть выполнены условия (18.3) и
<А‘< 0 (<pi > 0) i — 1, 2,... ,п — 1. Тогда решение задачи (18.1) — (18.2),
отличное от постоянной, не может принимать наибольшее положительное
(наименьшее отрицательное) значение во внутренних точках.
Доказательство. Пусть pi < 0. i = 1,2,..., п — 1. Предположим про-
тивное, что в некоторой внутренней точке i — р выполено ур = max yi =
М > 0. Так как yi congt, найдется внутренный узел к, в котором до-
стигает своего положительного максимума у к = М > 0. Тогда найдется
соседний узел, например к — 1, в котором yk-i < М. Может случиться,
что в нескольких узлах, соседних с узлом к, yi = М. Тогда в качестве
узла к возьмем меньший по номеру и запишем (18.1) в узле к :
ЬкУк+i — СкУк + dkyk-i = bk(yk+i — Ук) — dk($k — Ук-i) — (ск — bk — dk)yk < 0.
Но < 0, i=l, 2,..., п — 1, т. е. это выражение должно быть неотрицатель-
но. Получено противоречие. Второе утверждение доказывается аналогично.
Теорема доказана.
Следствием доказанной теоремы являются следующие утверждения.
Теорема 2. Пусть выполнены условия (18.3) и > 0. i — 1,2,..., п — 1
и Д1>0, ^2 > 0- Тогда yi > 0, i = 1, 2,..., п — 1. Если же
pi < 0, i—1,2,..., п—1 и /Xi < 0, /12 < 0. Тогда yi < 0, i = 1,2,..., п—1.
Доказательство. Пусть pi > 0, i — 1, 2,..., п — 1. и yi < 0 хотя бы в
одном внутреннем узле. Тогда yi достигает наименьшего отрицательного
значения во внутреннем узле, что невозможно.
Теорема единственности. Если выполнены условия (18.3), то краевая за-
дача
97
^Уг+1 ^гУг 4“	1 — 0	— 1, 2, . . . , П 1,
7/0 = 0, у = О,
имеет только тривиальное решение yi — 0 г = 0,1,...,п.
Доказательство. По теореме 2 из <Pi < 0, Mi < 0, м2 < 0 следует yi < 0 .
Аналогично, из >0, //1 > 0, м2 > 0 следует yi > 0. Следовательно,
yi = 0. Что касается задачи (18.1) - (18.2), то решение этой задачи един-
ственно, поскольку разность двух решений удовлетворяет однородному урав-
нению, которое имеет только тривиальное решение.
Теорема сравнения. Пусть yi — решение задачи (18.1) - (18.2), a yi —
решение задачи
lyi = -& i - 1,2,... ,п - 1, у0 = Р1, Уп — У>2>
коэффициенты которой удовлетворяют условиям (18.3) и выполнены усло-
вия
Ы < <А, г = 1,2,... ,п — 1, |mi|<Mi, Ы < М2-
Тогда \yi\ < |&|	i = 0,1, 2,..., п.
Доказательство. По теореме 2 \гр\ > 0. Заметим далее, что для
u^yi-yi и v =	+ yi
получаем задачу (18.1) - (18.2) с правыми частями соответственно (рц — рц >
0 и (pi + <pi > 0 и условиями на границе соответственно Ml — Ml > 0 и
Р>2 + М2 > 0- По теореме 2 щ > 0 и Vi > 0 , т. е.
Уг - Уг > 0,	+ yi > 0
и
~Уг <yi< Уг-
Теорема доказана.
Оценка решения первой краевой задачи
Ниже в качестве нормы сеточной функции будет использована первая
векторная норма
Iloilo = max |&|.
98
Теорема. Пусть выполнены условия (18.3), причем
9i = |с»| ~ IM ~ \di\ > 0 i = 1,2,	1.
Тогда справедлива оценка
НЛ < ||^|| +max(H,H).
Доказательство. Решение задачи (18.1) - (18.2) представим в виде у^ +
у(2\ где у^ — решение однородного уравнения с заданными граничными
условиями,
М+1 -	+ diVi-i = 0 г = 1,2,..., n — 1,
= Mb Уп} =
(2)
а у- — решение неоднородного уравнения с нулевыми граничными усло-
виями
М+1 - Ciy[2} + diyf\ = -<&	i = l,2,...,n — 1,	(18.4)
„(2) _ n	?.(2) _ q
У6	u,	yn u.
Покажем, что для y^ справедлива оценка
||?/(1)||с < max(|/zi|, |д2|).
В самом деле, для решения задачи
biyi+i - СгУг + diyi-x = 0 i = 1,2,..., п - 1,
£о = Д, Уп = &, Д = max(|/ii|, |д2|)
замечаем, что Уг > 0 и может достигать максимума только на границе, т. е.
при i = 0 или i = п. Но по теореме сравнения < |у,|.
Займемся оценкой г/г-2\ Пусть максимум модуля yi достигается в узле
г = р :
.Уг\ < |з/р|, i = 0,1,..., п. Перепишем (18.4) в виде
CjJ/i = ЬгУг^-1 + diyi-i +
В этом узле имеем
|ср| |?/р| = |^p?/p+i Н- dpyp-i + <£>| < (|&р| Т |^р|) |?/р| 4" l^pl-
99
Отсюда
w S 7T7 i'_ 'iPiLirf n S Й1-
Ucpl гя l“pl) 9
Теорема доказана.
Обсудим теперь вопросы аппроксимации, устойчивости и сходимости при
решении сеточным методом первой краевой задачи для дифференциального
уравнения
й + p(t)u 4- q(t)u = f(t)
w(0) = ^i,	и(1) = д2-
Как показано выше, этой непрерывной задаче можно поставить в соответ-
ствие сеточную задачу
ЬгУг+1 - СгУг + Ф^-1 = ~¥>i	2 = 1, 2. . . . , П - 1
УО = Ml Уп = М2,
где
bi — 1 + Pih, Ci = 2 + pih- qih2, di = 1 <Мг = fah2.
При p(t) > c > 0 условия (18.3) могут быть выполнены при любых огра-
ниченных q(t) за счет уменьшения h. При p(f) = 0 и q(t) > 0 условия
(18.3) не выполнены. Как известно, в этом случае и непрерывная задача
также может не иметь решения. Случай p(t) < 0 требует дополнительного
исследования, однако в приложениях он обычно не реализуется, поскольку
физический смысл этого коэффициента — трение, которое в приложениях
неотрицательно.
Подведем итог. Мы показали, что задача (16.4)- (16.5) аппроксимирует
задачу (16.1)- (16.2) с порядком 0(h). При выполнении условий (18.3) первая
краевая задача устойчива и, следовательно, решение первой краевой задачи
для сеточного уравнения (16.4) сходится к решению первой краевой задачи
уравнения (16.1).
100
§19. Разностные уравнения
Определение. Линейное уравнение относительно функции дискретного
аргумента (сеточной функции) yi = у(г) вида
®Ог?/г 4“ О,1гУг+1 4" • • • 4"ОчтУИ-т — fi О? il, i2, . . .	(19.1)
называется разностным уравнением т - го порядка. Если коэффициенты
ао, ai,..., ат не зависят от i, то уравнение называется линейным разностным
уравнением с постоянными коэффициентами. Если все fi = 0, то разностное
уравнение называется однородным. Рассмотрим разностные уравнения 1-го
и 2-го порядка.
Разностные уравнения и неравенства 1-го порядка.
Рассмотрим разностное уравнение 1-го порядка
4“ ^г?/г+1 ~ fi	О
ИЛИ
Уг+1 = ЦгУг 4- <&.	(19.2)
Если задано уо, то это уравнение становится реку рентным для определения
Уг. По индукции нетрудно получить формулу
г	г г
Уг+1 = УО П Ук 4- £(П ЧзУРк 4- <pi,
k=0	О s=fc+l
доказательство которой предлагается проделать самостоятельно. В случае
уравнения с постоянными коэффициентами qi = q эта формула имеет вид
г
Уг+1 = qi+Lyo 4- ^2 Ql~k(Pk-	(19.3)
к=0
Неравенства 1-го порядка.
Теорема. Пусть дано разностное уравнение 1-го порядка
щ+1 = qyi + fi q > о	(19.4)
и неравенство первого порядка
Щ+1 < qvi + fi	(19.5)
и пусть uq — Vo- Тогда имеет место оценка
Vi < щ	(19.6)
101
и
г-1
Vi <	+ £ r^fk-	(19.7)
fc=O
Доказательство. Вычтем из неравенства (19.5) равенство (19.4). Получаем
рекурентность
vi+1 - ui+i < q{yi - щ) < q2(vi-i - иг - 1) < ... < дг+1(и0 - v0) = О,
и неравенство (19.6) доказано. Неравенство (19.7) следует из (19.3) •
Разностные уравнения 2-го порядка с постоянными коэффициентами
Рассмотрим разностное уравнение
fa/i+i - cyi + dyi-r = fi г = 1, 2,..., n,
(19.8)
коэффициенты которого b, с, d не зависят от i и b 0, d 7^ 0. В § 17
при изучении метода стрельбы мы выяснили, что для нахождения общего
решения исходного уравнения необходимо найти общее решение однородного
разностного уравнения
byi+i - cyi + dyi-i = 0 г = 1,2...........п,
(19.9)
к которому следует прибавить любое частное решение исходного уравне-
ния. Займемся однородным уравнением (19.9). Последовательность урав-
ненй (19.9) можно рассматривать как систем}' п уравнений относительно
п + 2 неизвестных yo,yi,... ,уп+1- Нетрудно заметить, что ранг этой систе-
мы равен п, поскольку в матрице
—с d ... О
Ь —с d
О b —с d
О
О
0 \
О
О
\ 0	0	0 ... Ь —с d /
есть отличный от нуля минор размерности п. Следовательно, система имеет
2 линейно независимых решения — вектор = (у§\ у^\ . •., Уп^т и
вектор = (у^, у^\ ..., Уп^т и общее решение однородного разностного
уравнения имеет следующий вид
Уг = сху^} + с2г/-2).	(19.10)
102
Эти линейно независимые решения можно отыскивать так, как изложено
в методе стрельбы, но мы сейчас рассмотрим другой прием. Будем искать
ненулевые решения (19.9) в виде
У, = Ч\ 9 о.	(19.11)
Подставляя (19.10) в (19.9), получаем уравнение для определения д, которое
будем называть характеристическим
bq2-cq + d = 0.	(19.12)
Как обычно, рассмотрим три случая:
1.	Корни характеристического многочлена вещественны и различны:
/	т-ч	(1) i	(2) i
gi yt q2. В этом случае имеем два решения у[ — q{ и у> = q2, причем
легко видеть, что эти вектора линейно независимы, поскольку есть минор
2-го порядка, например
/ 1 i \
41 Q2
\ • /
отличный от нуля.
2.	Корни характеристического многочлена кратные: qi = q% = q- В этом
(2)
случае второе линейно независимое решение имеет вид у^ = iql. В самом
деле, подставляя у^ в (19.9) и группируя, получаем
b(i 4- l)gl+1 — cid1 + d(i — 1)дг-1 = qi~1[(bq2 — eg 4- d)i 4- (bq2 — d)].
Выражение в первой скобке равно нулю, поскольку q есть корень характе-
ристического многочлена. Выражение во второй скобке также равно нулю,
поскольку в случае кратного корня q = а поскольку корень кратный,
дискриминант D = с2 — 46 равен нулю и
Линейная независимость двух этих решений у^ — ql и у^ = iql также
легко проверяется.
3.	Корни характеристического многочлена комплексные:
с ± i \/ D _	9	4, ,
91,2 =---, D = с2 - 4bd.
103
Используя представление комплексного числа в виде
512 = ре±г{р = p(cos<£> ± г sine/?), /9==\/т>	— wretg——
V О	с
замечаем, что
5^2 = (cos кр ± i sin кр).
Нетрудно видеть, что если комплекснозначная сеточная функция являет-
ся решением линейного однородного разностного уравнения с вещественны-
ми коэффициентами, то тогда вещественная и мнимая части этой сеточной
функции также являются решением этого уравнения. Поэтому решением
являются функции
(1)	к 7	(2) к •	1
ук = р cos к<р и ук — р sin кр.
Для того чтобы убедиться, что эти решения линейно независимы, достаточ-
но в прямоугольной матрице
/ i °	\
р cos р р sin р
у pkCQskip рк sin кр /
найти минор 2-го порядка, отличный от нуля. Это предлагается проделать
самостоятельно в качестве легкого упражнения.
Пример 1. Найти общее решение уравнения
Ук+i ~Ук + 2ук-1 = ак.
Найдем общее решение однородного уравнения
Ук+1 ~Ук + 2ук-1 = 0.
Корни характеристического уравнения
q2 - q - 2 = О
равны Qi = — 1, 52 = 2 и общее решение однородного уравнения
Ук = ci(~l)fc + с22к.
Для отыскания частного решения неоднородного уравнения предположим,
что а 7^—1 и а^2. В этом случае решение можно найти в виде
104
у к ' = Аак. Подставляя предполагаемое частное решение в исходное урав-
1
нение и сокращая на а , получаем
а2 — а — 2
Если же а = 2, решение будем искать в виде = Ак2к. Подставляя в
исходное уравнение, получаем
А[(к + 1)4 - £2 - 2(& - 1)] = 2 А = 1/3.
Случай а = — 1 рассматривается аналогично. Таким образом, общее реше-
ние исходного уравнения имеет вид
Ук = ci (—1)А; + c22fc + у^.
В этом примере явно просматривается аналогия между методами отыска-
ния общего решения для линейных дифференциальных и линейных разност-
ных уравнений. Эта аналогия не случайна, поскольку разностные уравнения
есть дискретный аналог дифференциальных уравнений. В случае линейно-
го дифференциального уравнения с постоянными коэффициентами реше-
ние однородного уравнения ищется в виде ext . Если же время дискретно,
t = kh, то q = еЛЛ.
Пример 2. Пусть х — произвольное вещественное число, такое, что |ж| <
1. Рассмотрим уравнение
Уп+1 2хуп + Уп—1 = 0-	(19.13)
Будем искать решение этого уравнения, удовлетворяющее условиям уо = 1,
У1 = х. Поскольку |а;| < 1, для дальнейшего удобно обозначить cos <5 =
|оф 5 = arccos#. Характеристическое уравнение имеет вид
q2 — 2 cos Sq + 1 = 0.
Корни этого уравнения
<712 = cos 6 ± i sin 6 = e±lS
и общее решение
Уп = ci cos (n arccos х) + с2 sin (п arccos х).
105
Из условий 2/0 — 1) т/i — х получаем ci = 1, С2 = 0, и решение имеет вид
уп = cos (n arccos х).
В § 6 было доказано, что многочлены Чебышева Тп(х) — cos (n arccosх)
удовлетворяют рекурентному соотношению (6.2). Здесь мы заметили, что
решением разностного уравнения (19.13), которое и есть рекурентное соот-
ношение (6.2), при условиях 2/о = 1, yi = х является многочлен Чебышева.
Тем самым доказана эквивалентность двух определений многочлена Чебы-
шева, о чем говорилось в § 6.
106
§ 20. Собственные значения простейшего разностного оператора
Напомним, что собственными значениями линейного оператора А на-
зываются те значения Л, при которых существуют ненулевые решения
уравнения
Ах = Хх,
а сами эти ненулевые решения называютя собственными функциями, соот-
ветствующими собственному значению Л.
В этом параграфе мы рассмотрим задачу на собственные значения для
разностного оператора
uk+1 - 2ик + «t-i = Хи_ h=l_ А=1]2,	(20.1)
п/	п
uq = 0 ип — 0 •	(20.2)
Эта задача является разностным аналогом задачи на собственные значения
для дифференциального оператора
у"(х) = Ху(х) 0 < х < I
y(Q) = Q у(1) = 0.
Из курса дифференциальных уравнений известно, что собственные значения
этого оператора
__z^2
лт — I I ) >
а система собственных функций em(z) =	образует полную ор-
тонормированную систему в пространстве L2Qr), что позволяет применять
метод разделения переменных при решении задач математической физики.
Займемся собственными значениями разностного оператора. Перепишем
уравнение (20.1) в виде
Xh2
uk+i - 2сик + ик-1 - 0, с(Л) = 1 + — •	(20.3)
При |с(Л) | > 1 корни характеристического уравнения
g2 — 2cq 4-1 = 0	(20.4)
вещественны и различны и общее решение имеет вид
ик = cigf 4- ciq^.
107
Из левого граничного условия ио = 0 получаем
ci + с2 = О,
а из правого граничного условия ип = 0 следует
cig? + с2$ = 0.
Отсюда следует, что ci — С2 = 0, т. е. ненулевые решения невозможны при
|с(А) | > 1. При |с(Л)| = 1 корни кратные и общее решение имеет вид
uk = ciqk + c2kqk.
С учетом граничных условий получаем щ = 0, сг = 0, ив этом случае
также возможно только тривиальное решение.
Пусть |с(А) | < 1. В этом случае корни характеристического многочлена
комплексные. Как и раньше, обозначим
/1Ч	Xh2
c(A) = cosa cos о; = 14——.
Тогда корни уравнения (20.4) равны
<712 = cos а ± i sin а
и общее решение имеет вид
Uk = ci cos ка + С2 sin ка.
Из левого граничного условия uq = 0 следует щ = 0, т. е.
Uk = С2 sin ка.
Из второго граничного условия следует, что
С2 sin па = 0
и, следовательно, ненулевые решения возможны в случае, когда
Ш7Г
= •
П
Из (20.5)
—Xh2 = 2(1 — cos а;) = 4 sin2
108
Следовательно,
«	4.2
Ъп = ~T2Sm ~7Г-
№	2
Заметим далее, что |с(Л)| = 1 при ш=0 и при т=п и, следовательно,
т = 1,2... ,п — 1. Далее значения А повторяются. Таким образом, дока-
зана следующая теорема.
Теорема 1. Разностный оператор (20.1) - (20.2) имеет п — 1 собственное
значение	о
4 л Т777Г
Am = —— sin2——	т = 1, 2,... ,п - 1	(20.6)
h 2п
или
Ат = ~т sin2 т= 1,2, ...,п — 1,	(20.6х)
ГЪ	£1
которым соответствуют собственные функции
е™ = csink^- т — 1,2,..., п — 1	(20.7)
п
или, поскольку А = | и kh = х^
е™ = с sin т7ГХк т — 1, 2,..., п — 1.	(20.7х)
V
Заметим далее, что |Am| < |Am+i| в силу того, что функция
4 . о nmh
K2Sm IF
•2
монотонно возрастает при 1 < т < п — 1. Наименьшее по модулю собствен-
ное значение равно
.	4.2^^
1 = ”K2Sm 2Г‘
Обозначив
тгН
модуль Ai можно записать с виде
7г2 sin2 £
Элементарно показывается, что функция
sin£
109
и имеем минимум в окрестности £ = |, из чего следует, что Ai > Д. Отно-
сительно Ап_] имеем
|х |	. 27Г(п ~ 1)	4 /1	2
=	- = -(l-cos-) .
Без доказательства отметим, что если в пространстве сеточных функций
ввести скалярное произведение по формуле
п-1
(20.8)
(ui,^) = У j hUjVj,
г=1
то в пространстве с этим скалярным произведением собственные функции
el = csmk----
п
ортогональны. Для того чтобы норма собственных функций в этом про-
странстве была равна единице, константу с полагают равной
ПО
§21. Разностные методы решения кревых задач
для уравнения теплопроводности
Пусть задан стержень длины I и пусть u(x,t) — температура в точке
я (0 < х < /) стержня в момент времени t. Процесс распространения
тепла в одномерном стержне описывается следующим уравнением:
<211’
dt дх \ дх J
Функция /(#,£) связана с плотностью внешних источников тепла. Вели-
чины к — коэффициент теплопроводности, с — теплоемкость единицы
длины, р — плотность в общем случае могут зависить от х и t. В про-
стейшем случае, когда они константы, уравнение теплопроводности имеет
вид
ди >)д2и	.	/Л<
lTt=a +	(2L2)
После замены = |, <i = p-f, получим уравнение (индексы у t и х мы
опускаем)
ди д2и р/ .	/Л,
di = d^ + f(x,f>'	(21,3)
Задача Коши для уравнения теплопроводности ставится следующим обра-
зом. Требуется найти решение уравнения (21.3), удовлетворяющее началь-
ным
u(x,Q) = uq(x)	(21.4)
и граничным условиям
aiu(0, t) +	t) = 7i(<),
Of2w(l, t) + faux (1, t) =	(21.5)
Говорят, что задана первая краевая задача, если Д = О i = 1,2. Если
же Qfj = 0 i = 1, 2, то краевую задачу называют второй краевой задачей.
Если все эти коэффициенты отличны от нуля, то говорят, что задана третья
краевая задача.
Напомним, что в курсе математической физики решение уравнения (21.2)
для первой или второй краевой задачи ищется в виде ряда по собственным
функциям однородной краевой задачи, т. е. решение, вообще говоря, может
быть найдено лишь приближенно. Кроме того, для ряда краевых условий
111
задача нахождения собственных значений приводит к необходимости нахож-
дения корней некоторого нелинейного уравнения. Таким образом, задача
нахождения решения для уравнения теплопроводности методом разделения
переменных в ряде случаев оказывается нереализуемой.
Ниже будут описаны сеточные (разностные) методы решения краевых
задач для уравнений математической физики, которые мы подробно изучим
для уравнения теплопроводности. Пусть требуется найти значение решения
краевой задачи (21.3), (21.4), (21.5) при t = Т, т. е. u(xj,T). Разобьем
отрезок [0,1] точками i = 0,1,..., п с шагом h = а отрезок [О, Т]
точками tj, j = 0,1,..., m с шагом т =
Совокупность точек (a?j,tj) будем называть сеткой, а множество точек
при t = tk , т. е. точки (£&, хо), (^ь #1) • • • (^> будем называть слоем.
Разностная схема с весами для уравнения теплопроводности
При фиксированном t заменим производную по я в точке х = Xi раз-
ностным выражением
ди ~ Uj+i - щ
дх h 1
а при х = хп — выражением
ди Un Ufi—i
дх h
Вторую производную по х при х = Xi i = 1,2,... п — 1 заменим
выражением
д2и ~ ui+i -	+ Щ-i
дх2	h2
В результате вместо задачи (21.3), (21.4), (21.5) получаем систему
дифференциально-разностных уравнений
dyi(t)	- 2yi(t) + ?/i-l(*) , Г/.Ч 7	19	„	1	/91 «А
—=---------------Д2---------+ w) г = l,2,...,n-1	(21.6)
с начальными
г/г(О) = U0(Xi)	(21.7)
и граничными условиями
+	=71(«),	(21.8)
f L
+ 02y^)-v«-^ = 72(()	(21 g)
f L
112
Введем обозначения
Луц =
yi+i(t) - 2^(f) + j/i-i(t)
h2
В этих обозначениях (21.6) имеет вид
^- = byu + fi(ty
Заменим производную по t разностным выражением
dyi(tj) S/ife+i) ~ З/ife) = уГ1 ~ У1.	(21.10)
dt	т	т
Здесь и ниже у? = y(xi,tj). Правую часть в (21.6) заменим на линейную
комбинацию значений на j-м и j+1-м слое, т. е. на значения правой части
При t = tj И t = tj+i
J+1 - J . ,	.	.
' т ‘	+ (1 -	(21.11)
Здесь а — параметр, в зависимости от значений которого получается кон-
кретная разностная схема, а — некоторая правая часть, например f-. Ни-
же мы увидим, что в зависимости от а правая часть <р может принимать и
иные значения.
Вместе с начальными
У° = u0(xi)	(21.12)
и граничными условиями
«1У1 +	(21.13)
IЬ
с^ + 02У1п~^ =-Л	(21.14)
п
получаем задачу, которую называют разностной схемой с весами.
Методы решения разностной схемы с весами
Наиболее часто применяются разностные схемы при а = 0, а = | и
er = 1. При а = 0 схема называется явной, при а = j — симметричной, или
схемой Кранка — Никольсона, при а = 1 — неявной схемой. Любая разност-
ная схема с весами позволяет реализовать алгоритм определения решения
на неизвестном слое с помощью решения, полученного на предыдущем слое.
113
Рассмотрим этот алгоритм для явной схемы. Она имеет вид
Пусть j = 0. Решение yf на нулевом слое задано начальными условиями
(21.12.). Но для любой точки (г = 1, 2,..., п — 1) первого слоя имеем
у} =У? ++yLi)+ rfi i = l,2,...,n- 1.	(21.15)
Далее, зная ?/}, с помощью левого граничного условия (21.13) при j = 1
определим у^ и аналогично определим у„ с помощью правого граничного
условия (21.14). Таким образом, найдены решения во всех точках первого
слоя. Повторяя этот алгоритм, определим решения на втором слое и т. д.
Неявная схема. При а = 1 разностное уравнение (21.11) имеет вид
yi+1 -vi = vix - 2yi+i+
т	Л2	Ji
п — 1.
Ниже будет показано, что более точная аппроксимация получается при tpi =
f-+2. Обозначим через zi = y{+1- В этих обозначениях после приведения
подобных членов получаем
2f+i - (2 + — )zi + Zi-i = у? - f-+2h2 i = 1,2,... ,n - 1.	(21.16)
T	T
Правая часть в этом уравнении определена, поскольку у? — значения на
предыдущем слое. Вместе с граничными условиями
2 Zl Zq J-|-l
<*i*o + Pi—д— = 71 ,
(21.17)
anzo + fe2" ,гп~1=7Г1	(21-18)
f L
приходим к разностной задаче, которая решается методами прогонки или
стрельбы, которые изложены в §17.
Симметричная схема (схема Кранка — Никольсона). При а = | в (21.11)
получаем
У3*1 ~ У^ = 3//+11 ~ 2^+1 + У1-1	У3+1 ~ %у1 + У1-1
г	27г2	2h2	Ji
i = 1,..., п — 1.
114
Пусть, как и в неявной схеме, Z{ — y{+1. Приводя подобные члены,
приходим к уравнению
2h2
2<г4-1	(2 ~|“	)“F — l — pi i — 1, 2, . . . , 72	1,	(21.19)
Т
где
о l2	।
W	- feti - 2г4' + yU- (21.20)
Вспоминая о краевых условиях (21.17), (21.18), замечаем, что и эта разност-
ная задача решается методами стрельбы или прогонки.
115
§22. Аппроксимация, устойчивость и сходимость разностной схемы
с весами для уравнения теплопроводности
Напомним, что для того, чтобы решение разностной задачи сходилось к
решению задачи (21.3) — (21.5), надо, чтобы разностная задача аппрокси-
мировала исходную задачу с порядком <9(hfc, тр), где h и т — шаги сетки
по х и t соответственно, и была устойчива по правой части и начальным
условиям.
Порядок аппроксимации
Ниже для упрощения изложения мы будем рассматривать первую кра-
евую задачу, т. е. будем считать, что краевые условия в (21.5) имеют вид:
w(0,t) = 71(f),
= 72 (i).
Введем обозначения в безындексной форме
j	j+1	* yi+1 ‘tyi Уг—1
У = У3Ц У = у1 г &У =-------------^2-------’
7/+1 —
У, = ' т \ УМ = ^+1 + (1 -
В этих обозначениях разностное уравнение (21.11) имеет вид
у, = Лг/,’> + ^.	(22.1)
Начальные и граничные условия аппроксимируются точно
yf = u0(xi), y3Q = y1(tj'), у^ = ъ^).	(22.2)
Пусть и = и? — точное решение исходной задачи в узлах сетки и пусть
z = у — и. Тогда, подставляя у — z 4- и в (22.1), получаем разностную
задачу для погрешности z
zt — h.z^ 4- 0, z® = 0, Zq = z* = 0, для всех k,	(22.3)
где
0 = Л.и^ + — ut.
Заметим, что имеет место равенство

= аи 4- (1 — а)и =
(22.4)
116
В самом деле,
U + U ,	1.	U+U . 1Ч/Л ч ~	/,	\
—5— + (^ -	= —— + (а - -)(и - и) = аи + (1 - а)и.
Наша цель — выяснить порядок аппроксимации исходной задачи разностной
схемой с весами, а для этого необходимо получить оценку функции ip через
шаги сетки h и т. С этой целью разложим ip в ряд по степеням h и т в
окрестности точки (х{, t = tj + |т).
Пусть v = v(tj + |). Тогда
,	1 dv r2d2v r^d3v 4
v = t,(t> + T) = e + -T_ + __ + __ + O(r),
V
,	_	1 dv
= v(ij) = ® - 2Tat +
t2 d2v d3v
~ 48 dF +
и
Ui+r - 2Ui + Ui-L _ d2u	4
h2 ~ dx2 [
Следовательно,
u + и
2
= u + O(r2),
du
TUt = U-U = T-—
dt
+ О(Л
Поэтому
, A frr)	* /U + U ,	1. x
ip = Ли[ 1 + - ut = A(—---h (cr - -)rut) + <p - ut =
d2u du	1 d2 du 2	2
=	+ V + (<T - 2^dt} + O(h +T)-
Поскольку й есть решение исходной задачи в точке (xi,tj +1/2), то
,d2u	du. _ j
[dx2	dt’~ h
Отсюда видно, что если в качестве <р взять /, то ф будет иметь следующий
порядок:
ф = О(т + h2) при <т = 0,
117
ф — О(т + h2) при а = 1,
ф = О(т2 -f- h2) при а = 1/2.
В заключение заметим, что в случае краевых условий 2-го типа порядок
аппроксимации по h снижается на единицу, поскольку в граничной точке
приходится аппроксимировать 1-ю производную по х конечно-разностным
соотношением с погрешностью порядка h.
Устойчивость разностной схемы с весами
Мы будем исследовать на устойчивость первую краевую задачу с нуле-
выми граничными условиями для уравнения
J+1
*------= аЛу-" + (1 -	i = 1,2,..., п - 1	(22.5)
т
Уо = °	Уп = 0
и заданными начальными условиями
y® = uo(xi).	(22.6)
Решение этой задачи можно представить в виде
yi = Vj+Wi,
где Vj — решение однородного уравнения с исходными краевыми условиями
Vi = <тЛ^+1 + (1 - <т)А^ i = 1,2,..., п - 1	(22.7)
= 0	vn = °
и заданными начальными условиями
v? = u0(xi),	(22.8)
a Wi — решение неоднородного уравнения с исходными граничными
w>+1 _
—— ------- = crAw/+1 + (1 — cr)Aw- 4- pi i = 1, 2,..., n — 1	(22.9)
= 0	vn = 0
118
и нулевыми начальными условиями
W? = 0.
(22.10)
Будем говорить, что разностная схема с весами устойчива, если для
имеет место оценка
llz/ill < Мц'Ы! +
Очевидно, что для получения этой оценки необходимо получить оценки для
Vi и Wi :
1Ы| < Mi||u0(zi)|| и ||wj| < М2||<#||.
Если имеет место оценка для Vj, то говорят об устойчивости по началь-
ным данным, а если имеет место оценка для Wi, то говорят, что разностная
схема устойчива по правой части. Если Mi и М2 не зависят от шагов сет-
ки h и т, то разностную схему называют абсолютно устойчивой. Если же
эти оценки выполняются при определенном соотношении между h и т, то
говорят, что разностная схема условно устойчива.
Для исследования на устойчивость по начальным данным разложим v^
— значение сеточной функции на j-м слое по собственным функциям сеточ-
ного оператора
vi+i - 2vi + Vi-i
&Vi =---------------- v0 = vn = 0.
В §20 было показано, что собственными значениями этого разностного опе-
ратора являются
771
4 . 9 7rmh
K2Sm ~Г
m = 1,2,..., п — 1,
которым соответствуют собственные функции
ет =
. Ш7Г
2 sm г--
п
Эти собственные функции ортогональны в сеточном скалярном произведе-
нии
п— 1
(u, v) = UiVih
i=l
и нормированы на единицу в пространстве с нормой, порожденной этим
скалярным произведением. Пусть
п— 1
v[j> ~
к=1
(22.11)
119
Коэффициенты этого разложения, естественно, зависят от слоя j, что и от-
ражает верхний индекс. Разложение (22.11) подставим в уравнение (22.7).
Приравнивая коэффициенты при одинаковых е%, получаем
l(J+1) _ i/j)
------L = <тХ^+1} + (1 - <7)Xkb^.
T
Или
bl+' = 9^,
где
1 + (1 - a)AfcT
Qfc = —;----7----•
1 - (rrAfc
Таким образом,
6«+1> = ^+1ь‘0)
И
п— 1
Ч0) = Е'Й’Ч-
к=1
Но
п-1
^ = ч(0) = £б1°Ч.
к=1
В силу равенства Парсеваля
ни2 = E(<M’)2-
к=1
Пусть д = тах|дф Тогда
k
п— 1
ll<<92>E<’)2 = 922W-
к=1
Отсюда ясно, что для того, чтобы получить ограниченную оценку v\p для
любого слоя], необходимо, чтобы выполнялось неравенство |(?fc| < 1. Выяс-
ним, при каких а выполняется условие |g| < 1 или — 1 < qk < 1. Неравенство
qk < 1 означает, что
1 + (1 - cr)AfcT < 1
1 - атХк ~
120
или
^кГ
1 - arXk
0.
Напомним, что т > 0, так как это шаг сетки, а А& < 0, что следует из
результатов §20. Поэтому условие qk < 1 означает, что
1 - атХк > 0,
т. е.
1
Таким образом, неравенство qk < 1 выполняется при любом неотрицатель-
ном а. Рассмотрим неравенство % > — 1 или
,	2 - (2а - 1)тА*	„
9 + 1 =----3------2-1 >0.
1	— (УТ Afc
При а > 0 знаменатель дроби положителен и для выполнения неравенства
необходимо, чтобы
2	- (2а - l)rAfc > 0
или
Поскольку
Л 4 . 2 kirh
A‘ = ~/?sln —’
то
> 1 _ h2 > 1 _
“2 4т sin2 — 2 4т
или
а > | - ^ = а0.	(22.12)
Полученное соотношение позволяет сделать вывод об устойчивости по на-
кальным данным рассмотренных выше разностных схем.
1.	Явная схема а = 0. Неравенство (22.12) выполняется, если шаги сетки
удовлетворяют соотношению
< 0.5.	(22.13)
ПТ
Таким образом, явная схема условно устойчива.
Ее решение сходится к решению исходной задачи лишь при выполнении
условия (22.13).
121
2.	Схема Кранка — Никольсона а = Поскольку hur положительны, со-
отношение (22.12) выполняется всегда. Схема Кранка — Никольсона безусловн.
устойчива, и ее решение сходится к решению исходной задачи со скоростью
О(т2 + h2).
3.	Неявная схема а = 1 безусловно устойчива. Решение неявой схемы схо-
дится к решению сходной задачи со скоростью О(т + /г2).
Исследование на устойчивость по правой части проводится по той же
методике. Можно показать, что условие (22.12) является достаточным и для
устойчивости по правой части. Более подробно эта тема изложена в [11].
122
§ 23. Разностные схемы для уравнения колебаний струны
Пусть u(a?, t) — положение точки х однородной струны в момент времени
t. Тогда при малых отклонениях от положения равновесия положение этой
точки х описывается уравнением
д2и •>д2и г. .	_ Л
~д^ = а ~дх?	0 < х < Z, t > О,
или после соответствующих замен
д^и д2и
=	+	0<х<1, 0<t<T.	(23.1)
otz ox1
При этом задана начальная форма струны и начальное распределение ско-
ростей
и(ж,0) = v(x),	du^0) =	(23.2)
С/ L
Граничные условия имеют вид
tz(O,t) = Д1ВД, u(l,t) = ^2(0-	(23.3)
Введем в области D = (0 < х < 1, 0 <t < Г) сетку с с узлами (xi,tj)
х -- ih, h = tj — jr, т = и рассмотрим задачу нахождения решения
для разностной схемы с весами
7?+1 - 9?/ -I- 7?-1
= Л(^+1 + (1 - ia'jy’ + ауТ1) + f ,	(23.4)
т
причем
Л - т/9
?/o = Rfe), Уп = ^2^), y^v^Xi), ------------ = Wi.	(23.5)
т
Напомним, что
д ___ ^i+l 2Zj + Zj—i
' _ Л2 '
Краевые и первое начальное условие аппроксимируются точно, а функцию
w выберем так, чтобы погрешность аппроксимации прозводной по t разност-
ным выражением была минимальной. Поскольку
u(xi, т) = u(xi, 0) + тй(жг-, 0) + 0.5r2u(xj, 0) + О(т3)
123
из уравнения (23.1)
й(ху 0) = «"(ж, 0) 4- f(x, 0),
поэтому
u(xi, т) = u(xiy 0) 4- ru(xi, 0) + 0.5r2(u"(:ri, 0) 4- f(xi, 0)) 4- О(т3)
и
и(Ж,-,г)-п(^,0) =	+ 0 5^^,^.) +	0)) + о(т2у
т
Следовательно, если
Wi = w(xi) 4- 0.5т2(г/'(яг) 4- /(^, 0)),
то погрешность аппроксимации второго начального условия разностным вы-
ражением будет порядка O(h2}.
Таким образом, задаче (23.1) — (23.3) поставлена в соответствие разност-
ная задача (23.4) — (23.5), которая решается, как и разностная задача для
уравнения теплопроводности, методами прогонки или стрельбы. Применяя
те же методы, что и в предыдущем параграфе, можно показать, что раз-
ностная задача (23.4) — (23.5) аппроксимирует исходную задачу (23.1) —
(23.3) с порядком <9(т2 4- /г2).
Исследование на устойчивость разностной схемы (23.4)— (23.5)проводит-
ся так же, как и иследование на устойчивость разностной схемы для урав-
нения теплопроводности. Исследование на устойчивость по начальным дан-
ным приводит к следующим результатам.
При <7 = 0 (явная схема) для устойчивости по начальным данным необ-
ходимо, чтобы шаги сетки удовлетворяли
г! 1
h2 ~ 1 4- £1
т. е. явная схема условно устойчива.
При
. 14-е
а>—-
разностная схема безусловно устойчива.
Условие (23.6) позволяет для 1-й краевой задачи доказать устойчивость
по правой части при а > 0.
соотношению
£ > 0,
(23.6)
124
§ 24. Разностные схемы для уравнения Пуассона
Рассмотрим задачу Дирихле уравнения Пуассона
Aw =
д2и д2и .
fai + d^ = ~f(x’
(241)
Требуется найти решение этого уравнения внутри прямоугольника
Ст — С? U Г — {0 < х < Zj, 0 у < ^2})
принимающего на границе Г заданные значения:
u|r = ^(s).	(24.2)
Для численного решения задачи (24.1) — (24.2) разобьем отрезки [0, Zi]
и [О, I2] на П1 и П2 частей соответственно с шагами hi = 1/ni и h,2 = 1/п2,
определив тем самым сетку w с границей 7. Граничными узлами будем счи-
тать узлы (0, yj), (li, yj), (xi, 0), (xi, l2), i = 0,1,... ni, j = 0,1,... n2, t. e.
узлы на сторонах прямоугольника. Обозначим через Vij сеточную функцию,
заданную на сетке w. Аппроксимируем вторые производные выражением
0 U Uj+ij lUij + Wj—ij
~ ц +
д	 4"	1
Пользуясь этими выражениями, заменим исходное уравнение (24.1) разност-
ным уравнением

Vi,j+1 “Ь Vij—1
До
— fij •
(24.3)
С учетом граничных условий
Vio = v(xi,O), vin2 = n(xi,l2), wOj = д(0,2/;), vnij = j) (24.4)
получаем разностную задачу (24.3) — (24.4). Мы уже оценивали аппрок-
симацию 2-й производной разностным выражением. Поскольку граничные
условия аппроксимируются точно, легко можно получить, что при f =
разностная задача (24.3) — (24.4) аппроксимирует исходную задачу (24.1) —
(24.2) с порядком O(h2 + h2).
125
Если сетка квадратная, hi = h2, то уравнение (24.3) можно переписать в
виде
_	+ Vi+l,j +	+ vi,j+l + h?fij
vij -	4
Отсюда видно, что в случае однородного уравнения значение в каждой точке
есть среднее арифметическое от значений в соседних точках.
Среди методов решения задачи (24.3) - (24.4) разделяют прямые методы
и итерационные методы. Прямые методы применяются, когда легко можно
найти собственные значения и собственные функции некоторой разностной
задачи, что возможно, если исходное уравнение есть уравнение с постоян-
ными коэффициентами и область G есть прямоугольник. Если это не так,
то применяют итерационные методы, поскольку задача (24.3) — (24.4) есть
линейная система (ni — 1)(п2 — 1) уравнений относительно (ni — 1)(п2 — 1)
неизвестных Vjj i = 1, 2,..., щ — 1, j = 1,2,..., п2 — 1.
Прямые методы. Метод разделения переменных
Пусть Wij — сеточная функция, заданная на сетке w и равная нулю на
границе сетки, а вне границы удовлетворяющая разностному уравнению
где
Wi+ij - 2wij + Wi-ij	Wij+i - 2Wij + Wij-i.
) — t-Piji
hl
hl
(24.5) V
—75— при x = li — hi, 0 < у < l2
hr
1 V&y)	К	A .	.7
H--^2— при x = hi, 0 < у < I2
—^2— ПРИ У = *2 - h2, 0 < x < li
, m(^,o)
+ hl
при у = h,2, 0 < x
во всех остальных точках сетки.
Пусть
Л ...	_	‘Z'Wij + Wj-ij
lYxWij —
hl
Л _ w'J+i - 2wii + WiJ-l
li-yWij —
и
h2
126
В этих обозначениях уравнение (24.5) можно переписать в виде
~ tflij)	(24-6)
w|7 = 0.	(24.7)
Решение задачи (24.5) — (24.6) при фиксированном i будем искать в виде
разложения решения по собственным функциям оператора
k-yC-j — .Абу Cq —	— б.
В §20 для собственных значений и собственных функций этого оператора
получены формулы
.	4 . 2 ктгНъ
At=-^slnw	(24'8)
eT = \frs'n~r1' *: = 1,2,... ,«2 - 1.	(24.9)
V i2 h
Пусть
712 ~1
<24л°)
fc=l
Разложим функцию при фиксированном i по собственным функциям
оператора Ху :
П2-1
(24.11)
к=1
Как следует из рассмотрений §20, собственные функции ортогональ-
ны в пространстве со скалярным произведением
П2-1
(ц-, Wi) = 52 W™Mn)
fc=l
и нормированы на единицу в пространстве с нормой, порожденной этим
скалярным произведением. Это позволяет найти коэффициенты (3*.
$ = 52 h^ije(jk} =	52 /l2sin	(24.12)’
Подставим эти разложения в (24.6)
П2-1	719-1	712“1
-£e<‘>Aia<«-£Ate<V> = £^>.
/с2 = 1	&2 = 1	Д?2 = 1
127
+^-i + AtQW +	= 0
Подставим эти разложения в (24.6)
П2 — 1	По — 1	П2 ~ 1
- £	= £ 0^-
fc2 = l	fc2 = l	A?2 = l
Или
П2-1
(Л)
е}
к2=1
Следовательно,
^!±i----* +Хка^ + 0* = О.	(24.13)
При этом
а<‘> = а» = 0.
(k)
Если ввести обозначение z2 =	, то для гг- получаем краевую задачу
для разностного уравнения 2-го порядка
zi+1- (2-Xbh^Zi + z^i + hlff = 0 i = l,2,...,n1-l	(24.14)
z0 = zni = 0,	(24.15)
которая решается методами прогонки или стрельбы, изложенныьи в § 17.
В заключение заметим, что в случае, если область не прямоугольная или
уравнение с переменными коэффициентами, то изложенный метод обычно
неприменим. В этом случае для решения сеточной задачи применяются ите-
рационные методы решения линейных систем, которые изложены в § 3. Схо-
димость итерационных процессов обусловлена тем, что сеточный оператор
для уравнения Пуассона является симметричным и положительно опреде-
ленным в сеточном скалярном произведении. Более подробно эти вопросы
изложены в [9].
128
Лабораторная работа 1
Эта работа посвящена отысканию максимального по модулю собственного
значения матрицы А. Для нахожения максимального по модулю собствен-
ного значения необходимо задать любой ненулевой вектор Х^ и запустить
итерационный процесс (2.1). Далее, в случае сходимости последовательно-
стей (2.4') или (2.5х) или (2.7) — (2.8) вычислить максимальное по модулю
собственное значение и соответствующий ему собственный вектор, норми-
рованный на единицу в 1-й векторной норме.
Задание 1
/ 2 -i 1 \
I 1 2-1
\ 1 —2	2 /
Задание 2	-
/ 2.1	1 1.1 \
I 1 2.6 1.1 I
\ 1.1 1.1 3.1 /
Задание 3
/ 1 -2 -1 \
I -1 1 11
\ 1	0 -1 /
Задание 4
/ -3	4 -2 \
101
\ 6 -6	5 /
Задание 5
/1 -1 -1 \
I 1 1 0 I
\ 3 0 1 /
Задание 6
-3	4 -2 \
10 11
6-6	5 /
Задание 7
2 1	0 \
13-1)
-	12	3 /
Задание 8
-	3	2 2\
-3-11
-	1	2 0/
Задание 9
3 -3 1 \
3 -2 2 |
-1	2 0/
Задание 10
1-2	2 \
1	4 -2 I
1	5 -3 /
129
Задание 11
Задание 18
/ 2.4 1.2 -0.3
1.2 1.9	1.4
\ -0.3 1.4	0.8
Задание 12
Задание 19
2 -1 -1 \
1	0-1
3 -1 —2 /
1.6
1.2
-1.1
1.2
1.1
0.6
-1.1
0.6
0.8
Задание 13
/ ! —2	2 \
1 4-2
\ 1	5	- /
Задание 14
/ 1 —2—1 \
I -1 1 1
\ 1 о -1 /
Задание 15
/2.7	1	1.7	\
|	1 3.2 1.7
\ 1.7 1.7 3.7 /
Задание 20
/ 3.3	1	2.3	\
|	1	3.8	2.3	I
\ 2.3 2.3 4.3 /
Задание 21
Задание 22
2.6
1.2
-0.1
1.2 -0.1
2.1	1.6
1.	0.8
Задание 16
/ 1.4 1.2 -1.3 \
1.2 0.9	0.5
\ -1.3 0.4	0.8	/
Задание 17
/ 2.9	1	1.8 \
I 1 3.3 1.8
\ 1.8 1.8 3.8 /
Задание 23
/ 1.8 1.2 -0.9
I 1.2 1.3 0.8
\ -0.9 0.8	0.8
Задание 24
/ 3.4	1	2.4	\
I 1 3.9 2.4
\ 2.4 2.4 4.9 /
130
Лабораторная работа Jft 2.
Лабораторная работа посвящена итерационным методам решения линей-
ных систем. Для нахождения решения студенту предлагается в зависимости
от варианта реализовать итерационный процесс по методу простой итера-
ции, либо по методу Зейделя, либо по методу скорейшего спуска. В качестве
начального приближения к решению можно взять любой вектор Х^. Ите-
рационный процесс прекращается при выполнении условия
На печать выводится результат — вектор и вектор невязки
=	- Ь.
Кроме того, присдаче работы необходимо владеть материалом, изложенным
в §3.
Метод простой итерации
Задание 1.1 3.6xi +	+ 2.4х.з = 0.26 X] Н- 3.9^2 + 2.4хз = 0.31 2.4xi + 2.4x2 + 4.9хз = 0.76 Задание 1.2 3.1X1 + Х2 + 1.6хз = 2.16 Xi + 3.1x2 + 1.6хз = 2.21 1.6xi + 1.6x2 + 3.9хз = 1.46	Задание 1.3 3.5x1 + Х2 + 2.1хз = 0.56 Xi + 4x2 + 2.5хз = 0.61 2.1X1 + 2.5x2 + 4.7хз = 0.96 Задание 1.4 2xi + 0.9x2 + хз = 2.2 0.9x1 + 1.8x2 + О.Зхз = 2.31 Xi + 0.3x2 + 1.9хз = 2.06
131
Задание 1.5
Задание 1.10
1.3ж 1 4- 0.4хг + 0.5хз = 0.25
0.4a;i 4- 1.3x2 + О.Зхз = 0.51
0.5xi 4- 0.3x2 + 1.3хз = 1.44
3.1X1 + Х2 + 1.8хз = 3.2
xi 4- 3.3x2 4- 1.8хз = 2.1
1.8xi 4- 1.8x2 4- 3.9хз — 1.6
Задание 1.6
Задание 1.11
1.6x1 4- 0.7x2 4- 0.8хз = 1.2
0.7xi + 2.6x2 4- О.Зхз — 1.1
0.8x1 4- О.Зхг 4- 1.6хз = 2.3
3.6xi 4- Х2 4- 2.3хз = 1.26
xi 4- 3.8x2 4- 2.3хз = 0.11
2.3x1 4- 2.3x2 4- 4.7хз = 1.76
Задание 1.7
Задание 1.12
2.6xi 4- Х2 4- 1.2хз = 0.86
xi — 277x2 + 1.2хз = 2.2
1.2x1 4- 1-2x2 + З.2х3 = 1.44
3xi 4- Х2 4- 1.9хз = 3.2
Xi 4- 3.4x2 4- 1.9хз = 2.2
1.9xi 4- 1.9x2 4- 4хз = 0.6
Задание 1.8
Задание 1.13
3xi 4- Х2 4- 1.5хз = 3.2
Xi 4- 3.6x2 — 1.5хз = 1
1.5xi — 1.5x2 4- 3.5хз = 2.6
5.6xi 4- Х2 4- 2.3хз = 1.26
xi 4~ 7.8x2 4” 2.3хз = 0.11
2.3x1 + 2.3x2 4- 6.7х3 = 2.76
Задание 1.9
Задание 1.14
1.4xi + 0.5x2 4- О.бхз = 0.66
0.5xi 4- 1.4x2 4- О.Зхз = 0.37
0.6x1 + 0.3х2 + 1.4х3 = -0.26
4xi 4- Х2 4- 1.9хз = 3.2
xi 4- 4x2 4- 1.9хз = 5.2
1.9xi 4- 1.9x2 + Охз = 3.6
132
Метод Зейделя
Задание 2.1
3.1xi + Х2 + 1.62:3 = 2.16
#1 4- 3.1X2 4- 1-6хз = 2.21
1.6a?i 4- 1.6x2 4- 3.9хз = 1.46
Задание 2.2
2xi + 0.9x2 + хз = 2.2
0.9хх + 1.8х2 + О.Зх3 = 2.31
Xi 4- 0.3x2 4- 1.9хз = 2.06
Задание 2.3
1.6x1 4- 0.7x2 4- 0.8хз = 1-2
0.7xi 4- 2.6x2 4- О.Зхз = 1.1
0.8x1 4" 0.3x2 4- 1.6хз — 2.3
Задание 2.4
3xi 4- Х2 4- 1.5хз = 3.2
Xi 4- 3.6x2 — 1.5хз = 1
1.5xi — 1.5x2 4- 3.5хз = 2.6
Задание 2.5
3.1X1 4- Х2 4- 1.8x3 = 3.2
Xi 4- 3.3x2 4- 1.8хз = 2.1
1.8x1 4- 1.8x2 4- 3.9хз = 1.6
Задание 2.6
3xi 4- Х2 4- 1.9хз — 3.2
Xi 4- 3.4x2 4- 1.9х3 = 2.2
1.9xi 4- 1.9x2 4- 4хз = 0.6
Задание 2.7
4xi 4- Х2 4- 1.9хз = 3.2
Xi 4- 4x2 4- 1.9хз = 5.2
1.9xi 4- 1.9x2 4- 6x3 = 3.6
Задание 2.8
3.5xi 4- Х2 4- 2.4хз = 1.2
Xi 4- 4x2 4- 2.5хз = 2.2
2.4x1 4- 2.5x2 4- 4.5хз = 3.1
Задание 2.9
3xi 4- х2 4- хз = 0.2
Xi 4- 3x2 = 2.2
Xi 4- Зхз = 3.6
Задание 2.10
2xi 4- 0.9x2 4- хз = —0.2
0.9xi 4- 2x2 4- О.Зхз = 1.2
Xi 4- 0.3x2 4- 2хз = 2.6
Задание 2.11
4.9x1 4- 0.7x2 4- хз = —5.2
0.7xi 4- 5.1x2 4- хз = 2.2
Xi 4- Х2 4- 6.1х3 = 3.6
Задание 2.12
5.1X1 4- 2x2 4- хз = 8.2
2xi 4- 6х2 4- Зхз = 5.1
Xi 4- Зх2 4- 5хз = 3.6
133
Метод скорейшего спуска
Задание 3.1
5xi 4- Х2 4- 1.4хз = 1.2
xi 4- 4хг 4- 2.5хз = 2.2
1.4xi 4- 2.5хг + 4.5хз = 3.1
Задание 3.6
5.1X1 4- 2х2 4- хз = 8.2
2xi 4- 6x2 4- Зхз = 5.1
Xi 4- 3x2 4- 5хз = 3.6
Задание 3.2
Задание 3.7
3xi 4- х2 4- хз = 0.2
Xi 4- 3x2 = 2.2
xi 4- Зхз = 3.6
xi 4- х2 = 8.2
xi 4- 2x2 4- 2хз = 5.1
2x2 4- 5хз = 3.6
Задание 3.3
4.9x1 4- 0.7x2 4- х3 = -5.2
0.7x1 4- 5.1x2 4- хз = 2.2
xi 4- Х2 4- 6.1хз = 3.6
Задание 3.8
2xi 4- 2x2 ~ 2хз = 10.2
2xi 4- 5x2 — 4хз = —1.1
—2x1 ~ 4x2 4- 5хз = 3
Задание 3.4
4xi 4- Х2 4- 1.9х3 = 3.2
xi 4- 4x2 4- 1.9х3 = 5.2
1.9xi 4* 1.9x2 4- 6x3 = 3.6
Задание 3.9
5xi — 2x2 — 2х3 = 1.2
— 2X1 4- 6x2 = 1-1
—2x1 4- 4хз = 0.9
Задание 3.5
2xi 4~ 0.9x2 4” хз = —0.2
0.9x1 4- 2х2 4- О.Зхз = 1.2
xi 4- 0.3x2 4- 2хз = 2.6
Задание 3.10
6xi 4- 0.6x2 4- 0.7хз = 5.1
0.6xi 4- 3x2 — О.Зхз = 5.1
0.7xi — 0.3x2 4- 4хз = 2.9
134
Лабораторная работа Xе 3
Лабораторная работа посвящена итерационным методам отыскания ре-
шения нелинейных уравнений или систем вида
/1(х1,ж2,. • • ,пп) = О
/П(^1,ж2,...,пп) =0.
Или в векторной форме
F(x) = 0.
Для запуска итерационного процесса нулевое приближение к решению сту-
дент определяет самостоятельно, используя стандартные методы исследова-
ния функций, графические или иные соображения. Итерационный процесс
прекращается при выполнении условия
На печать выводится результат — вектор и вектор невязки
Кроме того, при сдаче работы необходимо владеть материалом, изложенным
в §4.
Методы нахожения решения для одного уравнения
1.	Используя метод секущих (§ 4, формула (4,4)), найти решение урав-
нения:
1.1	х — sin х = 0.25	1.7	1§0.58ж 4- 0.1 = ж2
1.2	у/х — cos 0.41а? = 0	1.8	ж3 — 6ж — 8 — 0
1.3	Зж — cos х — 1 = 0	1.9	ж2 4- 4 sin ж = 0
1.4	ж1пж — 1.5 = 0	1.10	ж3 + 4ж — 6 = 0
1.5	х 4- In ж = 0.5	1.11	ctgx = 0.5ж
1.6	ж3 — sin ж = 0.25	1.12	ж3 4- Зж2 4- 12ж 4-3 = 0
135
2.	Используя метод Ньютона (§4, формула (4.5)), найти решение урав-
нения:
2.1.	ж3 — Зж2 4- 9ж — 8 = О
2.3.	Зж — cos ж — 1 = 0
2.5.	ж2 + 4 sin ж = 0
2.7.	5 ж3 4- 2ж2 — 15ж — 6 = 0
2.9.	^(0.5ж + 0.1) = ж2
2.2.	у[х — соз(0.387ж) = 0
2.4.	ж2 — cos2 7гж — 1 = 0
2.6.	(ж - 1) 4- О.бе* = 0
2.8.	ж4 — 4ж3 + 6ж2 — 4ж + 1 = 0
2.10. 5ж — ех = 0
Итерационные методы для нелинейных систем
3.	Используя метод Ньютона (§4, формула (4.8)), найти решение систе-
мы:
3.1.	sin (ж + 1) — у = 1 2ж + cos у = 2	3.2. cos (ж — 1) 4- у = 0.5 ж — cos у = 3
3.3.	sin ж 4- 2у = 2 2ж + cos (у — 1) 4- ж = 0.7	3.4. cos ж 4- у = 1.5 ж — cos у = 3
3.5.	sin (ж 4- 0.5) — у = 1 cos у — 2 4- ж = 0	3.6. cos (т/ — 1) 4- ж = 0.8 у — cos ж = 2
3.7.	sin (ж 4- у) — 1.1ж = 0.1 ж2 4- у2 = 1	3.8. + 6 = 1 ж — cos у = 1
3.9.	sin (ж 4-1) — у = 1 2ж 4- cos у = 2	3.10. cos (ж 4-1) 4- у — 1 ех — у = 0
3.11.	sin (ж 4-1) — у = 0 ж2 4- у2 = 4	3.12. cos (ж — 1) 4- у = 0.5 ж — cos у = 3
136
4.	Используя метод скорейшего спуска (на выбор: либо формулы (4.14)
— (4.15), либо формулы (4.16) — (4.17) §4)) найти решение системы:
4.1. у = ех — 1
4ж2 4- у = 4
4.2. у = sinx + 1
у = х2 — 1
4.3.	sin (т + 1) — у — 1
2х + cos у = 2
4.5.	sin х — у = О
х2 4- у2 = 4
4.4.	cos (ж — 1) 4- у = 0.5
х — cos у = 3
4.6.	х2 4- у2 = 16
х — cos у = 3
4 7	£1 _1_ у! = 1
1 •	16	25
(х — З)2 4- у2 = 4
4.8. у = 1.5 4- ж3
х2 4- £ = 1
4.9. sin (х 4-1) — у = 1
х2 4- у2 = 9
4.10. cos (х — 1) 4- у = 0.50
х — cos у = 3
4.11.	у = х2	4.12. у 4-= 1
х = е~У	(х - I)2 4- £ = 1
4.13.	z-t/2 = 0	4.14. ^4-^ = 1
у — sin (х — 1) = 1 х — 1 = еу
137
Лабораторная работа. № 4
В работе предлагается вычислить значения функции
ь
а
в точках tj = с + jr, т = двумя способами:
а)	методом удвоения числа шагов для достижения заданной точности,
используя стандартные квадратурные формулы;
б)	используя составные квадратурные формулы Гаусса (11.8) с тремя уз-
лами при фиксированном п — 4.
Пусть
9j&) =
Тогда задача сводится к вычислению в цикле по j интеграла от функции
gj(x) по заданным квадратурным формулам.
Метод удвоения числа шагов состоит в следующем. Все квадратурные
формулы имеют вид
N
g(x)dx « ^ад&) = SN.
i=l
Разобьем отрезок интегрирования точками х-^...,хк и составим сумму S^.
Удвоим N и составим сумму Если величина
|S2W — Sn\ < е,
то считается, что заданная точность достигнута. В противном случае про-
цесс продолжается.
Варианты заданий
Квадратурные формулы
Обозначим через h = gi = gfa), gi+± = gfa + |), g$ = g(a), gN =
g(b). Начальное разбиение N выбирается произвольно, но при использовании
формул Симпсона начальное N четно.
А.	Квадратурные формулы прямоугольников
л-1
sn(9) = h^gi+i.
i=0
138
Б. Квадратурные формулы трапеций
*$jv(<7) =	0	+ 91 + 92 + ••• + 9N-1)-
В.	Квадратурные формулы Симпсона (N - четно)
SN(g) —	+ 9п + 2^1 + 4<?2 + 2#з + 4<у4 + ... + 2^jy_i).
О
Функция f(x,t).
1.	f(x, t) = sin( - t + 0.001т)
1 + Tz
2.	f(x, t) = cos(i - -2 + 0.001т)
3.	f(x,t) = e"(^+0-00b)
4.	/(t, t) =
5.	/(t, t) = ln( 1 +	+ 0.001т)
1 + Tz
6.	=
1. ж<) = н*+^+0-5)
2 + cos(t2 + 1)
8.	f(x = cos^ + ^2 + 0-5)
2 + sin(T2 + 1)
9.	f(xt) =_____Cos(* + x)__
’	2т2 + cos(t2 + 1)
10.	f(x,t) = sin2T
139
Значения a, b, с, d, m. Во всех вариантах е = 0.001.
Номер	а	b	С	d	ш
1	-1	0	0.5	1.5	20
2	0	2	0	1	20
3	0	1	2	3	25
4	1	2	1.5	2.5	25
5	-1	2	2	3	20
6	0.5	1.5	0	1.5	25
7	-2	2	2	3	25
Результаты следует оформить в виде таблицы с выводом значений функ-
ции вычисленной методом удвоения числа шагов и значением N, при
котором достигнута заданная точность, значениями, полученными с исполь-
зованием квадратур Гаусса.
При сдаче работы необходимо знать оценки погрешности квадратурных
формул и методы их получения, методы построения квадратурных формул
Гаусса, свойства ортогональных многочленов.
140
Лабораторная работа К- 5
В работе предлагается численно исследовать поведение решений системы
дифференциальных уравнений в окрестности особых точек и сделать вывод
о характере особой точки.
Пусть задана система	'
= Q&,y).
Напомним, что особой точкой (яо, Уо) называется такая точка, в которой
Р(%о, у0) = 0, Q(x0, уо) = 0.
Необходимо найти все особые точки системы, численно исследовать пове-
дение траекторий в окрестности особой точки и сделать вывод о том, явля-
ется ли данная точка устойчивым положением равновесия и к какому типу
(узел, седло, фокус, центр) относится эта особая точка. Результат числен-
ного исследования должен бать согласован с выводами, полученными ана-
литическим путем. Напомним, что аналитическое исследование проводится
по следующей схеме.
1.	Система линеаризуется в окрестности исследуемого положения рав-
новесия.
Это означает, что в системе делается замена
и = х - хо, v = y -уо
и в полученной системе отбрасываются все нелинейные члены.
2.	В полученной системе вида
du
— = аи + bv,
at
dv
— = си + dv
dt
находятся собственные значения матрицы
141
3.	Если собственные значения вещественные, различные и одного зна-
ка, то особая точка — узел, устойчивый или неустойчивый в зависимости
от того, отрицательны или положительны оба собственных значения. Если
собственные значения вещественны и разных знаков, то особая точка — сед-
ло. Оно всегда неустойчиво. Если же собственные значения чисто мнимые,
то особая точка будет центром в случае линейной системы. Для нелинейной
системы эта точка может быть как центром, так и фокусом. В случае, когда
собственные значения равны, особая точка — дикретический узел.
Численное исследование поведения интегральных кривых в фазовой плос-
кости XOY можно проводить по следующей схеме. Пусть найдено положение
равновесия (а,Ь). Используя предлагаемый вариантом задания метод Рунге
— Кутта, найдем несколько приближенных решений заданной системы, за-
пуская метод Рунге—Кутта с разными начальными условиями. Так, на-
пример, можно численно просчитать траекторию, задав начальные условия
а?(0) = а + е, 1/(0) = Ь, от 0 до некоторого to, затем проделать то же самое,
но уже с начальных условий #(0) = а — е, 1/(0) = 6, затем с начальных
условий я(0) = а, 1/(0) = b + е, затем х(0) = а, 1/(0) = b — е и так далее.
В вариантах задания задается шаг по t, а величины to и в выбираются по
своему усмотрению.
При сдаче лабораторной работы необходимо:
1.	Исследовать поведение траекторий в окрестности каждой особой точ-
ки. Если их бесконечно много, то в окрестности двух любых точек.
2.	Для каждой особой точки вывести на один и тот же экран масштабиро-
ванные графики не менее чем с четырьмя кривыми (каждая кривая разного
цвета) с обозначением начальной и конечной точки.
3.	Аналитически обосновать поведение решений в окрестности каждой
особой точки.
Варианты заданий
Методы Рунге — Кутта
Методы второго порядка. В методах второго порядка положить h =
0.001.
1.	Формула (15.5)
2.	Формула (15.6)
Методы третьего порядка. В методах третьего порядка положить h =
0.005.
142
3.	Формула (15.7)
4.	Формула (15.8)
Методы четвертого порядка. В этих методах h = 0.01.
1. Формула (15.9)
2. Формула (15.10).
Напомним, что эти формулы написаны для скалярного уравнения. В §15
показывается, как применять эти формулы для систем.
Системы дифференциальных уравнений.
х = х — у
у = х2 + у2 -2
х = (2т — у)(х — 2)
у = ху - 2
•	2
X = X — у
у = х2 - (у- 2)2
х = 2х + у2 — 1
у = 6т - у2 + 1
х = у2 — 4т
у = 4у - 8
х = 4 — 4т — 2у
у = ху
X = 1 — т2 — у2
у = 2х
х = 2 4- у — т2
у = 2т(т - у)
х = ху — 4
у = (х- 4)(т/ - х)
X = 1 — т2 — у2
у = 2ху
i = 2(х - 1)(у- 2)
•	9	9
у = у* -X*
х = (х + у)2 - 1
у = -у2 - х + 1
х = (2т — у)2 — 9
у - 9 - (z - 2у)2
х = (2х — у)2 — 9
у = (х- 2у)2 - 9
х = х2 + у2 — 6т — 8у
у = х(2у - х + 5)
2
X = Я — у
у = (х- у)(х -у - 2)
143
Лабораторная работа № 6
В лабораторной работе предлагается найти приближенное и точное
решение дифференциального уравнения
у" +ру' + qy = f(x)
а < х < b,
удовлетворяющего граничным условиям:
+ М(а) = 7^
Я2у(Ъ) + Ы(Ъ) = 72,
на сетке Xi = а + ih, i = 0,1, ...,n, h = найдя точное решение в
узлах сетки и построив и решив сеточную краевую задачу.
Построение точного решения задачи проводится элементарно, с исполь-
зованием стандартных приемов нахождения решений линейных дифферен-
циальных уравнений. Напомним, что общее решение дифференциального
уравнения второго порядка имеет вид:
Y = С\у\ + С2У2 + Уо,
где у\ и у2 — два линейно независимых решения однородного уравнения, а
?/о — частное решение неоднородного уравнения. Два граничных условия,
если определитель
&1У1(а) + Ы^а), а1у2(а) + fay^a)
О12У1 (а) -I- 02У1 (а), &2У2 (а) + /32У2 (а)
отличен от нуля, дают возможность определить константы С\ и С2 в фор-
муле общего решения. Таким образом можно определить конкретный вид
решения, и его следует вычислить в точках Xi.
Построение сеточной краевой задачи проводят по схеме, изложенной в
§16. При сдаче работы необходимо вывести на экран точное и приближенное
решение в точках Студент должен знать методы решения краевых задач
для линейных разностных уравнений 2-го порядка.
Варианты заданий
144
I. Методы решения краевых задач
1. Метод прямой прогонки.
2. Метод обратной прогонки.
4. Метод встречной прогонки.
5. Метод стрельбы.
II.	Уравнения и краевые условия
1.	у" + Чу' + у = —2х
У(0) = 0 т/(1) + 2т/'(1) = 0 п = 20
2.	у" + у = -2е-&
«/(!) +У(1) = 0	s'(2) = l п — 20
3.	у" + 4т/ + Зу = cos х
2/(1) = 1	т/(2) = 1 п = 20
4.	у" + у = -2х
»(0) + у'(0) = 1	у'£) = 0	« = 20
£
5.	у" + 4у' + Зу = —2а:
у'(1) = 0	у'(2) - у(2) = 1 п = 20
6.	у" - Зу' + 2у = sin(z)
Jz(O) + 2у'(0) = -1	1/(1) = 1 п = 20
7.	у" + у' - 2у = —2
у(~1)	= о у(1) — 1 п = 25
8.	у" + 4у' + Зу — —2х
у(0) = 1	у'(1) = 1	п = 25
145
9.	у" + 2у' - Зу = е4х
у(0) = 0	у'(1) = 1	п = 20
10.	у" - у = -2е“21
3/(0) = 0	1/(2) = 1	п - 20
11.	з/" - Зу' + 4у = -2е21
3/'(0) = 0 1/(1) = 1	п = 20
12.	у" - Зу' -4у = -2е2’
!/'(0) = 0	1/'(1) = 1	п = 20
13.	у” + 2у' — Зу ~ ех
3/(0) = 0	т/'(1) = 1	п = 20
14.	у" - 4у' + 8у = —2
3/'(0) = 1	з/(1) =0	п = 20
15.	у" -9у = -2е3х
3/(0) = 0	1/(1) = 1 п = 20
16.	у" + у = —2 siп(х)
у(0) = 0	т/'ф = 1	71 = 25
17.	у” - Зу' = —2х2
1/(0) = 1	з/(1) = 1	п = 20
18.	у" - Зу' + 4У = 1
3/'(1) = 0	1/(2) =0	71 = 20
146
Лабораторная работа № 7
Лабораторная работа посвящена изучению разностных методов решения
краевых задач для уравнения теплопроводности. Пусть задано уравнение
ди д2и ..	.
di= + Кх’
при начальных условиях
и(ж,0) = ф(х)
и граничных условиях
oiu(0,t) +ддц(°>*) = 71(£),
их
a2u(l, t) +	= 72^)-
ох
Требуется:
-	построить сеточную краевую задачу в соответствии с предложенной
разностной схемой,
-	вычислить при t = fo в точках хг-, i = 0,1, ...п приближенное решение,
полученное решением построенной сеточной краевой задачи. Для сравнения
результатов задается точное решение краевой задачи. Способы построения
разностных схем и методы их решения изложены в §21.
При сдаче лабораторной работы студент должен быть знаком с основ-
ными понятиями теории разностных схем (аппроксимация, устойчивость,
сходимость), уметь строить разностную схему для конкретной краевой за-
дачи.
Варианты заданий
I.	Типы разностных схем
1.	Неявная схема.
2.	Схема Кранка — Никольсона.
II.	Методы решения сеточной задачи
1.	Метод прямой прогонки.
2.	Метод обратной прогонки.
3.	Метод стрельбы.
147
III.	Значения n, m, to-
n	m	to
20	25	1
15	25	1
15	15	0.5
15	25	1.5
25	25	1.5
15	25	0.5
IV.	Уравнение, начальные и граничные условия.
ч ди д2и t	2.
1. -уг- =	+ е costtxII + тг )
dt дх1
и(х, 0) = COS7VX ^(О,1) — *4(1> t) = 0-
Точное решение u(x,t) = е* cos тле.
2.
ди д2и
~dt=d^ + X
— 2t
и(х, 0) = 0
Точное решение и(х, t) = x2t.
ди _ д2и д
dt дх2
и(х, 0) = simrx
u(Q,t) = 0 ufx(l,t) = 2t.
; _ I
- Г '~ l
u(0, £) — u(l, t) = 0.
Точное решение u(x,t) — e 7:21 simrx.
du _ d2u
dt dx2
u(x,0) = l-z u'x(^f) = t-l u(l,<) = 0.
Точное решение u(x,t) = (x — 1)(£ — 1).
_ du d2u 2
5-	di = d^~x +x+2t
u(z,0) = 0	u(0, t) = «(1, t) = 0.
Точное решение u(x, f) = t(x — x2).
148
„ du d2u _t .	( 2
6.	— =	4- e г«гп7гж(7г - 1)
dt dx2
u(a;,0) = sinivx	u(0,t) = u(l,t) = 0.
Точное решение u(x,t) = е~*з1п7г.х
_ du d2u 2 r. •
7.	— = 7—7 + x cost — 2smt
dt dx2
u(x,0') = Q u'x(Q,t) = 0	= sinf.
Точное решение u(x,t) = a?2sinZ.
л du d2u t
8. — = — - e-t
dt dx2
c\\ _ ~2
u'x(Q,t) = u'x(l,t) = 0.
Точное решение u(x,t) = х2е~*.
du _ d2u
' ~di = d^
= 0	— t.
Точное решение u(x,t) = tx2.
, _ du d2u
dt dx2
u(x, 0) = e-i
Точное решение u(x, t) = е*~х
u(0, £) = e* u(l,f) = et 1
—х
du _ d2u
’ ~di~d^
u(z,0) = 0	u(Q,t) = t w(l,t)=ie-1.
Точное решение u(x, t) = te~x.
du d2u л/ 9 л .
• 12' л=Г^ + 22: “2t
dt dx2
u(x, 0) = 0	u(0, i) = 0
149
Точное решение u(x,t) = 2trr2.
ди _ д2и
dt дх2
4- х3 — f-x — (i 4- 1)(6ж — 3)
и(х, 0) — х3 — |я2
w(0, /) = 0	u^.(l,i) = 0.
Точное решение u(x,t) = (t 4- 1)(я3 — |z2).
14.
ди д2и
~д1 = ~д^ + Х
и(х, 0) =	u(0, t) = 0
u(l, i) = t.
Точное решение и(х, t) = tx.
du d2u
15‘ di = d^
u(x, 0) = e~x
Точное решение u(x,t) = el~x.
w(0,t) = e<	= ef \
150
Лабораторная работа К2 8
Лабораторная работа посвящена применению метода разделения пере-
менных для прямоугольных областей, изложенному в §24. Пусть задано
уравнение
д2и д2и .	.
Требуется найти приближенное решение этого уравнения внутри прямо-
угольника
Ст — Ст U Г —	0 < у < 12у,
принимающего на границе Г заданные значения:
u|r = /z(s).
Для решения этой задачи необходимо
1.	Отрезок [0,4 и отрезок [ОД2] разбить с шагом hi = £ и /12 —
соответственно. Для Vij = v(xi,yj) построить сеточную задачу
-f- Vj—ij Vij+i 2vij + Vitj—i _
h2 + h2 ~ ~/ij'
vi0 = A(xj,0), vin2 = n(xi,l2),	= м(О,г/Д, vnij = p.(li,j)
2.	Рассмотреть эквивалентную задачу
Wi+ij - 2wij + Wi-ij Wi.j+i - 2wij + Wij-i
[ h2	+	h2 >
где
w|r = 0,
а функция <p определяется следующим образом

lAh,yj)
h2
при j = 1,2, . . . , 712 — 1
-fiJ + ~hT
..	- f _L
^1,712 — 1 — Л,П2~1 “1	^2
— f_i_ 0)
^i,! — Ji,l 4---^2---
при
при
j — 1? 2,..., n2 1
i — 1,2,..., ni 1
при i — 1, 2,..., ni — 1
151
во всех остальных точках сетки.
3.	При фиксированном i = iQ представить решение на слое «о в виде
разложения по собственным функциям оператора
Д _	—	+ WiOtj-i	/ к
Лу ~	Д2	’	V/
— Wi0,n — О,
т. е. в виде (24.10).
4.	По формуле (24.12) вычислить коэффициенты разложения правой ча-
сти tpij при фиксированном i = io по собственным функциям оператора (*),
что позволяет определить коэффициенты в разностном уравнении (24.13).
Решая это уравнение при нулевых граничных условиях методами §20, мож-
но определить коэффициенты в представлении (24.10). Наконец, в цикле по
i можно определить решение на каждом слое.
В результате работы на экран необходимо выдать значение приближен-
ного решения в узлах сетки и значение точного решения, которое дается для
контроля вычислений.
Варианты заданий
Методы решения сеточной задачи
1.	Метод прямой прогонки.
2.	Метод обратной прогонки.
3.	Метод стрельбы.
Уравнения и граничные условия
В скобках приведено точное решение
Задание 1
Ди = 2у
U|j;=o 0	— 0
«1^=0 = 0	U|y=2 = 2(я2 - х)
( и = у(х2 - х) )
0 < х < 1 hi = 0.1
0 < у < 2	/1.2 = 0.1
Задание 2
Ди = 2(х2 -|- у2) — 4
( и= (1 — ж2)(1 — ?/2) )
152
^[а:=0	1 У	^|х=1 — О
W|y=0 - 1 - X2 U|y=d = О
Задание 3
Au = 2(х2 4- у2)
«|х=о = 0 ^|®=i = У2
«|j/=o = 0 U|v=i = х2
Задание 4
О < х < 1
0<2/< 1
hr = 0.1
h2 = 0.1
( и = х2у2 )
О < х < 1	hi = 0.1
О < у < 1	h2 = 0.1
Au — О
( и = (х - 2)2 - (у - 2)2 )
и\х=о = 4у-у и\х=2 = -(у- 2)
^|з/=о ~~ х 2х U|y—2 — (х 2)
Задание 5
Au = Q(x — 2)(у — 1)
«|®=о = 8(1 - у) ик=2 = О
U|y=Q (2 ж) U|y=j — О
О < х < 2	hi = 0.1
О < у < 2	h2 = 0.1
( и = (х - 2)3(</ - 1) )
О < х < 2 hi = 0.1
О < у < 1 h2 = 0.1
Задание 6
Au = 6(х 4- у)
U\x=0 = У3 U|a;=1 = 1 + у3
U|y=0 = х3 U|y=i - 1 4- х3
( и = х3 4- у
О < х < 1
О < у < 1
hi = 0.1
h2 = 0.1
Задание 7
13
18
^_,У^
4 + 9
-
Uk=0 — g
X4
и\у=о =
_ 1	У^_
«|3=1 - 1 + 9
х*
и\у=1 = 1 4- —
О
hi = 0.1
h2 = 0.1
153
Задание 8
Ди = 2(у — 2)
гЬ|а:=О = 0 14^=1 = У - 2
U|y—q —	2х U|y=2 — О
Задание 9
Ди = 12ху
и\х=0 = О U^ = у3 + у
U|y=o = 0 U|y=i — х3 + X
Задание 10
Ди = 2(х2 + у2 — 2х — у)
^|а:=0 — 0 ^4|х=2 — О
U|y=Q — 0 U|y=j — О
( и = х2(у-2) )
0<z<l	Tii = 0.1
0<?/<2	Д2 = 0.1
и = х3у + ху3) )
О < х < 1	/ц = 0.1
О < у < 1	h2 = 0.1
( и = (х2 - 2х)(у2 - у) )
О < ж < 2	hi = 0.1
О < у < 1	h2 = 0.1
Задание 11
Ди = (2 — 6х)у
^|х=0 — 0	1 — О
U|y=0 = О и\у=1 = X2 - X3
( и = (х2 - х3)у )
О < х < 1 hl = 0.1
О < у < 1 h2 = 0.1
Задание 12
Ди = (бу — 2)х
u\x=Q = 0 u\x=i = у3 -у2
0	^|т/=1 — О
( и = (у3 - у2)х )
О < а; < 1 hi = 0.1
О < у < 1 h2 = 0.1
Задание 13
Ди = 4	( и = (х2 + у2) )
«|а;=0 = У2 U\x=2 = 4 + у2
2	-1	।	2
'44|у=О	14|у=1 — 1 “I” X
О < X < 2
0<т/< 1
hi = 0.1
h2 = 0.1
154
Литература
1.	Бабенко, К.И. Основы численного анализа / К.И. Бабенко. — М.: Наука,
1986. - 774 с.
2.	Бахвалов, Н.С. Численные методы. - 2-е изд./ Н.С.Бахвалов. - М.: На-
ука, 1975. — 632 с.
3.	Березин, И.С. Методы вычислений. - Т. 1/ И.С. Березин, Н.П. Жидков.
- М.: Наука, 1966. — 464 с.
4.	Березин, И.С. Методы вычислений. - Т. 2/ И.С. Березин, Н.П. Жидков.
- М.: Физматгиз, 1962. — 640 с.
5.	Воеводин, В.В. Численные методы алгебры/ В.В. Воеводин. - М.: Нау-
ка, 1966. — 208 с.
6.	Годунов, С.К. Разностные схемы. - 2-е изд. / С.К.Годунов, В.С. Рябень-
кий. — М.: Наука, 1977. — 440 с.
7.	Демидович, Б.П. Основы вычислительной математики. — 4-е изд./ Б.П.
Демидович, И.А. Марон. — М.: Наука, 1970. — 664 с.
8.	Крылов, В.И. Вычислительные методы высшей математики. — Т. 1/
В.И. Крылов, В.В. Бобкова, П.И. Монастырский. — Минск: Высшая
школа, 1972. — 284 с.
9.	Марчук, Г.И. Методы вычислительной математики. — 3-е изд./ Г.И.
Марчук. — М.: Наука, 1989. — 608 с.
10.	Самарский, А.А. Введение в численные методы. — 2-е изд./ А.А. Са-
марский. — М.: Наука, 1987. — 288 с.
И.	Самарский, А.А. Теория разностных схем. — 2-е изд./ А.А. Самарский.
— М.: Наука, 1983. — 600 с.
12.	Тихомиров, В.В. Некоторые вопросы теории приближений/ В.В. Тихо-
миров. - М.: Изд-во МГУ, 1976. — 98 с.
13.	Фадеев, Д.К. Вычислительные методы линейной алгебры. — 2-е изд./
Д.К. Фадеев, В.П. Фадеева.— М.: Физматгиз, 1963. — 734 с.
155
Учебное издание
Матвеев Владимир Николаевич
Методы вычислений
Учебное пособие
Редактор, корректор В.Н. Чулкова
Компьютерная верстка В.Н. Матвеева
Подписано в печать 27.04.2007 г. Формат 60x84/8.
Бумага тип. Усл. печ. л. 18,13. Уч.-изд. л. 9,0.
Тираж 100 экз. Заказ
Оригинал-макет подготовлен
в редакционно-издательском отделе ЯрГУ.
Ярославский государственный университет.
150000 Ярославль, ул. Советская, 14.
Отпечатано
ООО «Ремдер» ЛР ИД № 06151 от 26.10.2001.
г. Ярославль, пр. Октября, 94, оф. 37
тел. (4852) 73-35-03, 58-03-48, факс 58-03-49.