от со
#г с

И. Д. Протасов
ЛЕКЦИИ
ПО
ВЫЧИСЛИТЕЛЬНОЙ
МАТЕМАТИКЕ
Москва
«Гелиос АРВ»
2004

УДК 512.8
ББК 22.14+22.19
П58
Рецензент:
профессор В. И. Миллионщиков
Протасов И. Д.
П58 Лекции по вычислительной математике. Учебное пособие.— М.:
ГслиосАРВ, 2004.—184 с.
ISBN 5-85438-122-2
Целью курса является изложение рациональных вычислительных основ
(рациональных алгоритмов и, возможно, их программ) численных методов,
причем доказательства теорем анализа, высшей алгебры и теории вероятностей
оставлены фундаментальным курсам (однако необходимые для специализации
доказательства «вычислительных» теорем приведены).
Отличием настоящего курса является статистическая интерпретация
методов вычислений, что упрощает их обоснование и способствует усвоению
учебного материала.
Курс разбит на лекции, в которые включены упражнения практических
занятий. Предложенное разбиение по лекциям является попыткой автора наряду
с традиционным делением на главы и параграфы использовать в качестве меры
понятие лекции.
ББК 22.14+22.19
Актор благодарен академику И. С. Бахвалову за внимание к работе, автору
и полезные замечания и Д. И. Протасову за полезное обсуждение 17-й лекции.
Учебное издание
ПРОТАСОВ ИГОРЬ ДМИТРИЕВИЧ
Лекции по вычислительной математике
Заведующая редакцией Т. А. Денисова
Корректор Клитина Ε. Η.
Компьютерная верстка О. Ю. Самариной
Издательство «Гелиос АРВ».
Издательская лицензия ЛР № 066255 от 29.12.1998 г.
107140, ι. Москва, Верхняя Красносельская ул., 16., тел./факс: (095) 264-44-39,
e-mail: info^gelios-arv.ru, адрес в Internet: http://www.gelios-arv.ru
Формат 60x84/16. Печеть офсетная. 11,5 п. л. Тираж 2000 экз.
Заказ № Ц81
Отпечатано с готовых диапозитивов в ОАО «Чебоксарская типография № 1».
428019, г. Чебоксары, пр. И. Яковлева, 15.
<0 Протасов И. Д., 2004
ISBN 5-85438-122-2 О Оформление. Издательство «Гелиос АРВ», 2004

3
Памяти
Д.Α., Η.П. и В. Д. Протасовых
ЛЕКЦИЯ 1
Вычислительный эксперимент
ГЛАВА 0. Введение. Вычислительный эксперимент
Одной из основных задач прикладной математики является
построение моделей прикладных задач (адекватных опыту) с целью
прогнозирования поведения ее реального аналога в областях
состояний, недоступных или трудно доступных для эксперимента.
Например, возникает необходимость моделирования конфликтной
ситуации на Среднем Востоке для прогноза хода конфликта или
построения модели человека in vivo для изучения действия
экстремальных гравитационных полей на организм человека и др.
При построении математической модели следует помнить:
— модель тем лучше, чем точнее соответствует (говорят
адекватнее) опыту,
— усложнение модели может увеличить сроки разработки
модели и потребовать невозможных ресурсов от ЭВМ.
Аналитик-исследователь уже не ограничивается знанием
традиционных разделов высшей математики: "Анализ", "Высшая
алгебра", "Теория вероятностей и математическая статистика",
"Дифференциальные уравнения" и др. Это связано с тем, что
внедрение ЭВМ дало толчок развитию дискретной математики
(в связи с построением дискретных моделей исследуемых задач),
следовательно привело к развитию численных (машинных)
методов классической математики и тем расширило понятие
"эксперимент" до понятия "машинный эксперимент".
Вычислительный (машинный) эксперимент — это модельное
исследование естественно-научных проблем (задач) средствами и
методами вычислительной математики, включающее этапы:
а) концептуальное (содержательное) описание изучаемой
модели (с учетом ее ближайшего аналога) на основе данных опыта или
заказчика,

4
b) рациональное математическое описание модели,
c) алгоритм решения модели в рациональной форме (занимает
наименьший объем оперативной памяти, имеет наименьшее время
трансляции и счета и др.),
d) расчет на ЭВМ (программирование, отладка, расчет),
e) анализ полученных результатов для принятия решения:
— об уточнение модели и(или) исходных данных,
— о завершении вычислительного эксперимента.
Замечание 0.1. Существует мнение, что расчет на ЭВМ не
требует упрощения аналитических (например арифметических)
выражений модели и алгоритма. Такое мнение ошибочно:
математик обязан наиболее рационально представить
математическую модель задачи, а программист не только должен выбрать
наиболее рациональный алгоритм (по времени трансляции, времени
расчета и объему памяти), но обязан записать его рационально
(по числу и временной простоте используемых операций).
ПРИМЕР 0.1. Исследовать методом вычислительного
эксперимента задачу о движении тела, брошенного с высоты Η
вертикально вниз с начальной скоростью ν$.
РЕШЕНИЕ, а) содержательное описание модели включает:
— анализ литературы по теме с целью поиска опытных данных
и ближайшего аналога реальной исследуемой модели;
— организацию эксперимента или анализ литературных
данных исследуемой модели с целью установления количественных и
качественных характеристик исследуемой модели (в данном
случае анализ опытных данных позволяет установить пути,
проходимые телом в последовательные отрезки времени, а следовательно,
сформировать математическую модель задачи);
b) в соответствии с пунктом а) математическая модель
исследуемой задачи может быть представлена в первом приближении
уравнением ах2 + Ъх + с = О,
c) алгоритм модели в рациональной форме
Ж1|2 = (-Ь ± у/Ъ2-4ас)/{2а),

5
d) расчет на ЭВМ (программирование, отладка, расчет),
причем при отладке обычно используются опытные данные,
e) анализ данных вычислений позволяет, используя
допущение о пренебрежении сопротивлением воздуха, завершить
эксперимент или принять решения об уточнении модели и продолжении
вычислительного эксперимента. ■
Получение на ЭВМ численного результата не снижает вес
аналитической модели (задачи):
— аналитические ("физические") модели в физике, химии,
экономике, военном деле и др. при своем построении используют
адекватную эксперименту феноменологию и первые принципы.
Примером аналитических, "физических", моделей являются,
например уравнения состояния газов Клапейрона-Клаузиуса, Ван-
дер-Ваальса [16]. Построение таких моделей требует больших
ресурсов и времени и поэтому не всегда возможно;
— аналитические (аппроксимационные) модели обычно не
учитывают природу исследуемой модели или процесса, а полученное
дискретное "приближение" используют как аппроксимант, что
может оказаться не адекватным опыту, например при большом
шаге переменных сеточной функции. Аппроксиманту можно
улучшить, если использовать физически достоверное приближение с
неизвестными параметрами (в этом случае задачей аналитика
является достоверное вычисление параметров модели) [15].
Несмотря на отмеченные трудности, построение
аналитической модели (адекватной эксперименту) остается задачей номер
один для аналитика-исследователя, т. к.
а) позволяет исследовать особенности модели (например
критические и особые точки, асимптотику), с которыми может быть
связано выявление нового эффекта или особенностей реальной
системы,
б) допускает экстраполяцию решения с целью прогноза его
поведения в областях недоступных (трудно доступных) для опыта.
В настоящее время вычислительный эксперимент в полном объе-

6
ме невозможен. Это вызвано тем, что получение аналитической
модели задачи не всегда возможно, а иногда не нужно.
Поэтому курс вычислительной (численной) математики (часто
говорят о машинной арифметике) читают в узком смысле, как
курс машинных рациональных алгоритмов моделей классического
анализа, причем справедливо следующее
Определение 0.1. В вычислительной математике
рациональным выражением (в частности арифметическим) из нескольких
подобных назовем то выражение, которое
— содержит наименьшее число арифметических операций и
— занимает наименьший объем оперативной памяти ЭВМ.
ПРИМЕР 0.2. Записать рационально выражение
У = 0,125а2(8Ь-с). (А)
РЕШЕНИЕ. Запишем исходное выражение в рациональной
форме:
Υ = 0,125 · а2(8 · Ъ - с) = а2 (b - |) .
Выражение (А) содержало 5 действий, а рациональное — 4
действия; параметр 0,125 записал в форме целого 8: занимает меньший
объем оперативной памяти. ■
Вычислительная эксперимент обобщает определение 0.1.
Замечание 0.2. Математическое выражение (модели)
записано рационально (по сути в форме оптимального алгоритма), если
ее алгоритм допускает организацию последовательного процесса
вычисления без участия параллельных процессов. Так выражение
f(x) = ах3 + Ьх2 + сх + d = [(ах + Ь)х + с]х + d (0.1)
переписано рационально. Если не существует последовательного
алгоритма математического выражения (модели), то для
рационального представления выражения необходимо минимизировать
параллельные ветки алгоритма (каждое параллельное
вычисление представить последовательным подалгоритмом и т. д.).И

7
Построение оптимальных (машинных) алгоритмов связано с
дискретизацией задачи (переходом от функции f(x) непрерывного
аргумента к функции дискретного аргумента y(xi)), причем так,
что
Jradxy(Xi) = f(x), (0.2)
следовательно с оценкой погрешности полученного результата.
Вычислительными (численными) методами математики
называют алгоритмы вычисления методами машинной арифметики
рациональных моделей высшей математики, причем главной
задачей численных методов является выбор способа дискретной
аппроксимации и оценка качества приближения.
Традиционные разделы вычислительной математики:
— теория погрешностей,
— приближение (аппроксимация, интерполяция) функций,
— численное решение нелинейных уравнений (НУ),
— прямые и итерационые методы решения систем линейных
алгебраических уравнений (СЛАУ),
— численное дифференцирование и интегрирование,
— численное решение задачи Коши и краевой задачи для
обыкновенных дифференциальных уравнений (ОДУ),
— численное решение дифференциальных уравнений в частных
производных.
Приведены примерные темы индивидуальных занятий,
решаемые методом вычислительного эксперимента:
— вычисление определителя матрицы по Лапласу методами
матричной факторизации и блочной матрицы,
— обращение матрицы методами Ь?7-факторизованной,
присоединенной и клеточной матриц,
— итерационные методы решения СЛАУ (с плотной и
разреженной матрицей коэффициентов),
— анализ, сжатие и хранение плотных и разреженных матриц,
— методы сжатия текстовых данных (реферирование текстов),

8
— методы анализа и восстановления (числовых и текстовых)
данных с пропусками,
— анализ и методы решения задач линейного и динамического
программирования,
— игровые методы среднеквадратичного приближения
(игровая модель метода наименьших квадратов, стратегии, критерии
эффективности приближения и подмножества приближений),
— теоретико-игровые методы численного вычисления
определенных интегралов (модель игры, стратегии, критерии
эффективности квадратурных формул и др.),
— анализ и решение многошаговых игр (игры Риверси, шашки,
крестики-нолики 4x4, задача о разорении игрока и др.).

9
ЛЕКЦИЯ 2
Погрешности чисел. Значащие цифры. Округление
ГЛАВА 1. Основные сведения теории погрешностей
Основная задача теории погрешностей — указать область
неопределенности результата и способы его корректного
представления.
ПРЯМАЯ ЗАДА ЧА теории погрешности состоит в
вычислении погрешностей функции по погрешностям ее аргументов.
ОБРАТНАЯ ЗАДАЧА теории погрешности состоит в
вычислении погрешностей аргументов по погрешностям функции.
Основные источники погрешностей модели (задачи):
1) погрешность модели, связанная с содержанием физической
модели (системы или процесса),
2) погрешность, связанная с методом решения математической
модели задачи,
3) погрешность, связанная с обрыванием вычисления
(например вычисление ряда) или наблюдения (например бросание
монеты),
4) начальная погрешность — погрешность задания
коэффициентов, параметров, начальных или граничных условий
математической модели,
5) погрешность округления (например, число не помещается в
разрядную сетку ЭВМ),
6) неустранимая погрешность, связанная с операциями над
числами, заданными с погрешностью.
Итак, первые пять источников погрешностей — устранимые, а
шестой источник создает неустранимую погрешность, и
устранение этой погрешности может быть проведено только через
устранение погрешностей первых пяти источников.

10
1.1. Абсолютная и относительная погрешности числа
Приближенным называют число а, отличающееся от точного
А (обозначают а « А) и заменяющее последнее в расчетах (так
1,41 < у/2 — А < 1,42), причем Δα = А — а — погрешность
(точность) приближенного числа, причем Δα ^ 0. Поэтому
Определение (абсолютной погрешности |Δα| приближенного
числа а) 1.1.1:
|Δα| = \Α-α\ = Δ-+Α = α± |Δα|. (1.1.1)
Следует различать случаи:
— число А известно, тогда в зависимости от выбора α легко
определяется |Δα| = \А — а|;
— число А неизвестно (что бывает чаще), тогда вместо |Δα|
вводят оценку сверху предельной абсолютной погрешностью а:
|Δα| < α, (1.1.2)
причем а определяется способом измерения (вычисления). При
измерении линейкой предельная абсолютная погрешность по
модулю не превосходит ее минимальной единицы измерения (так на
медицинском градуснике написано: точность измерения 0,1 С).
Поэтому предельная абсолютная погрешность а представляет
теоретический интерес, а часто отсутствует вообще.
ПРИМЕР 1.1.1. Вычислить предельную абсолютную
погрешность приближенного числа α = 2,71, заменяющего е.
РЕШЕНИЕ. В примере А не задано, однако, вычисление
предельной абсолютной погрешности α связано с возможностью
выбора из литературы числа А с разным количеством значащих
цифр. Например,
— для А = 2,71828 получаем
|2,71828 - 2,71| = 0,00828 < 0,009 < 0,01< 0,02 и др.

11
(Здесь лучшее значение а = О,009.)
— для А = 2,7182818 получаем
|2,7182818 - 2,71| = 0,0082818 < 0,008283 < 0,009 и др.
(Здесь уже лучшее значение α = 0,008283.) ■
Зная а можно оценить точное число А:
Λ = α±|Δα| ζα±α. (1.1.3)
Несложно показать, что абсолютной погрешности недостаточно
для характеристики точности вычисления (эксперимента). Так,
если заданы числа 1\ = 10,2сле ± 0,1см и \<ι = 50, Зсм ± 0,1см,
то интуитивно очевидно, что второе число (как большее, но
измеренное с той же точностью) измерено точнее, чем первое.
Более сложно сделать вывод о точности измерения чисел
Ιι = 10,1см ±0,1см и ί2 = 30,2сл1±0,2сл.
С целью сравнения точности измерения (вычисления) вводится
Определение (относительной погрешности δ
приближенного числа а) 1.1.2:
SJM. („4)
Очевидно, чем меньше ί, тем выше точность измерения
(вычисления), причем δ выражают также в процентах.
С учетом определения 1.1.2 в последних примерах вторые числа
измерены точнее:
il = 10^"0'01>52=50^M-0·002'
5l = io^"°'01>52 = 3oJw"0·0067·
Учитывая, что А и (или) |Δα| могут быть неизвестны,
определение 1.1.2 можно переписать в виде
А _ ΙΔαΙ ^ А° - α И1М
'-"Щ"* "Η' ( ]

12
где δα называется предельной относительной погрешностью
вычисления (эксперимента). Зная δα = α/|α|, можно оценить
абсолютную предельную погрешность числа:
\Аа\ ζα = \α\·δα, (1.1.6)
а также точное число:
Α^α±α = α±\α\δα. (1.1.7)
ПРИМЕР 1.1.2. Какое равенство точнее: 17/3 = 5,6(6) или
^ = 7,6?
РЕШЕНИЕ. Находим значения сравниваемых выражений с
одинаково большим числом знаков: 17/3 = 5,666(6), ν58 = 7,616.
Вычисляем предельные абсолютные погрешности [-А — <х| ^ а\:
|5,6666 - 5,66| = 0,0066 ^ 0,0067 = аи |7,616 - 7,6| = 0,016 ^ 0,017 = а2.
Вычисляем относительные предельные погрешности δα = а/\а\:
St т °^Ζ я 0,0012 <«-^- 0,0021,
5,66 7,6
поэтому равенство 17/3 = 5,66 вычислено более точно, чем V58 =
= 7,6.1
1.2. Значащие цифры. Правила округления
Для приближенного числа а или числа, у которого в силу
причин отсутствует погрешность вычисления (измерения) с учетом
естественного представление в виде конечной десятичной дроби:
а = ±(αι · 10η + а2 · КГ"1 + ··.· + ат ■ 10n-m+1 + ...), (1.2.1)
вводится понятие значащей цифры.

13
Определение (значащей цифры) 1.2.1. Значащей цифрой
приближенного числа а называют:
— любую цифру, не равную нулю,
— нуль, стоящий между значащими цифрами,
— нуль, являющийся представителем (умышленно!)
сохраненного десятичного разряда.
Для безусловного выделения значащей цифры целесообразно
представление числа в виде конечной десятичной дроби, причем
умышленное сохранение десятичного разряда конечной
десятичной дроби подчеркивает значимость цифры этого разряда.
ПРИМЕР. По определению 1.2.1 при представлении числа в
виде
0,250100 = 2 · 10"1 + 5 · КГ2 + 0 · 10~3 + 1 · 10"4
последние два нуля — цифры, не значащие, и наоборот, при
представлении —
2·105+5·104 + 0·103 + 1·102 + 0·10 + 0·10°+0·10"1 = 250100,0,
все цифры — значащие.
Представление числа конечной десятичной дробью для
подчеркивания значащих цифр не всегда удобно, а для больших чисел
громоздко.
Представление чисел в форме с плавающей запятой (хотя не
единственно) позволяет компактно подчеркнуть значащие цифры.
ПРИМЕР. Представление 2,501000 · 105 = 25,01000 · 104
подчеркивает, что три последние нуля — цифры значащие, а в
представлении числа 0,2501 = 2,501 · 10~3 = 25,01 · 10~2 только цифры
{2,5,0,1} — значащие.
Пусть имеем число а (или А) и необходимо представить
(округлить) его корректно меньшим числом значащих цифр. Напомним
Правила (округления значащих цифр) 1.2.1. Чтобы
округлить число до η значащих цифр, отбрасывают все цифры,
стоящие справа от п, и заменяют их цифрами по правилам:

14
1) если первая из отброшенных цифр меньше 5, то оставшиеся
десятичные знаки сохраняются, например: 25700203 « 25700200,
2) если первая из отброшенных цифр больше 5, то к последней
оставшейся цифре добавляется 1, например: 25700267 « 25700300,
3) если первая из отброшенных цифр равна 5, а среди остальных
отброшенных цифр есть ненулевые, то к последней оставшейся
цифре добавляется 1, например: 2575002 « 2580000,
4) если первая из отброшенных цифр равна 5, а все остальные
отброшенные цифры равны нулю, то действует правило четной
цифры:
a) если последняя оставшаяся цифра четная, то она сохраняется,
например: 256500 « 256000,
b) если последняя оставшаяся цифра нечетная, то она
увеличивается на единицу, например: 257500 « 258000. При округлении
значащих цифр руководствуются правилом минимизации
разности между округляемым и округленным числами:
аъа*: \а - а*\ =Ф> min. (1-2.2)
Замечание 1.2.1. Повторное округление приводит к
нарушению условия (1.2.2), поэтому запрещено.

15
ЛЕКЦИЯ 3
Верные значащие цифры. Погрешности действий
1.3. Правила подсчета цифр. Верные значащие цифры
При массовых вычислениях с приближенными или точными
числами, а также с числами, у которых погрешность отсутствует,
используют
Правила (подсчета цифр) 1.3.1. Суть правил сводится к
следующему:
— промежуточные вычисления следует получать по крайней
мере с одной запасной цифрой по отношению к значащим цифрам
чисел, участвующим в промежуточном вычислении,
— окончательный результат вычисления содержит то
количество значащих цифр, которое имеет исходное число с
наименьшим числом значащих цифр.
ПРИМЕР 1.3.1. Вычислить выражение:
У = 0,125а2(8Ь-с),
где а = 18;Ь = 2,75;с = 3,232.
РЕШЕНИЕ. Так как погрешность чисел а, 6, с отсутствует,
то вычисления производим в соответствии с правилом подсчета
цифр.
Преобразуем исходное выражение к рациональному виду
(пример 1.3.1):
Υ = 0,125 · а2(8 · 6 - с) = а2 (ъ - |) .
Исходное выражение содержало 5 действий, а окончательное
выражение содержит 4 действия.
Далее последовательно производим необходимые вычисления и
(в соответствии с числом а = 18, у которого две значащие цифры)
записываем результат в форме с плавающей запятой:
Υ = 324 · (2,75 - 0,404) = 324 · 2,346 = 760 = 7,6 ■ Ю2.И

16
Для числа, заданого с абсолютной погрешностью Δ,
существует
Определение (верной значащей цифры) 1.3.2. Значащая
цифра числа называется верной, если абсолютная погрешность
этого числа не превосходит 1/2 единицы разряда, соответствующей
этой цифре.
ПРИМЕР 1.3.2. Задано число 3572,121 ±0,02. Определить его
верные значащие цифры.
РЕШЕНИЕ. Согласно определению 1.3.2
\Аа\ =0,02> J-HT3 = 0,0005,
Δ
|Δα| = 0,02> \· ΙΟ"2 = 0,005,
Δ
|Δο| = 0,02< \ · 10"1 =0,05.
Δ
Итак, в числе 3572,121 ±0,02 только первые пять цифр: 3,5,7,2,1
— верные значащие. Поэтому ответ имеет вид: 3572,1 ± 0,02. ■
В заключение параграфа полезно отметить:
— значащие цифры — это цифры числа, полученные в опыте
(при вычислении), когда погрешность числа отсутствует;
— верные (значащие) цифры это цифры, которые у точного
числа А и приближенного числа а совпадают;
— округляют только значащие или верные значащие цифры.
1.4. Погрешность действий над числами с погрешностью
При действиях с числами, имеющими погрешности,
необходимо вычислить результат действия и его погрешности
(абсолютную и относительную). Рассмотрим действия над двумя числами,
заданными с погрешностью, причем полученные результаты рас-
постраняются на действия с любым количеством чисел (прямая
задача теории погрешностей).

17
ТЕОРЕМА (о погрешности алгебраической суммы) 1.4.1.
Абсолютная погрешность алгебраической суммы слагаемых не
превышает суммы абсолютных погрешностей слагаемых.
Доказательство. Проведено для двух слагаемых χ ± \Ах\, у±
±|Лу|, а результат распространяется по индукции на любое число
слагаемых.
Погрешность суммы: Итак г±|Лг| = х±\Ах\+у±\Ау\, причем
ζ = χ + у, ±|Δζ| = ±\Ах\ ± \Ау\.
С учетом правила треугольника, Va,6 £ Μ : |a + b|<|a| + |6|,
запишем
| ± \Az\\ = | ± |Δχ|| + | ± |Лу|| =» \Αζ\ < \Ах\ + \Ау\.
Далее очевидно αζ ^ \Αζ\ ^ \Ах\ + \Ау\ ^ ах + ау, что допускает
представление αζ = ах + ау.
Погрешность разности: ζ± \Αζ\ = х± \Ах\ — (у± |Ау|), причем
z = x-y,. ±\Az\ = ±\Ax\T\Ay\.
Тогда с учетом правила треугольника можно записать
| ± \Az\\ = | ± |Дж|| + | Т |Ау|| =» \Αζ\ < \Ах\ + \Ау\ =ах + ау. Ш
Обобщая теорему, запишем:
α(Ι±ΙΙ) = αΙ + απ, δ(Ι±ΙΙ) = Ι'δΐΙ+±ΙΙΙΙδπ.
Далее приведем
а(1 ■//) = /· ΙΙ{δ! + δπ) = Π-αΙ+Ι· απ, δ(Ι ■ II) = <5/ + δη,

18
S(Im) = mSr, a(Im) = m-Im-laI.
ПРИМЕР 1.4.1. Вычислить при η = 3,0567 ± 0,0001, т = 5,72±
±0,02:
(п — 1)(т + п)
Υ =
-«\2
(т — п)
и определить погрешности результата.
РЕШЕНИЕ. С учетом теоремы 1.4.1 очевидно:
п- 1 = 2,0567 ±0,0001,
т + η = 8,7767 ±0,0201,
т-п = 2,6633 ±0,0201,
2,66332 7,0932
Последний результат получен с помощью
— правила подсчета цифр (число т содержит три значащие
цифры);
— последовательного округления.
Далее вычисляем
0,0001 0,0201 0,0201 Л ЛЛЛЛ,Л „ ^^~ Л Л ™„^„ ~ „,**
^ = ^+8^+2^=0·000049 + 0·002290 + 2·°·007547:=0·0174·
Следовательно, ау = 2,54 · 0,0174 « 0,0442.
Далее вычисляем верные значащие цифры результата:
0,0442 > i. ΙΟ""2 = 0,005,
0,0442 < i. 10-1 = 0,05.
Итак, цифры 2, 5 — верные значащие цифры, что позволяет
записать ответ: Υ = 2,5 ± 0,0442, 8У = 0,0174. ■

19
ЛЕКЦИЯ 4
Матрица: типы, действия, характеристики
ГЛАВА 2. Алгебра матриц
2.1. Матрицы: определение, типы, действия [8]
Матрица — линейное отображение векторных пространств, в
каждом из которых фиксирован базис; — прямоугольный массив
скаляров размером га х п, заключенный в скобки (круглые,
квадратные) и подчиняющийся правилам алгебры матриц:
ан ... ain
Атхп — {Л)тХп — [CLijJmxn ~ I : ац '. I / . ι ч
η Ι 0 = 1'η)·
(2.1.1)
Далее везде вектор χ = Xmxi — одностолбцовая матрица.
AnXn = Ап — квадратная матрица порядка п. Если Va^ = 0 : А =
= (0) — нулевая, и Vaij = 1 : -А = (1) — единичная
матрицы. Матрицу (θμ) = АТ называют транспонированной, причем
Ат = А — эрмитова (симметричная на R) матрица; —Ат = А —
кососимметричная матрица. Матрицу diagA = diag(an,... ,ann)
называют диагональной, если для \fi ^ j : сщ = 0, причем diagA —
Ε — единичная диагональная, если Мац = 1. Матрицу
называют верхней, правой, U (нижней, левой, L) треугольной, если для
Vi < J -* a».,· = 0 (Vi > j -> a^· = 0):
/flu «12 ain\
0 ... 0 an ... ain
\ 0 0 ann J
= Я = 17,

20
L =
\
N^nl απ2··· α(η-1)η ann'
/an
an
0 .
ац
0
0
0
(2.1.2)
Сложение, вычитание матриц: Χ = А ±В (или хц = ац ± Ьц),
возможно при одинаковых размерах матриц: η а — пв^тл = ^в-
Нулевая матрица (0) — "нуль" относительно операции сложения
(вычитания): А а: 0 = А.
Умножение матриц: А · В (или хц = ац* · bfcj), возможно,
если па = тв- Матрица Ε — "единица" относительно операции
умножения, А · Ε = Ε А = А.
ПРИМЕР 2.1.1. Умножение матриц, согласованных по размеру
2-1 + 3-2 + 4-1 2-0 + 3-2 + 4-4\ /12 22
Ы + 2-2 + 0-1 1.0 + 2-2 + 0-4/V5 4
Деление матриц: X = ^ = АВ~1 с обратной матрицей
рассматриваем как аналог деления матриц (обращение и возведение
А в степень в разделе 2.4).
Матрицы с действительными компонентами (коэффициентами)
сравнимы на множестве бинарных операций: {<,=,>,^,^}· Так,
например, А ^ В, если для Vi,j : ац ^ Ьц.
2.2. Характеристики матрицы
Существуют попытки охарактеризовать матрицу скаляром.
Скалярными характеристиками матрицы являются
определитель, норма, числа обусловленности. Определение нормы
матрицы вводится по аналогии с нормой вектора. Как известно, норма

21
вектора — это длина (модуль) вектора. Ограничимся нормами:
η τη ι
ι = max]T |ay|, ||j1||jj = max^Mi \\А\\ш = */Σα«·
j'=l «=1 у »,j
(2.2.1)
Эвклидова ||il||jj/ — аналог нормы вектора. Свойства норм:
1) || (Л) || = 0 <=> (А) = (0) — нулевая матрица,
2) ||a(i4)|| = a||(i4)||, где a — скаляр,
3)||(Λ)(Β)||<||(Λ)||.||(Β)||,
4)\\(Α) + (Β)\\ζ\\(Α)\\ + \\(Β)\\.
Замечание 2.2.1. По существу (2.2.1) является выборкой
нормы матрицы |\(А)\|, в качестве точечной оценки нормы которой
используем статистическое среднее и дисперсию.
Напомним [11], что математическая статистика для оценки
выборки использует статистические характеристики, например:
Σι#ι «а
Mw = (w) = —г——- — среднее (арифметическое), (2.2.2)
Ш = (w - (w))2 = ^ V 1 W/ дисперсия, (2.2.3)
YlWl\n\w-\
S = Μΐηώ = -(Inw) = -^г ·„ ' г| — энтропия, (2.2.4)
P*(wi) = р^ = \щ \/\W\ — вероятность, (2.2.5)
и др., которые называют выборочными характеристикам
(поведением) наблюдаемой в опыте случайной величины w.

22
Рассматривая нормы (2.2.1) в качестве выборки (реализации)
независимой случайной величины \\А\\, полезны статистические
оценки (2.2.2-2.2.5):
\\А\\ = М(\\А\\)±^Щ\Щ, (2.2.6)
где N — количество (выборка) осредняемых норм, причем
среднеквадратичное отклонение ^/£)(||Л||) можно рассматривать в
качестве абсолютной погрешности вычисления \\А\\.
ПРИМЕР 2.2.1. Вычислить нормы и Μ(||Λ||),Β(||Λ||) матрицы
/1 2 4 0\
А = {ац) =6 0 1 2 1.
\3 1 0 2/
РЕШЕНИЕ:
п
||A||j - max У^ |ο^ | = max{7,9,6} = 9,
г г
j = l
τη
\\А\\ц =тахУ^|а^| = max{10,3,5,4} = 10,
J ΤΞί j
\\Α\\ιη = JTAj = νΊ + 4 + 16 + 36 + 1 + 4 + 9 + 1 + 4 = ^ ^ 8>72
γ «,i
Итак, нормы матрицы близки по порядку величины, причем с
учетом (2.2.6) справедливы оценки:
М{\\А\\) = (9 + 10 + 8,72)/3 « 9,24,
В(||А||) = [(9 - 9,24)2 + (10 - 9,24)2 + (8,72 - 9,24)2]/3 « 0,302,
||A||=M(||A||)±v^pl = 9,24±0,55.1 (2.2.7)

23
ПРИМЕР 2.2.2. Доказать, если (А)(В) = (С) и {А), {В) φ (0),
то
||(А)||>||(ОЦ/||(В)||.
ДОКАЗАТЕЛЬСТВО. Из свойств нормы следует
||(C7)|| = ||(i4)(B)||<||(A)||.||(B)IHII(i4)||>
Н(С)Н
Н(Я)1Г
Характеристикой матрицы (А)п является определитель.
Определение (определителя матрицы) 2.2.1. Всякая
квадратная матрица (А)п характеризуется определителем
det An = det(A)n =
ап
«nl
«in
(2.2.8)
скалярное значение которого определяется аксиоматически.
Определение (минора матрицы) 2.2.2. Если в АтХп выбрать
к < min(ra, n) строк и столбцов, то элементы, стоящие на
пересечении этих строк и столбцов, образуют А*., определитель которой
называют минором fc-го порядка.
Определение (главного минора (А)п) 2.2.3. Главным
минором (А)п называют любой минор, построенный на непрерывном
фрагменте (или на всей) диагонали (А)п, начиная с ац, например
НД
«21 «22
ац
«21
«31
«12
«22
«32
«13
«23
«33
, . . . ,
«11
«nl
«In
(2.2.9)

24
Например, det(A)n — главный минор η-го порядка. Минор,
полученный вычеркиванием (удалением) г-й строки и j-ro столбца.
(Л)п, обозначают М^\
lvliJ
α(<-1)1
а(г+1)1
ai(j-i)
ai(j+i)
α(ί-ΐ)ϋ-ι) a(i-i)(j+i)
α(<+ΐ)0-ΐ) a(<+1)(j+1)
«In
a(i-l)n
a(i+l)n
«nl ··· an(j-l) an(j+l) ··· ann
(2.2.10)
Определение (ранга матрицы, rank(A)) 2.2.4. Рангом
матрицы А называют порядок ее наивысшего минора.
Определение 2.2.5. Определите ль detj42 вычисляют по
правилу
det(A)
2x2
ац αΐ2
а21 а22
«11^22 ~ «12^21,
(2.2.11)
что с учетом [8] позволяет переписать (2.2.11) в виде
2
det(A)2x2 = <
αιιΜί! - αι2Μ£ = ^-Ι^'α^,
JT (2.2.12)
M2^a22 - M^asi = J](-l)2+'a2jM£.
i=i
С учетом (2.2.11) и (2.2.12) справедливо представление (по
Лапласу):

25
det(;4)
ац αΐ2 αΐ3
«21 «22 «23
«31 «32 «33
(-1)1+1αιιΜί1 + (-1)2+1«21Μ2Α1 + (-1)3+1«з1М<
31
= 4£{-l)i+iailM£
г=1
= (-1)1+1«пМЙ + (-l)1+2a12M* + (-l)1+3a1+3Mj
13
(2.2.13)
и далее по индукции
аи ... а 1П
detA
Onl
ап
- ]Γ(-1)·+4ιΜ;? = ^(-lJ^ayMj},
»=1 j=l
(2.2.14)
где A{j = (—lJ'+^My — дополнение (минор со знаком) ау-го
элемента матрицы А.
ПРИМЕР 2.2.3. Вычислить по Лапласу определитель матрицы
А = (ay)
и оценить количество операций при вычислении det(A).
РЕШЕНИЕ. Приведенное подробно вычисление det(A)
|2 4 3
detA = 3 1 -2|
14 11 7
= 2(7 + 22) - 3(28 - 33) + 4(-8 - 3) = 58 + 15 - 44 = 29,
= 2
I -2
II 7
-3
4 3
11 7
4
4 3
1 -2

26
позволяет подсчитать число операций вычисления det (А),
которое целесообразно вычислять по (2.2.13) для каждого минора со
знаком, так для Ац:
— при вычислении Μ и число операций равно 3,
— умножение на αϊ ι,
— сложение 1 + г,
— логическая операция на четность (нечетность) 1 + г,
а число 6 умножить на количество миноров: (3 + 1 + 1 + 1) хЗ =
= 18.И
Свойства определителя:
1) detAT = detA.
2) det{An · Вп) = det{Bn · Ап) = det An · det J3n.
3) Определение 2.2.6. Матрицу называют особенной
(сингулярной, вырожденной), если detAn = 0.
4) Определение 2.2.7. Матрицу А называют строго
невырожденной (неособенной), если все ее главные миноры не равны
нулю.
5) Для положительно определенной матрицы A(Vx φ 0:χΑχτ >
> 0) необходимы и достаточны неравенства Сильвестра всех ее
главных миноров:
ам > 0,
ац αΐ2
α2ΐ а22
>0,
ап
а>п\
air
> 0. (2.2.15)
Наряду с нормой и определителем существуют характеристики
матрицы, связанные с достоверным решением систем линейных
алгебраических уравнений (СЛАУ) и более естественно
названные характеристиками матрицы (коэффициентов) СЛАУ
(например собственные числа и векторы, числа обусловленности),
рассмотренные при решении СЛАУ.

27
ЛЕКЦИЯ 5
Матрицы: разбиение, разложение (факторизация)
2.3. Разбиение матриц
Формирование миноров является частным случаем блочного
разбиения: матрица (А)тхп допускает разбиение на блоки
(подматрицы, клетки), согласованные между собой по размеру, что
часто упрощает матричные вычисления.
Блочное (клеточное) разбиение матриц (А) больших размеров
(порядков) помогает организовать параллельные вычисления и
тем ускорить вычисления, например определителя матрицы Η,
разбитой на четыре подматрицы:
det# =
А В
С D
ζ
= V
1 А В
|(0) D-CA~lB
\A-BD-lC (0)
С D
= \D-CA~lB\detA,
= \A-BD~lC\detD,
(2.3.1)
причем необходимы обратные матрицы A 1,D г (параграф 2.4).
Определение 2.3.1. Матрицу А называют разложимой, если
перестановка строк и (или) столбцов в принципе приводит А к
блочно-треугольному виду, например
\А2
0
А22
(2.3.2)
причем возможно Vi,j : Ац = a>iy
Определение 2.3Л позволяет иметь дело только с
неразложимыми матрицами.
2.4. Матричная факторизация
Разложением матрицы А называют ее представление в виде
произведения нескольких (обычно двух или трех) матриц, что

28
позволяет организовать рациональные преобразования над
разложением сравнительно с преобразованием исходной матрицы А.
В частности разложение строго невырожденной матрицы А
A = LU, (2.4.1)
где L, U (2.1.2) соответственно нижняя и верхняя треугольные
матрицы, в принципе не единственно: число неизвестных п(п +1)
больше числа уравнений п2.
С целью получения единственного разложения (2.4.1) все
диагональные элементы матрицы L\ (или U\) равны единице: в
теории матриц LU-факторизацией строго невырожденной матрицы
А называют ее представление (с точностью до диагональных
элементов) в виде:
А = Шг =LXU,
(2.4.2)
где
/1
О
L(2.1.2) и Ui =
Ul2-
1
^23-
Uln\
У>2п
О
О 1
\0.
(2.4.3)
О 1 /
или
/1 0....
hi 1 0.
ϋι =
0
hi
1 0
о
\lnl 1п2 hn
и 17(2.1.2), (2.4.4)
(n-l)n
1/
а термин "с точностью до диагональных элементов" означает:
на диагонали матрицы U (или L) стоят единицы, что объясняет
Э(!) LU-разложения А. Действительно, матрицы L, U
позволяют записать п2 уравнений: Σ£=1 'tfe^fej = α# Для нахождения η2
неизвестных коэффициентов: /21? · · · > '(n-i)n> wu» ^12, · · ·, wnn.

29
Теорема 2.4.1. Для существования LU-разложения
матрицы А необходимо и достаточно, чтобы А была строго
невырожденной, т.е. все главные миноры матрицы А отличны от
нуля:
«1,1 7^0,
ац αϊ2
«21 «22
7^0,
ац
«nl
«In
φ 0. (2.4.5)
Доказательство. Например, с целью вычисления неизвестных
компонент матриц L и U (2.1.2) очевидны комбинации матриц
главных миноров L,U и А, начиная с det(an) и кончая det(A),
что возможно лишь для ненулевых главных миноров:
(ии) = (ац) -> 1 · иц = ац -» ип = ац ^ 0,
1 0
/21 1
0 1*22
ац ai2
«21 «22
«И «12
a2i а22
/0,
Wll ^12 1X13
0 1*22 ^23
0 0 U33
«И «12
«21
«31
«11 «12
a2i а22
«31 «32
«13
«22 «23
«32 «33 ,
«13
«23
«33
Φ 0 и т. д.
Аналогичное доказательство можно привести для (2.4.4), что и
доказывает теорему. ■
ПРИМЕР 2.4.1. Провести LU-факторизацию матрицы:
Л =

30
РЕШЕНИЕ. Условия существования теоремы 2.4.1 выполнены:
αι,ι=2^0,
= -12 9*0,
4
1
11
Тогда с учетом (2.4.4) можно записать:
/2 4 3
A = LXU= 3 1 -2
\4 11 7
= 29 φ 0.
«π «12 «13
0 П22 «23
0 0 гх3з
Используя правило умножения матриц, последовательно
получаем:
1 · tin = «п =* «п = «п = 2,
1 · и\2 = αΐ2 =Ф «12 = αΐ2 = 4,
1 · «13 = «13 =* «13 = «13 = 3,
fel * «И = «21 =* bl = «2l/«ll = 3/2,
^21 * «12 + «22 — «22 => «22 = «22 - hi * «12 = 1 — (3/2) · 4 = —5,
'21 ' «13 + «23 = «23 =* «23 = «23 - hi ' «13 = -2 - (3/2) · 3 = -6, 5,
fei * «и = «3i => hi = «3i/«n = 4/2 = 2,
«32 - Lm -пи 11-2-4
«22
-5
hi * «12 + ^32 * «22 = «32 => ^32 =
= -0,6,
'31 ' «13 + '32 * «23 = «33 = «33 =* «33 = «33 ~ 'з1 * «13 ~ 'з2«23 =
= 7-2·3-0,6·6,5 = -2,9.
Полученные результаты позволяют записать
1 0 О1
Li = ( 1,5 1 0
2 -0,6 1
и
2 4 3
0 -5 -6,5
0 0 -2,9
detA = det(i · U) = deti · det 17 = 1 · 29 = 29.
ПРОВЕРКА. Полученные матрицы L и U следует проверить:
Ίι«ιι = 12 = ац = 2, /ц«12 = 14 = αι2 = 4, /цг^з = 13 = а13 = 3
и др. ■

31
Наряду с LU-факторизацией матрицы (А) существуют иные,
полезные в вычислительной практике, факторизации матриц:
полярное и сингулярное разложение, скелетное разложение и др.[14].
Полезно дальнейшее разложение LU-факторизации, например
A^LxV^LxDUu (2.4.6)
процедура которого естествена, так как выделяет главную
диагональ U в диагональную матрицу D = diag(dn,... ,dnn). Для
симметрической матрицы А разложение (2.4.6) принимает вид
А = UlDUx = V\V, (2.4.6)
где V = Dll2U\, причем последнее равенство справедливо для
\/ац > 0. В заключение приведем QR-факторизацию матрицы
(А) А = QR с ортогональной (унитарной) матрицей Q и правой
верхней треугольной матрицей R — U~l.

32
ЛЕКЦИЯ 6
Обратная матрица. Обращение матриц
2.5. Обратная матрица. Обращение матриц
Определение Ε позволяет для неособенной матрицы А
определить обратную матрицу А~{:
А А~1=А~1 А = Е, (2.5.1)
причем структуру А~1 обосновывает следующая
ТЕОРЕМА 2.5.1. Всякая неособенная матрица А имеет
обратную матрицу:
Ац А12 ...
A-l _ (а. Λ-1 _ 1 Ι Α*1 ^22 ...
A -[агз> "detA Ац
>Дч An2 ... Ann,
Ац A2i
Au A22 ... An2 \ _ {A)
detA\ Aji ... j detA'
* Ain A2n
(2.5.2)
где (A) = adj(A) — присоединенная матрица, Aij = (—1)г+^ Mij
— дополнение (минор со знаком) а^-го элемента матрицы А.
ТЕОРЕМА 2.5.2. Если матрица А обратима, то А" —
единственна.
Доказательство (от противного). Пусть
ЗВ -+ АВ = В А = Е^
ЗС-+АС = СА = Е
Л:В = ВЕ = В{АС) = {ВА)С = ЕС = С.

33
Теоремы 2.5.1 и 2.5.2 позволяют единственным образом обратить
всякую неособенную матрицу.
ПРИМЕР 2.5.1. Обратить методом присоединенной матрицы
матрицу:
'2 4 3
А = ( 3 1 -2
4 11 7
РЕШЕНИЕ. Вычисляем определитель матрицы:
2 4
3 1
4 11
3
-2
7
= 2
1
11
-2
7
-3
4 3
11 7
+ 4
4
1
3
-2
detA =
Далее вычисляем обратную матрицу:
421
422
^23
АзЛ
А32 =
Азз/
29 5
= -^ I -29 2
29 * 29 -6 -10.
= 29 φ 0.
1 |
29
Л
\
1 -2
11 7
|3 -2
|4 7
3 1
4 111
4
3
11 7
2 3
4 7
2 4
4 1
1
4 3
1 -2
|2 3
|3 -2
2 4
3 1
\
/
ПРОВЕРКА. По определению А~х А = Ε проверяем матрицу
А-1:
(Anan+A2ia2i+A31a3i)/ detA = (29·2+5·3+(-11)·4)/29 = 1 =

34
(Anai2 + A21a22 + A3ia32)/detA = (29-4 + 5-l + (-ll)-ll)/29 =
= О = ei2,
(Aua13+A21a23+A31a33)/det A = (29·3+5·(-2)+(-11)·7)/29 =
= 0 = ei3,
{A12an + A22a21 + Λ32α3ι)/ det A = (-29 · 2 + 2 · 3 + 13 · 4)/29 =
= 0 = e2i,
(Αι2α12 + Α22α22 + Α32α32)/άβίΛ = (-29·4 + 2·1 + 13·11)/29 =
= 1 = e22,
(Ai2a13 + A22a23 + A32a33)/ det A = (29 · 3 + 5 · (-2) +13) · 7)/29 =
= 0 = e23,
(A13an+A23a2i+A33a31)/detA = (29·2+(-6)·3+(-10)·4)/29 =
= 0 = e3i,
(A13a12 + A23a22 + A33a32)/detA = (29 · 4 + (-6) · 1 + (-10) x
xll)/29 = 0 = e32,
(Aisais + A23a23 + A33a33)/ det A = 29 · 3 + (-6) · (-2) + (-10) x
x7)/29 = 1 = e33.
Действительно,
/1 0 0\
E= I 0 1 0.1
Другой распространенный метод обращения матрицы связан с
LU-факторизацией'. А — LU и последующим обращением
произведения треугольных матриц: A~l = (LU)~~l = U~lL~l, причем
последнее равенство подтвержает
Лемма 2.5.1.
{LU)'1 =U~1L-\ (2.5.3)
Доказательство. Домножим 2? = i-1i слева на Ϊ7"1 и справа
на Ϊ7:
υ~λΕυ = U~lU = E = U~lL-lLU = U^L^A.
Если домножить справа на А"1 выражение Ε = U~1L~1A,
получим
ЕА'1 = U-lL~lAA~l -> А-1 = {LU)'1 = [Г1//"1. Ш

35
ПРИМЕР 2.5.2. Используя ££/-факторизацию, обратить
матрицу:
/2 4 3
А = 3 1 -2
\4 11 7
РЕШЕНИЕ. Результаты примера 2.4.1 позволяют записать
1 0 0\ /2 4 3 \
Χι = | 1,5 1 0 , U= 0 -5 -6,5 ,
2 -0,6 1/ \° ° -2'9/
det A = det(L · U) = deti · deti/ = 1 · 29 = 29.
Вычисляем миноры матрицы L~l,U~l методом присоединенной
матрицы:
г I ! °Ι ι. τ - I ° °Ι-η· г -1° °1-п
Ln= -0,6 ι Г ' 21-~ -0,6 ι Г0,131- ι о -°'
Ll2 = -1
1,5 0
2 1
£32 =
-1,5; L22
1 0
2 1
= 1;
1 0
1,5 0
0,
Li3 =
1,5 1
2 -0,6
= —2,9; L23 = —
О
-0,6
= 0,6;
£зз =
1 О
1,5 1
= 1,
ип =
-5
О
-6,5
-2,9
14,5; C/i2 = -
0 -
Ιο -
6,5
2,9
= 0;

36
tf2i = -
4 3
0 -2,9
13
0 -5
0 0
0,
= 11,6; E/22 =
2 3
0 -2,9
= -5,8;
c/23 = -
2 4
0 0
= 0,
U3l =
4
-5
3
-6,5
= -12, Us» = -
2 3
0 -6,5
= 13;
^33 =
1 -1,5 -2.9
= -10,
L~l = f 0
1
0
14,5 0
U~l = U\ И'6 ~5'8 °
zy v -11 13 -10
1 0 0
1,5 1 0
2,9 0,6 1,
14,5 11,6 -11'
0 -5,8 13
0 0 -10
что позволяет вычислить А 1 (аналогично примеру 2.5.1):
14,5 11,6 -11\ /1 0 0Л
А'1 = U-'L-1 =
0
0
29
-29
29
-5,8 13
0 -10/ \
5 —11\
2 13 ) .■
-6 -10/
-1,5 1 0
^-2,9 0,6 1
29
29
Для прямоугольных (А)тхп и особенных матриц определена
обобщенная обратная (псевдообратная) матрица А® = {А~[ ,

37
A^1}, причем Л^1 — левая, А^1 — правая обратные матрицы
[14]:
т>п: AAl1 φ {I) = Al1 A -+AQ= Al1 = (ATA)~l Ατ,
n>m: AA-RX = (/) φ A~KXA -> 4Θ = A~Rl = AT (AAT)~l,
(2.5.4)
где AA°A = A, A°AA° = А°, a AA°, A°A — эрмитовы
матрицы. Для неособенных матриц (А) псевдообратную матрицу
называют обратной, Α~ι = ΛΘ.
ПРИМЕР 2.5.3. Псевдообратить матрицу:
А =
РЕШЕНИЕ. Псевдообращение (2.5.4) позволяют записать
последовательно для случая т> п:
что подтверждает проверка
.::)■(: s:)i(si)-C!)

38
Замечание 2.5.1. С обращением неособенной матрицы А
связана операция возведения А в целую отрицательную степень
А~п = {Апу1 = (д-А--А-...-4[ , (2.5.5)
\ праз }
однако не следует забывать
Определение (идемпотентной матрицы) 2.6.1. Матрицу А
называют идемпотентной, если А · А = А.
Для идемпотентной матрицы понятие возведения в целую
степень (2.5.5) теряет смысл, так как
Ап = А-А-А-...-А = А-+А-п~ А"1. (2.5.6)
η раз
Доказательство справедливости (2.5.6) для идемпотентной
матрицы очевидно и оставлено читателю в качестве упражнения.

39
ЛЕКЦИЯ 7
Характеристики матрицы СЛАУ
ГЛАВА 3. Решение систем линейных алгебраических
уравнений
Классической задачей вычислительной математики является
решение систем линейных алгебраических уравнений (СЛАУ)
АтхпХ = *>, где хТ = (жьЖ2,...,жп),Ьт = (ЪиЪ2,...,Ът).
Возможно, что А — это отображение χ в Ь, А : χ => Ь.
В последние годы интерес к СЛАУ возрос в связи с
интенсивным развитием численных методов математики. Интерес
представляют СЛАУ больших размеров с разреженными матрицами
специальной структуры (например πι φ η). Условно численные
методы решения СЛАУ подразделяют на
— прямые методы (корни в принципе допускают точное
вычисление, например в простых дробях),
— методы итерации (корни СЛАУ вычисляют с заданной
погрешностью ε), причем решению предшествует анализ
совместности СЛАУ (rank A, rank A\b) и устойчивости решения
(обусловленности -А), причем представляет интерес существование 0,1 или
оо решений СЛАУ:
3.1. Анализ СЛАУ (совместность СЛАУ,
обусловленность А)
Совместность СЛАУ устанавливает
Теорема (Кронекера—Копелли) 3.1.1. СЛАУАтхпх = Ь
совместна, если гапкА = гапкА\Ь. Ш
Далее полезен анализ СЛАУ:
Г10ж1 + 9ж2 = Щ

40
Г 10χι+9χ2 = 18,99ϊ
{ 9*1+8*2 = 17,0ΐ)^ = 1'17;*2 = °'81;
10 · 10"6Ж1 + 9 · 10-6ж2 = 19 · КГ6!
й Л ^det(A) = (80-81)10"12 =
9 · 10_6Ж1 4- 8 · 10~6а;2 = 17 · КГ6 I ;
= -10
-12
Существуют СЛАУ, незначительные погрешности вычисления
(измерения) которых, приводят к неустойчивости решения СЛАУ.
Определим скалярную характеристику, отвечающую за
устойчивость решения СЛАУ Ах = 6 при возмущении А и (или) 6:
(А±АА)(х±Ах) = Ъ±АЪ,
Ах + А(±Ах) + х{±АА) + (±АА){±Ах) =Ъ±АЪ,
А(±Ах) + х(±АА) = (±Δ6),
А(±Ах) = -х{±АА) + (±АЬ),
А~1А(±Ах) = (±Ах) = А'^-х^АА) + (±Δ6).
Последнее равенство перепишем в нормах:
Ι|Δχ|| =
= ЦА-'К-х^АА) + (±Δ6)]|| < ΙΙΑ-^ΚΙΝΙΙΙΔΛΙΙ + ||Δ6||),
INI <Ι|Λ Ч" ||+ \Н \\ъ\\
причем с учетом
||6|| = ||Л*||<||Л||.||*||->М<||Л||

41
последнее неравенство принимает окончательный вид
ι|Δχ" < 11Д-Ч
\х\
\\А\\ +т \\Ц\\
и позволяет следующее
Определение (числа обусловленности condA матрицы А)
3.1.1.
condA = \\Α-λ\\ \\Α\\. (3.1.2)
Для качественной оценки cond(A) используем определение
обратной матрицы (2.5.1):
ί1,
[ос,
, _,еслиАА г=А гА = Е,
cond(A) = { (3.1.3)
- если det(A) = 0.
Замечание 3.1.1. Оценка числа обусловленности cond(A)
полезна, так как позволяет в принципе получить представление о
"качестве" вычислительной устойчивости (А) к возможным
погрешностям вычислений, а также необходима при расчетах, так
как предвидит возможные ошибки, связанные с погрешностью
вычисления ЭВМ "некачественной" А, причем вычислительная
практика рекомендует оценки:
1) cond(A) « 1 - матрица (А) "отлично" обусловлен,
2) cond(A) < 10 матрица (А) "хорошо" обусловлена. ■
ПРИМЕР 3.1.1. Вычислить число обусловленности
приведенных матриц:
(л\ (W $\ /т /10-Ю-6 9·10-6\
{А)={9 8У ' (jB) = V 9. Ю-6 8.10-6J·
РЕШЕНИЕ.
det(A) = -1, det(B) = -ΙΟ"12, \\Α\\Τ = \\А\\П = max{17,19} = 19,

42
Τ
-9 10 у V 9 -1
= (-8)(-10)-92 = -1,
^"1 = fii-9 ΤοΊ =17 -;o):det^_1
\\Α~% = \\Α~ι\\π = max{17,19} = 19 : cond{A) =
= \\α-1\\ι>π\\α\\ι<π = μι,
\\Β\\Τ = \\Β\\Π = max{17 · ΙΟ"6,19 · 10~6} = 19 · ΙΟ"6,
ίΒ4-ι_ ΙΟ"6 / 8 -9\Τ_/-8 9 Αιηβ_,,.„-ι
ίο-12 ν-9 ίο у ν 9 ~10,
= ЦВ^Цл = max{17 · ΙΟ6,19 · ΙΟ6} = 19 · ΙΟ6 : cond{A) =
= ||B-1||ilJ/||B||ilff = 361,
следовательно, матрицы (А), (В) плохо обусловлены.
Замечание 3.1.2. Для симметричных А = АТ положительно
определенных матриц:
cond{A) = ^, (3.1.4)
где Amax,Amin — соответственно максимальное и минимальное
собственные числа матрицы (А).
3.2. Собственные числа и собственные векторы
матрицы
Важной характеристикой матрицы являются собственные
числа и собственные векторы матрицы.

43
Определение (собственных чисел λ матрицы А) 3.2.1.
Собственным значением (числом) квадратной матрицы А называется
такое число λ, что для некоторого ненулевого вектора ж,
называемого собственным вектором, выполняется равенство:
Ах = \Ех. (3.2.1)
По сути, нахождение собственных чисел матрицы А сводится к
решению уравнения
det(A - \Е) = 0 (почему)?, так как χ φ 0. (3.2.2)
Задача вычисления собственных чисел и векторов А включает:
— вычисление собственных чисел,
— вычисление собственных векторов по (3.2.1).
Задача вычисления соственных чисел матрицы А достаточно
сложная и решается различными методами, например методом
непосредственного развертывания, методом итерации и др.
Метод непосредственного развертывания: развертывание
характеристического определителя det(A — ХЕ) = 0, т. е.
приведение его к полиномиальному виду (характеристическое
уравнение).
Множество всех собственных чисел матрицы А называют ее
спектром, причем наибольшее по абсолютной величине
собственное число называют спектральным радиусом р(А) матрицы А:
max{|Ai|,|A2|>...,} = p(il). (3.2.3)
Спектральный радиус играет фундаментальную роль в
исследовании сходимости итерационных процессов на матрицах.
ПРИМЕР 3.2.1. Вычислить методом непосредственного
развертывания собственные числа и собственные векторы матрицы:
Ч-< 1)·

44
РЕШЕНИЕ. Итак det А = 11 φ 0. Далее:
1) развертывание
det(A-A£)
3-λ -2
-4 1-λ
= (3-λ)(1-λ)-8 = λ2-4λ-5 = 0,
2) решение характеристического уравнения λ2 — 4λ — 5 = 0
с целью вычисления собственных чисел матрицы А: А^г = 2±
±V4+~5 = 2 ± 3 ==► λι = 5, λ2 = -1,
3) вычисление собственных векторов А:
для собственного числа Αι = 5
з-λ! -2 ν*η_0
-4 1-aJUJ
Решая данную СЛАУ, получаем систему уравнений
-2ал- 2*2=0Ί (1)_ (ц
}=>*ί
. - , «ί/ο —- Х\
—\х\ — 4x2 =0]
Итак, любой вектор вида ( 1/1л I, где х\ φ 0 является
собственным вектором матрицы А. Например, вектор ( 1 I —
собственный.
Для собственного числа Аз = — 1
3-А2 -2 \ /χ?5
-4 ι-λ2Α42)
Решая данную СЛАУ, получаем систему уравнений
4ж(2) - 2Ж(2) - 0 Ί
4а?! 2аг2 -и,^ _^а.(2)=2я.(2)>
-4ж(12)+2х22) = 0)>

45
Итак, любой вектор вида
242)
, где х\' φ 0 является
собственным вектором матрицы А. Например, вектор

собственный.
Далее max{|Ai| = 5, |λ2| = 1,...,} = р(А) = 5. ■
ПРИМЕР 3.2.2. Вычислить методом непосредственного
развертывания собственные числа и собственные векторы матрицы:
А =
РЕШЕНИЕ. Вычисляем по Лапласу det<4 = 324 φ 0. Далее
1) развертывание
det(A - ХЕ) =
11-λ
-6
2
-6
10-λ
-4
2
-4
6-λ
= (11-λ)[(10-λ)(6-λ)-16]+6[(-6(6-λ)+8]+2[24-2(10-λ)+8] =
= (11 - λ)(10 - λ)(6 - λ) - (176 - 16λ) - 168 + 36λ + 8 + 4λ =
= (11 - λ)(10 - λ)(6 - λ) - 56(6 - λ) = (54 - 21λ + λ2)(6 - λ) = 0.
2) решение характеристического уравнения
(54 - 21λ + λ2) (6 - λ) = 324 - 180λ + 27λ2 - λ3 = 0,
1) λι,2 = (21 ± V441 - 216)/2 = (21 ± 15)/2 -> λι =
= 18, λ2 = 3; '
2) (6 - λ) -»■ λ3 = 6.
3) вычисление собственных векторов А: ограничимся случаем
Ао = 3
42)\
с(2)
„(2)
0.

46
Решая данную СЛАУ, получаем систему уравнений
8a?i -6х2 + 2х3 = (П
-6χι + 7х2 - 4х3 = 0 I =* χί2) = х22)/2; х22) = х32).
2χι - 4х2 + Зжз = О
Итак, любой вектор вида
/42) = 42)/2
42)
I _(2) 1 _(2)
\ х3 — х2
ж2 ' I, где х2 7^ 0 является
/0,5
собственным вектором матрицы А. Например, вектор I 1
V ι
собственный.
Далее max{|Ai| = 16, |А2| = 6, |λ3| = 3} = р(А) = 16.
Замечание. При вычислении на ЭВМ настоящего примера
необходимо оценить число обусловленности cond(A). Исходная
матрица (А) симметрична и хорошо обусловлена:
cond{A) = Amax/Amin = 18/3 = 6 < 10. ■

47
ЛЕКЦИЯ 8
Прямые методы решения СЛАУ
3.3. Решение СЛАУ (методы Крамера и Гаусса)
Традиционно СЛАУ имеет представление
ацх\ + αΐ2#2 ~\ V ο>ιηχη — £>ιν
«21^1 + α22#2 ^ 1" α2η#η = &2 .
> = Αχ = Ь,
(3.3.1)
ania?i + hn2#2 Η l· annxn = bn)
причем Ап — заполненная матрица. К прямым методам относят:
Метод Крамера. В методе при det(A) φ 0 можно записать
χι = \ΑΧι\/\Α\Λ
Xi = \A%i\l\A\,
> det(ax.) =
Xn = \AXn\/\A\,j
an «12
a<i a*2
flnl an2
ain
ь*
bn ... ar
(3.3.2)
Если det(A) φ О и 3 del^a^) ^ 0, то СЛАУ имеет единственное
решение.
Если det(A) φ 0 и для Vj Ε η : det^j) = 0, то СЛАУ не имеет
решения.
Если det(A) = 0 и для Vj € η : det^a^) = 0, то СЛАУ имеет
бесчисленное множество решений. Недостаток: вычисление det(Ax.)
требует п! операций.
ПРИМЕР 3.3.1. Решить методом Крамера СЛАУ:
2а?1 + 4ж2 + Зж3 = 4,
За?1 + а?2 — 2жз = "~2,
4а: χ + Иж2 + 7жз = 7.

48
РЕШЕНИЕ. Проводим анализ СЛАУ и затем вычисляем
неизвестные.
det(A) =
2 4
3 1
4 11
3
-2
7
= 2
1
11
-2
7
-3
4 3
11 7
+ 4
4
1
= 29^0,
= 2
I -2
II 7
-3
4 4
11 7
+ 4
4 4
1 -2
2 4 4
det(A|6) = 3 1 -2
14 11 7
= 58 ^ 0 —)· rankA = ranfcA|6 = 3 —»■ (χι, #2,#з) — независимые
переменные ->· существует единственное (!) решение СЛАУ,
1 -2
11 7
1
4 3
11 7
+
4 3
1 -2
4 4 3
det{AXl) = \-2 1 -:
I 7 11 7
= 2(7 + 22) - 3(28 - 33) + 4(-8 - 3) = 58 + 15 - 44 = 29
2 4 3
-> χι = \АХ1\/\А\ = l.det^) = '
3 -2
4 7
= 2
-2 -2
7 7
4 3
7 7
+ 4
+ 4(-8-6) = 0-21-
4 3
-2 -2
1 = -29
= 2(-14 + 14)-3(28-21)+
х2 = \АхЛ/\А\ = -1,
det(AX3)
2 4
3 1
4 11
= 2
I -2
II 7
-3
4 4
11 7
+ 4
4 4
1 -2
= 2(7 + 22) - 3(28 - 44) + 4(-8 - 4) = 58 + 48 - 48 = 58 -> х3 =
= |АЖЗ|/|А| = 2.

49
Решение хт = (χ ι = 1, χ<ι — — 1, ж3 = 2) обращает СЛАУ в
тождество. ■
Метод Гаусса (последовательного исключения
неизвестных)
Решение СЛАУ (3.3.1) методом Гаусса (последовательного
исключения неизвестных) включает прямой и обратный ходы.
Прямой ход приводит первое уравнение СЛАУ (где ац φ 0) к виду
• αΐ2 ι αΐη *ι ν ι ; ; ; τ
Χ\ λ Χ2 Η 1 Χη = > \Χ\ + Cl2%2 ~\ l· CinXn = d\
ац ац ац '
(3.1.3)
Затем из (3.1.3) выделяют х\:
Ь\ ai2 ain
a?i = Ж2 — · · · #n
ац ац ац
и исключают ж ι из всех остальных уравнений исходной СЛАУ. Из
полученного СЛАУ (Α)η_ι(:Ε2,#3> · · · ·>χη)Τ = (Ь) последовательно
берут первое уравнение СЛАУ (где α^ι Φ 0) и приводят к виду
«23, , «2п &2
Х2 + -Ζ—ΧΖ ~\ 1" -—Хп = >
«22 «22 «22
#2 + ^23^3 Η l· С2пХп — ^2,
(3.3.4)
откуда выделяют Х2\
&2 «23 «2п
#2 = Ζ 1—#3 —Жп
а22 «22 «22
и исключают ж 2 из всех остальных уравнений СЛАУ:
(i)n_i(z2,z3,...,Zn)T = (Ъ) и т. д.

50
В результате такой (прямой) процедуры исходная СЛАУ(3.3.1)
преобразуется к виду с верхней треугольной матрицей
Х\ + С12Х2 Η 1" Cin#n = d\ У
Х2-\ V С2пХп = d>2 t
%n — **n
/1 tXi2.
0 1
->ffi
0
^23-
0
U2n
1
щ
(3.3.5)
\0 0 1 /
Обратный ход включает последовательную процедуру
вычисления неизвестных a?i,a;.2>· · · >Яп СЛАУ (3.3.5) треугольной
формы, причем исключение начинают с последнего уравнения СЛАУ:
хп = ап —У %п—1 —f * · * ==Ф> #ι·
В символической форме решение методом последовательного
исключения Гаусса для прямого хода имеет вид:
Ах = LU\X — Ь ==> U\X = L~lb — прямой ход,
_1 -1 _1 _1 _1 (3.3.6)
Ux U\x = υλ L b -» χ = E/\ L b — обратный ход.
ПРИМЕР 3.3.2. Решить методом последовательного
исключения СЛАУ:
2a?i + 4x2 + Зж3 = 4,
Зж1 + Х2 - 2х3 = -2, (А)
4а?х + 11ж2 + 7ж3 = 7.
РЕШЕНИЕ. Анализ СЛАУ и det А = 29 φ 0 (пример 3.3.1).
Первое уравнение СЛАУ (А) приводим к виду (разделив на αϊ ι φ
0):
χι + 2χ2 + 1,5ж3 = 2 =φ· χι = 2 - 2ж2 - 1,5χ3· (Β)

51
С помощью выражения для х\ системы (В) исключают последнее
из второго и третьего уравнений исходной системы (А):
Г -2 = 3χι + х2 - 2х3 = 3(2 - 2х2 - 1,5х3) + х2 - 2х3 = 6 - 5х2 - 6,5х3,
\ 7 = 4xi -Ь Наг + 7хз = 4(2 - 2х2 - 1, 5х3) + Ия2 + 7х3 = 8 + Зх2 + хз-
Последняя система принимает окончательно вид
5ж2 + 6,5ж3 = 8,
{
Первое уравнение системы (С) приводим к виду:
Х2 + 1,3ж3 = 1,6 -> я2 = 1,6 - 1,3а:з, (D)
что позволяет исключить χ<ι из последнего уравнения системы
(С):
3(1,6-1,Зжз) + жз = 4,8-3,9ж3 + а;з = 4,8-2,9жз = -1,
2,9жз = 5,8-»жз = 2. (Е)
Система уравнений (В), (D), (Е) после прямого хода принимает
треугольный вид:
х\ + Ίχι + 1,5жз = 2,
*2 +1,3яз = 1,6, (F)
ж3 = 2.
Обратный ход позволяет последовательно решить систему (F),
начиная с последнего уравнения: жз = 2,Ж2 = 1,6 — 1,3жз = 1,6—
-2,6 = -1,
χι = 2 - 2х2 - 1,5ж3 = 2 - 2(-1) - 1,5(2) = 1. ■

52
Метод Гаусса (с выбором главного элемента
исключения):
Определение 3.3.1. Главным элементом А матрицы А
называется максимальный по модулю ее коэффициент (возможно не
единственный), причем существуют (возможно не единственные)
главный элемент г-й строки (обозначим А%) и (или) столбца
(обозначим Aj) в указанном смысле. Строку с главным элементом
называют главной строкой этой матрицы.
Целесообразно (доказательство оставлено читателю).
УТВЕРЖДЕНИЕ. Если VAij главные элементы %-х строк
матрицы Ап одновременно главные элементы j-x столбцов (не
обязательно \/г = j), то 3(!)РП — матрица перестановок, что
РпАп = An : VAj -» An, где An ^ Д22 · · · ^ Ann, причем An ~
матрица с диагональным упорядочением главных элементов.
Отличие настоящего метода от метода последовательного
исключения состоит в том, что на каждом этапе исключаются
неизвестные при главном элементе матрицы ее коффициентов.
По сути предварительная реализация диагонального
упорядочения исходной А позволяет при решении СЛАУ вернуться к
методу последовательного исключения.
ПРИМЕР 3.3.3. Решить методом Гаусса с выбором главного
элемента матрицы коэффициентов СЛАУ:
( 2х\ +Х2 + 4жз = 16,
< 6χι + 2х2 + Зж3 = 19,
[ χι + 5x2 + 2х$ = 17.
РЕШЕНИЕ. Определитель detA φ 0 -» rankA = rankA\b —
3 -* (#1,Ж2,#з) — независимые переменные. Реализуем
диагональное упорядочение главных элементов, где А = 6 — главный

53
элемент (Л):
(
л=[
\
2
3
1
1
2
[5J
Rh
3
2У
к /
Η (Д) =
1 V
"б]
1
2
2
Ξ
1
3
2
И
6ж1 + 2ж2 + Зж3 = 19,
χι 4- 5ж2 + 2ж3 = 17,
2χι 4- Х2 + 4жз = 16.
Исключаем неизвестное ж ι при при главном элементе:
19
19
χι = τ ~ ъХ2 - бжз = τ ~ Г2 ~ Г3
из всех уравнений СЛАУ с (Л):
19 1 1 Е 19 14 3
"б" ~ З^2 ~ 2ЖЗ + + 2жз = у + -д-ат2 + -ж3 = 17 =►
14 3 83
^у^-Н-жз^-g-,
п,19 1 1 \ 19 1
У ~ 3Х2 ~ 2ХзУ + Ж2 + Υ + З^2 + = 16 =»·
1 29
=>> зЖ2 + Зж3 = у.
После исключения система имеет вид:

54
14
Главным элементом этой матрицы является коэффициент —, по-
о
этому исключаем неизвестное χι при этом коэффициенте:
Х2 ~ 6 14 2 ЫХз ~ 28 28Хз"
Исключаем с помощью последнего уравнения хъ из первого
уравнения последней системы:
1 /83 9 \ „ 83 81 29
\\ 28 ~ 28*8) + ЗЖЗ = 84 ~ 28*8 = Τ *
81 812 - 83
=> 28ЖЗ = ~W~ * Х3 = 3·
Собирая все уравнения с главными элементами, запишем систему
с верхней треугольной матрицей:
XI +
Х2 +
3*2 +
9
28Жз =
2*3 =
_ 83
= 28'
6'
хз = 3.
Совершая обратный ход, последовательно находим неизвестные:
83 J 56 0 19 1 А 6 л
хо = 3, Жо = о— = — = 2, хл = 2 о- = — = 1.
3 ' 2 28 28 28 ' 1 6 3 2 6
Решение жт = (а^ = 1, ж2 = 2, жз = 3) обращают СЛАУ в
тождество. ■

55
3.4. Решение СЛАУ (метод предобусловливающей
матрицы)
Символическая форма Ах — Ь полезна при решении СЛАУ.
При численном решении СЛАУ с понятием предобусловливающей
матрицы связано использование такой (невырожденной и легко
обратимой) матрицы UT, которая позволяет улучшить число
обусловленности матрицы коэффициентов СЛАУ:
Ах = Ь -> К*1 Ах = К-гЬ : cond{K~гА) -> 1. (3.4.1)
Метод обращения матрицы. Символическая запись СЛАУ
Ах = b позволяет при detA ψ 0 с учетом А~гА = Ε записать
Ах = Ь -* Α~Μ* = А^Ь -+ а? = А"1*, (3.4.2)
где А = К — предобусловливающая матрица: cond(A~lA) — 1.
Решение СЛАУ методом обращения матрицы имеет развитие
с учетом LU-факторизации А"1 (2.5.3):
Ах = Ь => χ = А~гЪ = {LU)-1 Ь = U~lL-\ (3.4.3)
рациональное при больших порядках матрицы Ап:
ПРИМЕР 3.4.1. Решить СЛАУ, используя (3.4.2) и (3.4.3):
2a?i + 4ж2 + Зжз — 4,
3^1 + Х2 - 2жз = -2,
4a:i + 11ж2 + 7ж3 = 7.
РЕШЕНИЕ. Анализ СЛАУ и detA = 29 φ О (пример 3.3.1),
вычисление обратной матрицы методом присоединенной матрицы

56
(пример 2.5.1) и χ = А ХВ:
5 -11\ / 4"
χ = А~1В = ^- | -29 2 13-2
-6 -10 / V 7у
29-4- +5· (-2) + (-11) -7 = 29
(-29)-4 + 2-(-2)+ 13· 7 = -29
(29) · 4 + (-6) · (-2) + (-10) · 7 = 58,
χι = 1
Ж2 = —1
хз = 2
P. S. Решение методом обращения LU-факторизованной
матрицы (примеры 2.4.1 и 2.5.2) позволяют записать:
Lrx=| -1,5 ι ο|, ff_1 = ^
1
1,5
2,9
0 0
1 0
0,6 1
14,5
0
0
11,6
-5,8
0
-11
13
-10
что позволяет вычислить χ (3.4.3):
x = U~1Li1b =
14,5 11,6 -11
0 -5,8 13
0 0 -10
29 5 -11'
-29 2 13
29 -6 -10

57
ЛЕКЦИЯ 9
Решение СЛАУ (методы итерации)
Напомним, что традиционно итеративные методы вычисления
корня х* £ [а, Ь] непрерывной функции у = f(x) связаны с ее
представлением в форме χ = φ(χ), причем процесс итерации
xt+i =ψ{χι) (t = 0,1,2,...)
(А)
сходится с абсолютной погрешностью вычисления б к ж* — φ(χ*) =
= /(#*) = 0 при любом начальном жо> если ψ{χ) определена,
дифференцируема и
3ς<1, 4Τθ\φ'{χ)\ζ4, (Β)
причем фактически условие (В) вытекает из (А).
3.5. Решение СЛАУ методом простой итерации
Для СЛАУ {А)(х) = (Ь) (с учетом заданной абсолютной
погрешности вычисления ε), представленной в итеративной форме:
οι п о>\2 a>in
Xl = и Х2 — · · · Хп
an an ац
&2 a21 A «2n
χ2 — χι _ (J _ . . . χη
0,22 a22 a22
Q>nn
0>nl
(Ал) = -
-χι -
ι
ί О
Q21
α22
ι β*1'1
^ arm.
«η2
Or;
-ж2 -
Qi2
ац
0
Qn2
Ann
> ξ (χ) = (ьл) + Ил)И -+
Qin \
ац
Q2n
«22
0
, (Ы
ац
a>22
^ α,ηη '
(3.5.1)

58
необходимым и достаточным условием (х) —» (ж*) £ [а, Ь] является
3||Д4||г < 1j что гарантирует матрица (А) исходной СЛАУ с
диагональным преобладанием.
Определение 3.5.1. Матрицей с диагональным
преобладанием элементов по строкам называют матрицу, если
η
4Ai:\Aii = aii\>J2\aij\. (3.5.2)
г=1
Определение 3.5.2. Матрицей с диагональным
преобладанием элементов по столбцам называют матрицу, если
η
VAi : \Лц = ац\ > 22 Ιαύ'Ι· (3.5.3)
i=i
Определение 3.5.3. Матрицей с диагональным
преобладанием элементов называют матрицу, у которой УЛп преобладают по
i-й строке и г-му столбцу одновременно.
Определение 3.5.4. Пусть
Г ЗХ* €М\
\ 3XW е Rn,k = 0,1,2,... .
Тогда говорят, что последовательность {Х^} сходится к X*,
если
lim ||Χ^-Χ*||=0.
к~*оо
Замечание 3.5.1. Сходящиеся процессы итерации
обладают свойством самоисправляемости: погрешности округлений при
итерационных вычислениях не накапливаются. Поэтому
погрешность окончательного результата, удовлетворяющая заданной
абсолютной погрешности вычисления ε, — это погрешность
округления на последнем шаге итерации: (x)k+i = {Ьа) ~ {Ал){х)к·

59
Число шагов к итераций, обеспечивающих заданную точность
ε вычисления корней СЛАУ, определяется соотношением
<ε.
(3.5.4)
Из соотношения (3.5.4) следует: число шагов итерации к
уменьшается вместе с ЦД4Ц, причем для ускорения процесса итерации
естественно брать наименьшую норму.
ПРИМЕР 3.5.1.Решить методом итерации с точностью е < О,01
СЛАУ:
χι + 2х2 + 4х3 = 31,
Ъхг +Х2 + 2х3 = 29,
3xi — Х2 + ХЗ = Ю.
РЕШЕНИЕ. Проверяем определитель:
\А\
1 2 4;
5 1 2
3 -1 1
= 1(1 + 2) - 5(2 + 4) + 3(4 - 4) = 2 - 30 + 0
= 1
1 2
-1 1
-5
2 4
-1 1
+ 3
2 4
1 2
-28 φ 0,
Χι + 2x2 + 4х3 — 31"
5ж1 +Х2 + 2х3 = 29
3χι -Х2 + #з = 10 j
=» <
Г а?1 = 31 - 0 - 2х2 - 4жз
х2 = 29 - 5ж1 - 0 - 2х3
[ хз = Ю — 3χι + Х2 — 0
^АА = -1
=*·
5
.3
2 4
0 2
-1 0
Очевидно, что У||-Ал||< > 1, поэтому исходная СЛАУ не реша-
бельна методом итерации. Переформировываем систему с учетом
диагонального преобладания

60
5Jxi + X2 + 2жз = 29
Sx{ — X2 +X3 = 10
Ж1 + 2x2 +14|ж3 = 31
Ж1 = 14,5-0-0,2x2-0,4x3^
> =Ф> x2 = -10 + 3χι - 0 + X3 >
х3 = 7,75 - 0,25χι - 0,5х2
=>АЛ =
О 0,2 0,4'
-3 0 -1
0,25 0,5 О
В последней СЛАУ второе уравнение системы "плохое":
именно оно нарушает существование диагонального преобладания и
"портит" IIA^Ht СЛАУ. Тогда, объединяя
5χι + х2 + 2х3 = 29 (-2,8)'
=* χι - 7,8х2 - 0,6х3 = -31,2.
3χι - Х2 + хз = Ю (5)
Получаем СЛАУ с диагональным преобладанием:
5^ι+χ2 + 2χ3 = 29Ί
χι
7,8 Ь2-0,6x3 = -31,2
χι +2х2 +[4Ьз = 31;
а?1 = 14,5 - 0 - 0,2х2 - 0,4х3 ^
Х2 = 4+—Х!-0-—Х3
х3 = 7,75 - 0,25χι - 0,5х2 J
0 0,2 0,4
-1/7,8 0 0,6/7,8 j ,
0,25 0,5 0
> =Ф·
Ъа

61
Тогда
1 f\
||Лл||/ = тах{0,6,-^-, 0,75} = 0,75 < 1,
7,8
ν, м г 1 1112 0,6Ί η „ Λ
||Ал|Ь = тах{- + Г,- + - + - + ^} = 0,7<1.
Полученные нормы позволяет сделать вывод о возможности
решения последней СЛАУ методом итерации (подробное решение
методом итерации рассмотрено в примере 3.6.1). Ограничимся
вычислением с учетом (3.5.4) числа шагов к итераций для
достижения заданной точности е ^ 0,01:
^J5^1^5 < 0,01 => 0,75fe+1 < 0,000172 =*
1 — 0,75
, ^ In 0,000172 , -8,668 , пп лп , ял
ОТВЕТ, к « 29. ■
Замечание. В принципе возможно такое формирование
матрицы (Аа — Е), что к = 1, т. е. (А) — диагональная матрица.
Решение СЛАУ методом Зейделя
По существу решение СЛАУ методом простой итерации
происходит по схеме (x(t + 1)) = (Ъа) + {Аа)(х^)), причем Уж* £ (я(*))
вычислены на £-м шаге итерации:
η
Vxi{t) e {x{t)) : Xi(t + 1) = Y^aijXjit) + bi(t). (3.5.5)
i=i
Если по мере вычисления компонент χ £ (x(t + l)) учитывать их
в (x(t)) на t + 1-м шаге итерации, то
г η
Xi(t +1) = Y^aijXj{t +1) + ]Г ααχΛ*) + Μ*)> (3·5·6)

62
причем такой итеративный алгоритм решения СЛАУ называют
методом Зейделл.
ПРИМЕР 3.5.2. Решить методом Зейделя СЛАУ примера 3.5.1:
х\ +2ж2 +4жз = 31,
5^1 +х2 + 2хг = 29,
Зхг - Х2 + ХЗ = Ю.
РЕШЕНИЕ. Используя решение примера 3.5.1, выпишем
det(A) = —28 /Ои СЛАУ с диагональным преобладанием, уже
подготовленную для итерации
XI
J5ja?i + х2 + 2х3 = 29
7,8]ж2-0,6жз = -31,2
a?i + 2ж2 +1 4 Ьз = 31
хг = 14,5-0-0,2ж2-0,4жз'
1 л с
Жз = 7,75 - 0,25ari - 0,5ж2,
0 0,2 0,4 \ /14,5'
=» Ад = -|-1/7,8 0 0,6/7,8 , ЬА = 4
0,25 0,5 0 У \7,75
Тогда
||Лл||/ = шах{016,^|,0,75} = 0>75<1,
и, и г 1 1112 0,6, п „ ,
l|A^IUz = max{- + i,- + - + --f^} = 0,7<l,

63
на примере которой рассмотрим первые два шага решения
методом Зейделя:
t = h {
{ ая(1) = 14,5[14,5],
4 + 1 45 7
х2(1) = V-^14,5 = -^ « 5,859 [4],
7,8
7,8
45,7
х3(1) = 7,75 - 0,25 · 14,5 - 0,5 · -~ ю 1,19551 [7,75],
V 7, о
хх(2) = 14,5 - 0,2 · 5,859 ~ 0,4 · 1,19551 и 12,85 [10,6],
οι О 1 Π Λ
t = 2: < ж2(2) = —^- + 12,3714- — 1,19551 « 5,5555 [5,2628],
7,8 7,8 7,8
[ х3{2) = 7,75 - 0,25 ■ 12,3714 - 0,5 · 5,55548 » 1,75976 [2,125].
причем в квадратных скобках указаны для сравнения результаты
расчетов методом простой итерации. ■

64
ЛЕКЦИЯ 10
Решение специальных СЛАУ
Определение. Специальными назовем СЛАУ {А)тХп(х) = (6),
у которых пфт или (и) (А) — разреженная матрица.
3.6. Методы решения СЛАУ (η φ т)
Линейная алгебра рассматривает решение СЛАУ
η
Y^aijXj = bi (i = l,m/n) = (A)mxn(x) = (ft).
i=i
Решение определенной (η = га) неоднородной СЛАУ единственно
при det(.A) φ 0. Применяя к переопределенным (га > п) СЛАУ
метод Гаусса или определяя ее ранг, в принципе приходим к
совместной системе, для которой η = га' или η > га''. Случай η = га' тра-
диционен. Рассмотрим решение недоопреленной (п > га') СЛАУ.
Определение 3.6.1. Если определитель (минор),
составленный из га коэффициентов при переменных СЛАУ (п > га),
отличен от нуля и удовлетворяет теореме Кронекера—Копелли (3.1.1),
то такие переменные и их минор называют базисными
(основными), а остальные (п — га) переменные называют неосновными.
Определение 3.6.2. Базисным решением СЛАУ (п > га)
называют всякое решение X основных (базисных) переменных через
неосновные переменные, причем базисные переменные называют
допустимыми, если они удовлетворяют ограничениям на
переменные СЛАУ.
Итак, СЛАУ (А)тхп(х)т= (Ь) совместна (3 ОДР СЛАУ),
если ранги матрицы (А) и расширенной матрицы {{А\Ь) совпадают:
г{А) = г{А\Ъ), причем СЛАУ
— имеет единственное решение, если г (А) = п,
— имеет множество решений, определяемое решением базисной
подсистемы (базисного минора) СЛАУ, если г(А) < п.

65
Число способов выбора т < η переменных из η вычисляют по
формуле сочетаний С™ — п!/[га!(п — т)!].
Решение СЛАУ (η φ га) методами линейной алгебры
использует теорему Кронекера—Копелли (3.1.1), а также элементарные
преобразования (перестановка двух строк, умножение строки на
ненулевое число, прибавление к строке другой строки,
умноженной на число) и дополнительные базисные переменные, которые
приводят исходную матрицу (А) СЛАУ к ступенчатой форме:
Определение 3.6.3. Каноническая ступенчатая форма
{Uech)mxn матрицы {А)тХп определяется следующей структурой
[14]:
(i) всякая нулевая строка (состоящая из нулей) расположена
ниже всех ненулевых строк;
(ii) первая единица, стоящая на г-м месте всякой j-й
ненулевой строки при условии, что остальные элементы г-го столбца —
нули, называют ведущей;
(iii) все ведущие единицы образуют ступенчатую
конфигурацию слева направо (по структуре верхней квазидиагональной
матрицы), например, ведущая единица г-й строки всегда расположена
левее ведущей единицы г + 1-й строки. Например, матрица (В)
представлена в канонической ступенчатой форме:
(*) =
/О \Т\ 0 -1 0 3\
I О О |Т] 2 0 3 1
0 0 О О \Т\ 3 I
\о о о о о о/
= (Uech)4x6,
где ведущие единицы обведены рамкой.
Замечание 3.6.1. Очевидно, что (Uech)mxn — единственна.
Для рациональных вычислений достаточно пользоваться
ступенчатой формой (Aech)mxn матрицы (-A)mxn, определение которой
исключает (ш)-й признак определения 3.6.3, поэтому (Aech)mxn
не единственна.

66
Замечание 3.6.2. Подматрицей из ведущих единиц (Uech)mxn
является (E):det(E) = 1, что обосновывает существование
базисных переменных СЛАУ (А)тоХп(а;) = (Ь) (определение 3.6.1). В
принципе допустимо использование —1 в качестве ведущих
элементов (Uech)mxn·
ПРИМЕР 3.6.1. Найти общее решение СЛАУ:
( —1х\ - 6x2 + Зхз + 4ж4 = 5,
4a?i + 12ж2 + 2хз — 8ж4 = 3,
4ж1 + 12х2 + 2жз + 5^4 = 2.
РЕШЕНИЕ. Здесь п = 4>т = З^С$ = 4-3- 2/(3!)
Последовательно формируем миноры второго порядка:
= 4.
-2 -6
4 12
1) (жьх2) ··
2) (x2,x3):det(B)
= 0,
-6 3
12 2
-48 φ 0.
Итак, гапк(В) = 2, причем для расширенной матрицы
-6 3 15
(В\Ь)
12 2 |3
->· гапк(В\Ь) = 2:
например
3 5
2 3
= -1^0.
Выбираем матрицы, окаймляющие (В) (метод окаймления):
-6 3
12 2
(С) =
2 -6
12
12
-6 3
12 2
12 2
3
2
2
4
8
5
= 0 -> Ж1,Ж2,Яз,
φ 0 -> Ж2,#3,#4,

67
причем гапк(С) = rank(C\b) = 3, т. к. для расширенной
матрицы
-6 3 4 |5
(С\Ь) = J 12 2 8 |3 | -»■
12 2 5 12
3 4 5
2 8 3
2 5 2
= -19?έ0.
Записываем СЛАУ с базисными переменными Х2>жз>ж4 в форме
-бжг + Зжз + 4х4 = 5 + 2xi(i),
12i2 + 2хз — 8x4 = 3 — 4xi(w),
12^2 + 2жз + 5x4 = 2 — 4жх(ш).
Последовательно решая базисную систему:
(И) - (г): -13x4 = 1 -)·
2(г) + (ш): 8ж3 = 13
Х4 =
1
13
хз
13
8
13 —1 402
(it): 12x2+ 2 j-8— = 12x2 + — =3-4х! -+
Х2 =
15
"208
χι
3
найдем общее решение: X — (жι, —щ — ^, ^, — ^) .
Аналогичное решение СЛАУ получено, используя элементарные
преобразования:
-6 3 4 |5'
12 2 8 |3
12 2 5 12
'0 0 8 13 112
= | 0 0 0
к4 12 2
0 0 8 0
0 0 0 -13
16 48 0 20
-2-6 3 4 |5'
0 0 0 -13 |1
4 12 2 5 |2
/0080 |13'
= 0 0 0-13 |1
\4 12 2 5 |2
0 0 8 0 113
0 0 -13 |1
208 624 0 0 1-45

68
дополнительные базисные переменные х$ = 0,^6 = Ο,χγ = О:
-»
-2х\ - 6^2 + Зхз + 4^4 + \хъ I + 0^6 + 0^7 = 51
4χι + 12^2 4- 2#з - 8^4 + Ожб + ГЁб| + Οχγ = 3
4a?i + 12^2 + 2жз + 5^4 + 0a?5 + Ожб + Г#7| = 2 J
> -»
1
0
0
0 0
1 0
0 1
^о->
4 1 0 0|5
-8 0 1 0|3 | -+
5 О О 112
хь — 5 + 2a;i + 6ж2 — Зхз — 4^4,
аг6 = 3 - 4а?1 - 12ж2 - 2хз + 8ж4, ->· —
Ж7 = 2 — 4х\ — 12x2 — 2жз ~ 5а?4·
7 2
^2
3?4
жз
15
" 208
1-
_ 13
~ 8
1 1
~ Г1 + 24Ж5
- Xq+X7
13
1 1
4 8
26-24
#6
39
ж7)
< Х4 —
3.7. Решение СЛАУ методом прогонки
В связи с интенсивным развитием вычислительной техники
возрос интерес к разреженным матрицам, т. е. матрицам,
большинство элементов которых — нули.
Разреженные матрицы порождают проблемы
вычислительного характера, а также проблемы по их размещению и хранению
в памяти ЭВМ. Примером разреженных матриц служат трех-
диагональные матрицы (А)п уже при η > 6. Трехдиагональные
матрицы возникают, например, при численных решениях
дифференциальных уравнений (см. Численное решение диффереци-
альных уравнений). Метод решения СЛАУ с трехдиагональной

69
матрицей А аналогичен методу Гаусса (с последовательным
исключением), однако существует особенность: каждое уравнение с
исключаемым неизвестным подставляется не в СЛАУ, а только в
одно последующее уравнение (прогонка), что связано с
отсутствием исключаемого неизвестного во всех других уравнениях СЛАУ,
кроме следующего.
Метод прогонки (аналогично методу последовательного
исключения Гаусса) включает прямой и обратный ходы, причем в
символической форме прямой ход (собственно прогонка) имеет вид:
Ах = Ь ^ LUx = Ь^ L^LUx = L^b =>Ux = 1Г\ (3.7.1)
причем матрица [/-верхняя двухдиагональная.
Решение СЛАУ с трехдиагональной матрицей (А) методом
последовательного исключения называют прогонкой.
ПРИМЕР 3.7.1. Решить СЛАУ рационально:
[1)х\ - 2х2 = 5,
{П)х1+х2-3х3 = 2,
(III) - Зх2 + 2х3 - хА = 1,
(17//) - Ъхъ + χα = 8.
РЕШЕНИЕ. Матрица коэффициентов исходной СЛАУ трех-
диагональная
/1-2 О О
[А) ~ I 0 -3 2 -1
\0 0 -5 1
что позволяет решение методом прогонки. В соответствии с
методом Гаусса последовательного исключения выражаем χ ι из (I)
уравнения СЛАУ:
(Ι):\χι - 2х2 = 5 -> χι = 5 + 2х2 -> (II);

70
(И): 2 = 5 + 2х2 + Х2 ~ Зж3 = 5 + Зх2 - Зж3
Х2 - #з = -1 -> #2 = -1 + хз -* (HI);
(III): 1 = -3(-1 + х3) + 2хг - х4 = 3 - х3 - ха ~+
ж3 + а?4 = 2 l·-^ а?з = 2 - ж4 -> (IV);
(IV): 8 = -5(2 - хА) + х4 = -10 + 6ж4 -> 6ж4 = 18 ->
Х4
В методе прогонки прямой ход (собственно прогонка)
завершает реконструкцию исходной СЛАУ в форму с верхней двухдиаго-
нальной матрицей (U):
(U)(x) = (L)-l(b)
Χχ — 2X2 = 5,^
{П)х2 -хз = -1,
{ΠΙ)χ3 + χ4 = 2,
(ΙΙΙΙ)Ζ4 = $)
(U)
При обратном ходе решаем полученную СЛАУ (U)(х) = (L) 1(Ь),
начиная с последнего уравнения:
(IV):x4 = 3^
(III): 2 = ж3 + #4 = #з + 3 -> жз == -1
(II): -1 = х2 - хз = Х2 + 1 -> #2 == -2
(1):5 = Ж1 -2ж2 = a?i - 2(-2) = я-1 + 4-*Я1 = 1,
(*) =

71
что в символической форме эквивалентно последовательному
представлению
(с/)(х) = щ-\ь) -► (uy'iu^x) = (игНЦ-Чъ) -> (χ) =
= (U)~\L)-\b).
Проверка подтвердила справедливость полученного решения:
(I):si-2a2 = l-2(-2) = 5,
(И): χι + х2 - Зж3 = 1 + (-2) - 3(-1) = 2,
(III): -Зх2 + 2хг -х4 = -3(-2) + 2(-1) - 3 = 1,
(IV): -5аг3 + χα = -5(-1) = 3 = 8. ■
Возможно исходная СЛАУ (А)(Ь) = (х) записана в форме, к
которой метод прогонки непосредственно не применим (например, (А)
не является треугольной матрицей, хотя в принципе допускает
такое представление). Поэтому над (А) необходимы
преобразования (например, переход к новым переменным), которые позволят
переписать СЛАУ в форме пригодной для решения методом
прогонки.
ПРИМЕР 3.7.2. Решить рационально СЛАУ:
Г-32/1+г/2+ 2/4 = 2^
2/2 - 2у4 = 5
-5yi + 2/з = 8
22/1 -2/3-32/4 = 1,
(А)(у) = (Ь).
РЕШЕНИЕ. Матрица (А) исходной СЛАУ
(А) =

72
формально записана не в трехдиагональной форме.
Переименование переменных с учетом структуры (A): y*i ~> #χ, у 4 —> #2,
У\ —* #з> Уз "~* ж4 позволяет в новых координатах переписать
СЛАУ в форме:
( (I):a?i-2ж2-5,
I (И):жх +х2 - Зжз = 2,
| (Ш):-Зя?2 + 2я:з-я?4 = 1,
[ {W):-5x3 + x4 = 8.
для которой матрица коэффициентов приведена к
трехдиагональной форме (пример 3.7.1).
Для матрицы (А) существует другая схема переименования
переменных: у2 -* Z4, У4 -* *3, У1 -* Z2, УЗ -> *1-
1(1): 21-5*2 = 8,
(II):-*i+2*2-3*3 = l,
(П):-3*2 + *з + *4 = 2,
(П):-2*з + *4 = б.
которая читает СЛАУ примера 3.7.1 с конца. ■

73
ЛЕКЦИЯ 11
Решение нелинейных уравнений (метод секущих)
ГЛАВА 4. Численное решение нелинейных уравнений
Точное вычисление корней нелинейного уравнения достаточно
сложной структуры не всегда возможно, а часто не нужно,
например, если параметры уравнения заданы приближенно.
Задача приближенного (численного) вычисления корней нелинейного
уравнения f(x* £ Щ методом итерации с точностью ε включает:
1) (аналитическую или графическую) локализацию (отделение)
корней: выбор множества [а^, Ь{ : г = 1,2,...], каждому из которых
принадлежит лишь один корень, х* £ [а*, 6*], причем для удобства
изложения существа метода /(ж* £ [α»,6ί]) особенностей не имеет
(нет точек разрыва функции, f'(x) и f"{x) не меняют знак и др.);
2) вычисление корней х\ £ [a^fy] с заданной точностью ε.
4.1. Метод секущих (хорд)
Итеративным методом численного решения нелинейных
уравнений и их систем служит метод хорд, интуитивная идея
построения изначально эмпирической модели которого ясна из рис. 4.1.1
Рис. 4.1.1. Метод секущих (хорд)

74
Процедура построения хорд АВ, СВ, Ό В и других для
вычисления f(x* Ε [α, b]) = 0 формирует сходящийся процесс итерации
метода хорд: условие монотонности f(x € [α, b]) необходимо и
достаточно, чтобы х(0) = а, и точки ж(1),ж(2),ж(3),...
пересечения хорд с осью Ох образовали сходящуюся монотонную
последовательность х(0) = а < х(1) < х(2) < · · · < x(fe+1) < · · · ^ х* < Ъ
доказывает
УПРАЖНЕНИЕ 4.1.1. Записать рекурсивно алгоритм метода
хорд вычисления корня функции f(x* € [α, Ь]) (рис 4.1.1).
РЕШЕНИЕ. Дополним рис. 4.1.1 построениями, что позволяет
записать
1 шаг: из подобия треугольников АРЕ и ΑΒΑ' следует
АЕ^_АА^_ х{1) - ж(0) _ Ь - х(0)
^Р~А7В~* -/(*(0)) " f(b) - /(α) ~*
->«(')-*№-дДод/Мо».
2 шаг: из подобия треугольников CKF и СВС следует
С^_СС_ а?(2)-а?(1)_ ft-s(l)
FK'C'B^ -/(x(l)) "/(ft)-/(I)-»1
3 шаг: из подобия треугольников DNM и DBD' следует
DA£_DD[ ж(3) - а?(2) _ 6 - а?(2)
Ш~В~&^ -/(х(2)) ~ /(b) - /(2)
^Ж(з)^(2)-/(ь;:^(2))/и2)),...
-»

75
что по индукции допускает обобщение
/(а) < 0,/(b) > 0 : f(b)f(a) < 0) φ + 1) = х(к) - ^"^УУчУ =
/(«.)/"(*)< О-►*(<)) = Л , /(*)-/(«(*))
» I f{b)x(k)~f(x(k))b
fW»(X) > О - х(0) = aj = /(»)-/(,(*)) '
(4.1.1)
причем
lim ж(А; 4-1) — lim a;(fc) = 0 =
fc—>оо fc—юо
6 — х(к) ,. ,. ....
= - lim —7- ,. .... hm f[x(k)) ->■
fe-»oo f(b) - f{x{k)) fc->oo' v v ;;
-> 0 = lim /(x(fc)) = /(**),
I™ w,4 ~ У/, чч ^ 0, т. к. b > x*. ■ (4.1.2)
fc-+oo f(b) - f(x(k)) Г У '
Замечание 4.1.1. Возможен иной вариант итерации (вывод
аналогичен (4.1.1))
/(a) >0,/(Ь)<0-*/(*)/(<>)< 01
/(a)/"(х) > 0 -+ *(0) = Ь \х(к + 1) = а- Х[> * /(a), (4.1.3)
/(b)/»(b)<o^(o) = J /(»(*))-/(«)
ПРИМЕР 4.1.1. Вычислить методом хорд с точностью ε ^
^ 0,001 корни функции /(х) - χ3 - 27х2 + 180х ~ 324 (пример
3.2.2).
РЕШЕНИЕ. Вычисляем производные функции: /'(х) = Зх2 —
54х + 180, f"(x) = 6х — 54. Предварительно аналитически
оцениваем количество и границы корней /(ж) = х3 — 27х2 + 180х — 324
по изменению знака f(x):
Зх2 - 54х + 180 = 0 -» х^>2 = 9 ± 4,58 ->
xj 6 [-оо; 4,42],
xi « 13,58,
х2 ~ 4,42,
«5 €[4,42; 13,58],
хз€ [13,58;оо]

76
и результаты заносим в таблицу:
X
№
/'(*)
f"(x)
0
-324
180
-54
2
-64
84
-42
4
28
12
-30
5
26
-15
-24
10
-224
-60
6
14
12
Используя f(a)f(b) < 0, приведем уточненный расчет отрезка
существования корня х\ £ [—оо;4,42]:
f(x = 4) = 43 - 27 · 42 + 180 · 4 - 324 = 28 > (Г
f(x := 2) = 23 - 27 · 22 + 180 · 2 - 324 = -64 < 0
-»
-»ж* е [а = 2;Ь = 4],.
причем /"(2) = 6 · 2 - 54 = -42 и /"(4) = 6 · 4 - 54 = -30,
следовательно
f(a)f"(x) > 0 -)· х(0) = 6 : /(2)/"(2) = (-64)(-46) > θ)
f(b)f"(x) < 0 -* х(0) = Ь : /(4)/"(4) = (28)(-38) < 0J ~*
-» я?(0) = Ь = 4 -*· /'(ж = 4) = 4(3 · 4 - 54) + 180 = 12,
что позволяет вычислить итеративно х* при х(0) — Ь = 4 и е ^
< 0,001:
η хШ - α _ *(°) ~ α f (a) - 2 _ (4-2)(~64) „ 3 3913.
1,ΐ(Ι,"β /(χ(0))-/(α)/(α)~2 28-(-64) ~3'3913'
2) /(3,3913) = 14,9125,
3) /(3,12838) = 5,48254;

77
x(2) — a
Ж(3) = а-/(ж(2))-/(а)/(а^3'°39345;
4) /(3,039345) «1,74272;
5) /(3,011794) = 0,52823;
Ж(5)-а-/(х(4))-/(а)/(а)"3'0°351;
6) /(3,00351)^0,157723;
s(6) = а - *^j ~ ° /(а) * 3,0010 < β = 0,001;
Да?(5)) - /(а)
что удовлетворяет ε ^ 0,001. Полезен вычислительный
эксперимент с расчетом х* € [2,4] по (4.1.1) при х(0) = α = 2:
2) /(3,3913) = 14,9125,
ж(2) - х(1) - /(»W)(b-Hl)) и 2 6977.
Ж(2)~Ж(1) /(b) -/(х(1)) ~2'Ь977'
3) /(2,6977) = -15,2759;
χί3) - Х(2) - f{x{2)){b ~ {Х{2)) » 3 15739·
Х(3)-Ж(2) f(b)~f(x(2)) ~3>15739>
4) /(3,1579) = 6,6406;
5) /(2,8954) = -4,90486;

78
{)~ () /(6)-/(x(4)) ~3'0601'
6)/(3,0601) = 2,6397;
x(6) = *ί5ϊ - /(*(5))(*-(*(5)) „ %22
*W X{5} f(b)-f(x(b)) ~2>96227'
7)/(2,96227) = -1,72355;
8) /(3,0224) = 0,99898;
9) /(2,9862) = -0,6244;
10)/(3,00831) = 0,3727;
И)/(2,99493) = -0,2286;
.(η)^(,ο,-^^-Η;;)), .cos,,
12) /(3,00307) = 0,13798;

79
13) /(2,99813) = -0,08421;
χίш - *fш - /И12Жь" И12)) „ з nous-
x(13)-z(12) /(6) _/(χ(12)) -3,00113,
14) /(3,00113) = 0,050827;
15)/(2,99931) = -0,031059;
Wi«rt r(u\ /Н14))(^-И14)).,?ппп.0
«(15) = *(14) - /(ь) _ /(ж(14)) » 3,00042;
что с заданной погрешностью соответствует точному значению
х* — 3 и интуитивно обосновывает использование итеративного
алгоритма метода хорд (4.1.1) как для расчетов при х(0) = а, так
и при х(0) = Ь. ■
Уже анализ полученного численного решения показывает, что
при указанной модификации метода хорд происходит:
— последовательное сужение отрезка корня как на левой, так и
на правой границах уточненного отрезка корня χ* [α, Ь], что
естественно увеличивает число шагов итерации для получения значения
корня с заданной ε;
— левая и правая последовательности сужения отрезка корня
χ* Ε [α, b] монотонны:
а = ж(0) < αϊ < α2 < α3 < · · · ^ ж* ^< · · · < Ь3 < Ъ2 < h < Ъ.
(4.1.4)
Для целей дальнейшего анализа воспользумся графической
реализацией приведенного численного решения (рис. 4.1.2):

80
У=Ях)
РИС. 4.1.2. Модифицированный метод секущих (хорд)
Геометрический анализ показывает, что происходит
построение хорд по несколько иной схеме, что не нарушает итеративный
процесс монотонного сужения отрезка корня х* (с обоих границ)
вплоть до вычисления х* с заданной точностью ε.

81
ЛЕКЦИЯ 12
Решение нелинейных уравнений (метод касательных)
4.2. Метод касательных (Ньютона)
Идея метода Ньютона следует из определения 3.5.1 и рис. 4.2.1,
на котором изображен фрагмент функции f(x) и итеративная
процедура поиска корня f(x) на [а, 6], используя касательную /'(ж):
/{х*) = ф{х*)-х* = 0->
х* е [а,Ь],
f(a) < 0, f(b) > 0.
Рис. 4.2.1. Метод касательных
ТЕОРЕМА 4.2.1. Если f{a)f{b) < 0, причем f'(x) φ 0, /"(ж) φ Ο
и сохряняют знак при \/х 6 [а, Ь], то с учетом нулевого прибли-
жения ж(°) Ε [а, Ь], удовлетворяющего условию ff(x^)f^(x^ >
> 0, можно вычислить методом Ньютона
x(k + l) = x(k)
/(*(*))
Ж*))
корень х* функции f(x) с любой степенью ε.
(4.2.1)

82
Доказательство. В соответствии с рис 4.2.1 можно записать:
(4.2.2)
(О /(а) < 0, f(b) > 0;
(it) Vz € [а,Ь] : f'{x) > 0,/"(х) > 0;
(m) /(ж(0)) > 0, например ж(0) = 6.
Тогда для х* < Ck < х^к' по формуле Тейлора следует выражение
о = /м=/(х«) + *Цг^л**')+^р^гы +...,
которое с учетом f"(ck) > 0 принимает вид утверждения теоремы
ж < а;
(*)
/(*<*))
= χ
(fc+l)
причем всегда ж^4"1) > ж*, т. е. процесс итерации сходится. ■
Следствие 4.2.1. Замечание 3.3.1 позволяет рассматривать
f'(x(k)) в Качестве абсолютной погрешности (вычисления)
корня х* G [а, Ь] исследуемой функции f(x), причем если х(к +1) —
точное значение корня, то
/(*(*))
/'(*(*))
= |Аа?*|<е.
(4.2.3)
СЛЕДСТВИЕ 4.2.2. Существует упрощенный метод
Ньютона, связанный с заменой f'(x^) на f(x^) в алгоритме (4-1.1):
x(k + l) = x(k)
/(*(*))
ГШУ
(4.2.4)

83
Рис. 4.2.2. Упрощенный метод касательных
Идея упрощения (4.2.4) очевидна (рис. 4.2.2) и полезна (если
f(x), /'(#) достаточно сложны для вычислений), хотя и
увеличивает число итераций сравнительно с (4.2.1).
ПРИМЕР 4.2.1. Вычислить корни функции f(x) = хг - Tlx2 +
180ж — 324 примера 4.1.1 методом касательных с точностью ε ^
< 0,001.
РЕШЕНИЕ. Начальная процедура соответствует примеру 4.1.1
(приведена для удобства):
— вычисляем f'{x) = Зх2 - 54ж + 180, f'(x) = 4ж - 54;
— оцениваем границы корней f(x) = χ3 — 2Ίχ2 + 180ж — 324,
используя f'(x):
Зх2 - 54ж + 180 = 0 -* x*h2 = 9 ± 4,58 -*
χι « 13,6,
Х2 « 4,4
х\ е [-оо; 4,4],
*$ €[4,4; 13,6],
а?з £ [13,6;оо],
и заносим в таблицу:

84
χ
! /(*)
/'(*)
I /"(*)
0
-324
180
-54
2
-64
84
-46
4
28
12
-38
5
26
t · ·
10
Уточненный расчет отрезка существования корня х\ € [4,4;
13,6] с учетом f(a)f{b) < 0:
f(x = 5) = 53 - 27 · 52 + 180 · 5 - 324 =)
= (5(5(5 - 27) 4-180) - 324 = 26 > 0
f{x = 4) = 43 - 27 · 42 + 180 · 4 - 324 = 28 > 0
f(x = 2) = 23 - 27 · 22 + 180 · 2 - 324 *= -64 < 0
-*ж* € [а = 2;Ь = 4],
причем /"(2) = 4 · 2 - 54 = -46 и /"(4) = 4 · 4 - 54 = -38,
следовательно
/(α)/"(α) > 0 -)· х(0) = α : /(2)/"(2) = (-64)(-46) > 0Ϊ
W(fr) > 0 -> *(0) = Ь : /(4)/"(4) = (28)(-38) < 0J ~*
-> а?(0) = β = 2 -»· /'(ж = 2) = 2(3 · 2 - 54) + 180 = 84
позволяет вычислить итеративно значение корня х(0) = а = 2 с
точностью ε ^ 0,001:
2) /(2; 7619) = -11,7484, /'(2,7619) = 53,7417,
ж(2)=х(1)~/Ш"2'9805;

85
3) /(2,9805) = -0,88435, /'(2,9805) = 45,70314,
*(3) = x(2)~jr^ = 2,99985;
4) /(2,99985) = -0,0067504, /'(2,99985) = 45,00540,
,(4) = ,(3)-М«2,99999;
5) /(2,99999) = -0,00045, /'(2,99999) = 45,00036,
я?(5) w 2,99999 + 0,00001 = 3.
Значение абсолютной погрешности 0,0001 ^ ε = 0,0001 позволяет
прекратить процесс итерации на пятом шаге итерации.
Полученное значение соответствует значению х* = 3, что подтверждает
проверка.
Сравним процессы итерации метода Ньютона и его
модификации:
1)ж(1)=ж(0)-^ш=2-7619;
2) /(2,7619) = -И, 74845; х{2) = х{\) - ^^ « 2,90176;
3) /(2,90176) = -4,59545; х{Ъ) = х{2) - ^(0)) = 2'95647;
4) /(2,95647) = -1,99306; х{4) = ж(3) - $Ш « 2,98020;
/'(ж(0))
5) /(2,98020) = -0,89809; *(5) = а(4) - f„ ,ff4 « 2,99089;
6) /(2,99089) = -0,41143; ж(6) = я?(б) - ^Ш « 2,99579;
/ \ж\"))
7) /(2,99579) = -0,189787; х(7) = я(6) - {, ,?) « 2,99805;
/'(ж(0))

86
8) /(2,99805) = -0,087824; х{8) = х{7) - &^г ~ 2,9991.
Значение абсолютной погрешности ε = 0,0001 позволяет
прекратить процесс итерации в модифицированном методе Ньютона на
восьмом шаге итерации сравнительно с пятым шагом в методе
Ньютона. ■
Замечание 4.2.1. Предварительным сужением х* € [а, Ь]
можно уменьшить число шагов итерации в методе для достижения
заданной точности ε вычисления корня. Однако предварительное
сужение отрезка существования корня часто оказывается более
рутинным процессом.
Замечание 4.2.2. По существу численные методы
(касательных и секущих) решения нелинейных уравнений связаны с
выбором моделей (рис. 4.1.1, 4.2.1 и 4.2.2), которые достоверно и
наглядно обеспечивают сходимость процесса итерации.
Замечание 4.2.3. Итеративный алгоритм метода хорд
вычисления корня f(x) на [а, Ь] является частным случаем алгоритма
Ньютона, если воспользоваться определением производной
(правой):
ГШк)) ~ /И*0) - /И* -1))
mfc))~ x(k)-f(x(k-l)) '
полагая всегда х(к — 1) = хс(0), например,
<к+1) - *> - Шт"х(к) - /(*wWw)-/(°(o))·
(4.2.5)
Но процесс итерации метода хорд можно улучшить (пример 4.1.1),
если строго использовать определение производной
х(к + 1) = х(к) - , " « х(к) - }(х(к))·
f'(x(k)) v y ^ v "f(x(k))-f(x(k-l))'
(4.2.6)
Выражения (4.2.5), (4.2.6) связаны с разностным представлением
производной, к изучению которых переходим.

87
ЛЕКЦИЯ 13
Разности. Разностные уравнения
ГЛАВА 5. Разностные уравнения
5.1. Сеточные функции
В вычислительной математике возникает задача
дискретизации, связанная с переходом от функции непрерывного
аргумента f(x) к функции дискретного аргумента — сеточной функции
y(xi £ X : г = 0,1,2,..., |Х|) — функции, определенной на
множестве узлов сетки X = {хг,Х2>..., ж»,...., #|Х|}5 которое
называют также сеткой, например
х%
Уг
XI
г/1
Х2
У2
%г
Уз
Хп
Уп
где Xi+i — Xi = hi = Axi — шаг сетки. Если hi = h = const, то
сетку называют равномерной, иначе неравномерной.
Замечание 5.1.1. В задаче дискретизации непрерывной
функции у = f(x) при выборе шага hi = Axi = Χ{+ι — Χχ необходимо
lim y(xi : г = 0,η) = f(x). (5.1.1)
η-»οο
Замечание 5.1.2. Возможен переход от Xi —» г = 0,1,...,
причем у(#г) = 2/(0 Ξ 2/г — сеточная функция целочисленного
аргумента.
5.2. Разности
Полезно (аналогия с производной)
Определение (разностей) 5.2.1

88
Δϊ/ί = yi+i — yi — правая разность,
^Уг = Уг — Уг-ι ~ левая разность, (5.2.1)
8yi = (Ауг+ι + Vyi)/2 — центральная разность,
где Δ, V, δ — линейные операторы.
Из определения (5.2.1) сразу следует Syi = ^O/t+i-J/t-i), &Уг =
= Vi/t+i. Приведем разности
А2у{ = А(Ауг) = Д(и+1 - й) = Δκ+ι - ΔΚ =
= 2/i+2 - 2/i+l - Ι/i+l + Уг = Уг+2 ~ 2j/i+i + Jfi = V2yi+2,
= 2/г - 2/г-1 ~ Уг-1 + Уг-2 = Уг ~ tyi-Ι + Уг-2 = &?Уг-2,
ДУЙ = Δ(!/< - 2/t-l) = &Уг ~ &Уг-1 = г/г+1 ~Уг~Уг + Уг-1 =
= й+1 - 2уг· + й_1 = ЧАу{ = Δ22/ί_ι = V22/i+i
(5.2.2)
и др. Весьма полезна следующая разность
Δ(ΐ/ί*ί)(= Δ(^)) = J/i+i*i+i - 2/i^i = (j/t + Ayi)(Zi + Azi) =
= 2/i^i + 2/tAzi + ZiAyi + Ay{Azi - y{z{ = (^ + Ay{)Azi + ZiAy{.
(5.2.3)
Замечание 5.2.1. Последнюю разность целесообразно
сравнить с дифференциалом d(yz) = з/dz + zdy во избежания ошибок
при переходе от дифференциальных уравнений к разностным [16].
ПРИМЕР 5.2.1. Построить таблицу разностей сеточной
функции y{xi):

89
Xk
Ук
0
1
1
4
2
15
3
40
4
85
РЕШЕНИЕ. Строим таблицу разностей (правых) для заданной
У (я?*):
Хк
0
1
2
3
4
Ук
1
4
15
40
85
Ау
3
11
25
45
Δ2?/
8
14
20
А*у
6
6
Δ4ί/
0
5.3. Разностные уравнения
Линейным разностным уравнением га-го порядка называют
Q>o{xi)y{xi) + a>i{xi)y{xi+i) + ··· + am(xi)y(xi+m) = φ(χ{), (5.3.1)
где ao(xi) φ 0,... ,am(xi) φ 0,ф(хг) — заданные сеточные
функции. При ф(х{) — 0 разностное уравнение называют однородным,
иначе неоднородным. В частных случаях имеем разностное
уравнение га-го порядка:
га = 1: a0{xi)y{xi) + ai(xi)y(xi+i) = ф{хг)>
га = 2: a0(xi)y(xi) + ai(xi)y(xi+i) + a2{xi)y{xi+2) = Φ{Χί)·
Перепишем линейное разностное уравнение (5.3.1), например,
неоднородное разностное уравнение 1-го порядка в разностях:
a>o(xi)y(xi) + a.i{xi)y{xi+i) =
= a0(xi)y{xi) + ai{xi)y{xi) + аг(хг)Ау(хг) =
= [a0{xi) + ai(xi)]y(xi) + ai{xi)Ay(xi) = φ{χί). (5.3.1)

90
5.4. Решение разностных уравнений
Решение разностных уравнений рассмотрим на примере
однородного разностного уравнения 1-го порядка с целочисленным
аргументом:
а0{г)у{г) + ai(i)y(i + 1) = 0, (5.4.1)
где αο(ί),αι(ΐ) — известные сеточные функции, г = 0,1.
Перепишем (5.4.1)
у(г + 1) = !_— = ~гу(г).
αϊ (г) αϊ (г)
Тогда последовательно можно записать
i _ i. v(i , χ) _ _M0 fi) = _ао(г>о(г-1)...ао(1)ао(0) ( ,
г-г. гДг + IJ- ^Ш а1фв1(* _ ι).. .αι(1)αι(0)У[ h
(5.4.2)
Рассмотрим решение неоднородного разностного уравнения 1-го
порядка с целочисленным аргументом:
а0{г)у{г) + аг{г)у(г + 1) = <p(i), (5.4.3)
где ao(i),ai(i),(/?(i) — заданные сеточные функции Vi.
Перепишем (5.4.3) в виде
(. + 1} = ^)-ao(W) = Щ _ «о£)
аг(г) аЦг) аф)

91
Тогда последовательно можно записать при г — О,1,2,...
м=т.^т
У(2)
у{г + 1) =
αι(0) αι(Ο)1
φ(1) - ао(1)у(1) φ(1) α0(1)
αι(1)
αι(1) αι(1)
>(0) οο(0)
βι(0) αι(0)
2/(0)
у>(») _ <*ο(*)
αϊ (г) αϊ (г)
ар(г - 1) Γ
αι(ϊ'-Ι) L
ν(< -1)
αι(ί-Ι)
>(1) οο(1)
L«i(l) αι(1)
¥>(0) οο(0)
αι(0) ' αχ(0)
2/(0)
(5.4.4)
ПРИМЕР 5.4.1. Решить разностное однородное уравнение 1-го
порядка
t/f+i - 3t/j = 0 ->· α0 = Ι,αι = -3.
Чему равно j/4> если уо = Ю.
РЕШЕНИЕ. Решение имеет вид t/j+i = 3j/j:
г = 0: j/i = Зу0,
г = 1: у2 = Зг/ι = 3 · 3j/0,
г = 2: уз = 3t/i = 3 · 3 · 3j/0,
j/i+i = 3j/» = 3t+1yo -> J/4 = 34j/o
81 · 10 = 810.
Приведенные решения разностных уравнений полезны, если г не
велико. Существует решение, связанное с построением, используя
подстановку у(к) = \к, характеристического уравнения
разностного уравнения m-го порядка. Так разностное однородное
уравнение 2-го порядка аоу(к) + а\у(к +1) + а,2у(к + 2) = 0 с помощью
подстановки у(к) = \к ψ 0 запишем в характеристической форме
a0Xk+aiXk+1+a2Xk+2 = 0 -> α0+αιΑ+α2λ2 = О, решение которого

92
λι,2 = (—αϊ ± ι/αι — 4αοα2)/(2α2) зависит от знака D = aj — 4αοα2:
1) D > 0 : yi(i) = Xly2(k) -> »(fc) = ciAf + c2A^,
2) £> = 0 -> λι = λ2 -> y(fc) ■= ciAfe + c2 · fcA*-1, (5.4.5)
3) D < 0 -> λι, λ2 — комплексные корни.
ПРИМЕР 5.4.2. Решить разностные однородные уравнение 2-го
порядка
y(fc + 2) - Ъу{к + 1) + 4y(fc) = 0, y(fc + 2) - Ау{к + 1) + 4y(fc) = О,
используя характеристический многочлен. Вычислить у(6), если
у(0) = 3,у(1) = 6.
РЕШЕНИЕ. Решаем характеристическое уравнение λ2 — 5λ+
+4 = О -¥ λι = 4, λ2 = 1:
y(fc) = cx · 4fc + c2 ->· y(0) = ci + C2 = 3,t/(l) = ci · 4 + c2 = 6 ->·
->-ci = l,c2 = 2-H
y(fc) = 4* + 2 , y(6) = 4b + 2 = 4098.
Решаем характеристическое уравнение λ2 — 4λ + 4 = Ο —»· λχ =
= λ2 = 2:
у(к) = d · 2* + c2k2k~l -> y(0) = С! = 3, y(l) = a ■ 2 + c2 · 1 · 2° =
= 6 ->· ci = 3, c2 = О -Н
y(fc) = 3 · 2fc Ly(6) = 3 · 26 = 192.1

93
ЛЕКЦИЯ 14
Приближение функций (полином Ньютона)
ГЛАВА 6. Приближение и дифференцирование функций
6.0. Введение. Понятия
Одной из основных задач численного анализа является
задача приближения (аппроксимации) функцией φ(χ) непрерывного
аргумента некоторого множества пар точек, заданных сеточной
функцией y(xk : к = О,1,2,..., η}, χ £ [α, Ь], полученной
— дискретизацией известной непрерывной функции у = /(ж),
— экспериментально, причем в этом случае ее реальный
(говорят теоретический) аналог у = f(x) a priori неизвестен.
Функцию φ(χ € [a,b]) используют в расчетах на χ £ [c,d\ и
— называют интерполянтом, если [с, d] С [а, Ь];
— называют жстраполлнтой, если [с, d] D [а, 6].
Основными вопросами при выборе приближения являются:
— что выбрать в качестве приближения (аппроксиманты) φ(χ),
— как количественно оценить качество приближения φ(χ).
Традиционно выбором является
η
J=0
где (j)j(x) — линейно-независимые функции, а постоянные Cj =
= const:j = 0,1,2...,η определяют решением неоднородной
СЛАУ
η
Vk = УЫ) = Σ Cj<l>i{xk), (6.0.2)
j=o

94
если определитель матрицы коэффициентов (6.0.1) не равен нулю:
Φθ(χθ) Φΐ(χθ) ··· Φη(χθ)
Φο{χι) Φι{χι) ··· Φη{χι)
φ 0. (6.0.3)
Φθ(χη) Φΐ{χη) ··· Φη(χη)
В качестве ф(х) используют
— степенные функции, φ^ (χ) — χ*;
— тригонометрические функции, ф^х) — {cos jx,sinjx} и др.
В качестве меры качества φ(χ) используют
— понятие близости φ(χ) и у(хк), если последние совпадают
на множестве узлов сеточной функции, Vfc —» \(f{xk) ~ у{хк)\ = 0,
— условие минимальности их среднеквадратичного
приближение Σ* Ых) ~ УЫ)? -* min.
Познакомимся с простейшими методами приближения сеточной
функции, заданной на равномерной и неравномерной сетках.
6.1. Приближение (полином) Ньютона
На равномерной Xk+i — Хк = &xk = h — const сетке {ж*. : к —
= 0,1,2,...,п} для аппроксимации используют полином
Ньютона, который пррст при вычислениях и строится, используя
следующий алгоритм.
Теорема 6.1.1. Для η + 1 узлов Xk+i — хк + h — хо + {к+
+l)h сеточной функции у(хк '· к = 0,п) с h — const приближение
(полином) Ньютона (например, в правых разностях) имеет вид
fc + 1 А (к + 1)к а2
Ук+ι = Уо + —у-&уо + —2\—Δ Уо+
+ (fc + 1)3fc,(fc"1)A4 + - + Afc+V (6-1.1)
Доказательство (метод полной математической индукции). С
учетом таблицы разностей

95
ж0
Х\
Х2
хз
ХА
2/о
2/1
2/2
Уз
2/4
Δνη
Ayi
Δ2/2
Δί/з
Δί/4
Δ22/ο
Δ22/!
Δ22/2
Δ22/3
Δ32/0
Δ32/ι
Δ32/2
Δ42/ο
Δ4|/ι
Δ52/ο
где χι = χο + h, Χ2 = χι + h = χο + 2h,..., Xk+i
= χο + (к + l)h, запишем последовательно:
Хк + h
xi=x0 + h : t/i = уо + Aj/0,
a?2 = xq + 2h :y2 = yi + Ayi -yo + Ayo + A(y0 + Ay0) =
= 2/0 + 2&y0 + A2y0,
хз = xq + 3/t: уз = 2/2 + Δι/2 =
= 2/o + 2Δ2/0 + Δ22/0 + A(y0 + 2Ay0 + A2y0) =
= y0 + ЗЛуо + 3Δ22/ο + Δ32/ο,... и т. д.
Пусть
кА к(к-1)л2 fc(fc-l)(fc-2)A3
2/fc = № + γ,Δ^ο + V 2, ;Δ22/0 + - ^ -Δ32/0 + · · ·
fc!Afe
Докажем,что
fc + 1
(fc + l)fcA2
2/fe+i = 2/0 + ~3ΡΔ^> + 2!
(fc + l)fc(fc-l).3
А2у0+
+
3!
д^0 + ... +
(fc + 1)!
(fc + 1)!
Δ*+12/ο,

96
Для этого запишем
к к(к — 1)
Ук+ι = Ук + &Ук =Уо + γϊΔ2/ο + —2! &2у0+
fc(fe-l)(A-2) ., *.' fc
+ -* £ '-А3уо + ■■■ + йАУ°+
. ( к . fc(fc-l) л2
+ Δ ί j/o + jjAyo + 2, yA2j/o+
fc(fe-l)(fc-2)A3 fc»fc \
+ -* ^ L^y0 + ■■■ + ^Akyo) =
Δ2ί/ο+
fc + 1 A
yo + —η—Δι/ο +
fc(fc-l) _fc
2! +Ϊ!
+
fe(fc-l)(fc-2) fe(fe-l)
3!
■· +
2!
2/0 · Уо +
A3y0 + · · ·
fc(fc - 1)... 1 . *., fc + lA
-* гт- Δ*+1?/ο · i/o + -^τ-Δι/ο+
+
fe! "w *" ' 1!
(fe + l)feA, j (fc + l)fc(fe-l).3
2!
-Δ^0 +
3!
А'и, + ...
(fe + 1)!
fe + 1
+ λ . /wAfc+1yo = 2/o + '-Чт^^Уо + ν- ■ -'--д*уо+
(fc + l)feA,
(fe + 1)
2!
что соответствует утверждению теоремы. ■
Результат теоремы буквально означает, что для известных
яо, Уо = 2/(жо), Ау0 = Ау(я?о), · · ·, Afe+1</o = Ак+1у{х0)
можно использовать 2/^_|_ι(6.1.1) в качестве приближения у(хк)'·
1!
</>(a;fc+1) = yo + —π—Δι/ο + v" ^,'y"Aayo+
+ №ϋΙν+...+4«,„.
fc + 1
£fc+l - Жо

97
причем в качестве хо, в принципе, может быть любой узел сетки.
Полезно переписать (6.1.1) для χ fi {xk : к = 0, η}.
Запишем Xk — xq + kh в виде χ — Xq + uh -ϊ и = (χ — xo)/h.
Тогда (6.1.1) принимает вид
+ ^-1К,-2)д3^ +
■■4tt(tt"1)"H(tt"t + 1)AV (6.1.2)
ПРИМЕР 6.1.1. Найти приближение сеточной функции y(xi),
заданной таблицей, а также у(х = 2,5), з/(ж = —0,5)
соответственно для хо = 1 и хо = 2.
Zfc
Ук
0
1
1
4
2
15
3
40
4
85
РЕШЕНИЕ. Строим таблицу разностей (правых) для заданной
уЫУ-
Хк
0
1
ί 2
3
4
Ук
1
4
15
40
85
Ау
3
11
25
45
А2у
8
14
20
А3у
6
6
6
Д"4у
0
0
0
Сетка равномерная: Axk = Л = 1. При ж0 = 0 -> и =
= (ж — жо)/1 = ^· Далее по таблице разностей вычисляем полином

98
Ньютона (6.1.1):
^> = 1+ϋ3+τ8+χ(*~3!(χ~2)6 =
= 1 + Зх 4- 4х(ж - 1) + х{х - 1)(ж - 2) =
= 1 + Зх + 4ж2 - 4ж + х3 - Зж2 + 2ж = 1 + ж + ж2 + х3 =
= (1 + х)(1 + х2).
(А)
Для полинома ψ(χ)(Α) вычисляем:
7 29
φ(χ = 2,5) = (1 + 2,5)(1 + 6,25) = - — = 25,375,
ψ{χ = -0,5) = (1 -0,5)(1 + 0,25) = \\ = 0,625.
С другой стороны, для ψ(χ = 2,5)|ro=i —»· и = (ж — χο)//ι = (2,5—
-0)/1 = 2,5:
ф = 2,5)U0=1 = 1 + ^3 + 2'5(2^"1)8+
+ 2,5(2,5^)(2,5-2)6 = 1 + 7|5 + 53+
+ У^= 23,5 +1,875 = 25,375.
Для у?(ж = -0,5)|Жо=1 -► и = (х - аг0)/Л = (-0,5 - 0)/1 = -0,5:
+ (-0,5K-l,S)(-2,6)6 = 1_15 + 3_
-yi-M-^-o,·».

99
Для вычисления у(х — 2,5),у(х — —0,5) для xq = 2 следует
достроить таблицу разностей до Д3у включительно. Вычисляем
аналогично для φ(χ = 2,5)|а;0=2 -* и — {х — xo)/h = (2,5 — 2)/1 =
= 0,5:
Φ = 2,б)|.ова = 1б+М25+0,б(0^"1)20+
+ o.g(Q.g-m-oV-»)6a8lB + 12t5.2tB+
+ ^ = 25 + 0,375 = 25,375.
Δ Δ Δ
Для ν(χ = -0,5)|.о=2 -► и = (ж - яо)/Л = (-0,5 - 2)/1 = -2,5:
+ (-2,5)(-3,5)(-4,5)6=15_62|5 + 87|5_
_"» 40-^=0,625,
222 8
что соответствует полученным ранее. ■

100
ЛЕКЦИЯ 15
Приближение функций (полином Лагранжа)
6.2. Полином (интерполянт) Лагранжа
Полином Лагранжа используется для приближения сеточной
функции y(xk)i заданной на неравномерной сетке. Полином
Лагранжа φ(χ) изначально формируют в форме:
η η
ψ{χ) -^2сзФз{х) = Σ Ск{х - хо){х - xi) - - ·
... (χ - хк-г){х - Xk+i)... (ж - жп) =
= Co(x - x\){x - X2) · · · (x - Xn) + Ci(x - xo){x - x2)...
... (x - xn) Η + СаДж - x0)(x - a?i)...
... (χ - χ*-ι)(χ - я?л+1)... {χ - xn)+
+ Ск+\(х - xo)(x - χι)... (ж - xfe)(x - хк+2)...
... (χ - xn) + · · · + Сп(х - χι)(χ - х2)... (ж ~ Xn-i),
(6.2.1)
причем постоянные Cj (6.2.1) являются решениями СЛАУ:
χ = х0 -+ <£>о(хо) = 2/о = С0(х0 - хг){х0 - х2)... (хо - юя))
χ = χι -» φο{χ\) = 2/1 = Ci(a?i - £o)(xi - s2)... (χι - xn)
X = Xn-+ ψθ{Χη) = Уп = Cn(Xn - Χθ)(Χη - Xl) . . . (xn - Xn-l) J
(6.2.2)
при условии (6.0.3), принимающем для СЛАУ (6.2.2)
диагональную форму
ФоЫ 0 ... О
о
о
Φι{χι)
0
О
Фп{Хп)
= Φο{χο)Φι{χι) · · · ФоЫ) т^О,
(6.2.3)

101
что позволяет записать
с Уо
(х0 -Xl){Xb ~Ж2).
с - У1
(xi -x0){xi -ж2).
. (ж0 -хпУ
.(χι -хпУ
сп = -, w ^г—-, г
(6.2.4)
(#п - Я())(Ж„ ~~ Χύ · · · (Жп "~ #n-l)
Подставляя постоянные Cjrj = 0,1,...,η (6.2.4) в φ(χ) (6.2.1),
получаем приближение сеточной функции по Лагранжу:
(χ-χι)(χ-Χ2)...(χ-χη)
Ψ\χ) = 7Г ΖΎ7Ζ ζ\ 7Ζ ГТ2/0+
(я?о ~ #l)(#0 - Я?2) . . . («О - жп) *
|(ж ~ ж2)...(ж
(я?1 - ж2)... (а
(ж - х0)(х - χι)... (χ - a?n_i)
(а;-жр)(а? - а?2)>>> (a? ~a?n) +
(я?1 ~ Ж0)(^1 - Ж2) · · · (ffl - Я?п)
■ +
:Уп. (6.2.5)
(хп - Ж0)(ЖП ~ Xl) · . - («η ~ Яп-l) *
Замечание 6.2.1. Условие (6.2.3) нарушается
ФоЫФЛхг) · · · ФоЫ = 0 & ЩФкЫ = 0 (6.2.6)
при совпадении (неоднозначности) узлов сетки приближаемой
сеточной функции y(xk)·
ПРИМЕР 6.2.1. Построить интерполянт Лагранжа для y(xi):
Xk
Ук
0
1
2
3
2
4
3
6
РЕШЕНИЕ. Вычисляем определитель (6.2.3):
ФоЫФМФгЫФзЫ = (0 - 2)(0 - 2)(0 - 3)·
•(2 - 0)(2 - 2)(2 - 3) · (2 - 0)(2 - 2)(2 - 3) · (3 - 0)(3 - 2)(3 - 2) = 0,

102
что исключает прибижение Лагранжа для данной у(х{). Ш
ПРИМЕР 6.2.2. Построить интерполянт Лагранжа для у(х{),
заданной таблицей, и вычислить φ(2.5):
Хк
Ук
0
2
2
3
3
5
6
6
РЕШЕНИЕ. Вычисляем определитель (6.2.3)
ФйЫФМФъЫФъЫ = (0 - 2)(0 - 3)(0 - 6).
.(2 - 0)(2 - 3)(2 - 3) · (3 - 0)(3 - 2)(3 - 6) · (6 - 0)(6 - 2)(6 - 3) ф 0,
что позволяет последовательно записать φ(χ) для заданной y(xk)i
= (*-2)(s-3)(*-6) (χ - 0)(s - 3)(s - 6)
ΨΚ ' (0-2)(0-3)(0-6) (2~0)(2-3)(2-6)
(g-0)(g-2)(s-6) (i-0)(i-2)(»-3)
(3 - 0)(3 - 2)(3 - 6) (6 - 0)(6 - 2)(6 - 3)
ЗЗж3 - 174ж2 - 99ж2 + 522ж - 44ж3 + 352ж2 - 528ж + 12х2 - 96а; + 144
72
ж(11ж2 -91x4-102)
72
+ 2
и вычислить
У(2,5)^-2'5(11-6'25:.91-2'5+102)+2=
2,5(-56,75)
72
72
+ 2 = 3,97135.
Вычисляем φ(χ = 2,5), используя полином
ф = 2,5) J0.5)(-0,5-3,5)2 + (2,5)(0,5)(3,5)3+
+
-36 " ■ 8
(2,_5)(0.5)(3L5)5 + (2,5)(0,5)(-0,5)
6 =
(-9) " ' 72
= -0,0486 + 1,6406 + 2,4305 - 0,0521 = 3,97135.

103
6.3. Обобщение приближений по Ньютону и Лагранжу
Анализ приближений Ньютона и Лагранжа позволяет
Следствие 6.3.1. В соответствии с доказательством
теоремы 6.1.1 наличие Δу о возможно, если в сеточной функции
у{хк) сущестует [к + 1)-й узел (см. таблицу разностей),
причем полиномы Ньютона и Лагранжа'с учетом данных сеточной
функции у(хк · к = 0,1,... , п) после преобразований принимают
структуру приближения (6.0.1):
η
<р{х)=^Саф;(х), (6.3.1)
J=0
содержащее Cj'.j = 0,1,..., (η + 1) неизвестных.
Постоянные Cj можно вычислить методами линейной алгебры
по данным сеточной функции y(xk · к = 0,1,... ,п), если
выполнено условие (6.0.3):
ίψ(χο) = Σ?=ο СзФзЫ) = 2/0,
<Ρ(*ι) = Σ?=ο СзФзЫ = Уъ (б 3 2)
Ч>Ы = Σ?=ο СзФз(хп) = Уп,
что позволяет построить для данной у(хк) канонический полином,
который проходит через все узлы сетки, т. е.
Vfc-> Mxk) - у{хк)\ = 0.
Следовательно необходим анализ
— адекватности данных сеточной функции у(хк) и модели φ(χ)
эксперименту: \φ(χ) — /(ж)| -> 0, а также

104
— погрешности приближения.
Для заданной сеточной функции алгоритмы построения
глобальных приближений по Ньютону и Лагранжу не содержали
рычагов управления структурой приближения в степенной или
тригонометрической формах. Сплайн-приближение для заданной
сеточной функции уже не глобально и допускает управление
локальной структурой, однако локальное (кубическое или линейное)
приближение единственно.
Такая "жесткость" приближения φ(χ) заданной сеточной
функции по существу исключает управление приближением, однако
позволяет оценить качество приближения φ(χ) и опытные
данные сеточной функции у(х{) реальной физической модели f(x):
\<p{xk) ~ У{хк)\ = \<p{xk) - У{хк) + f{x) ~ f{x)\ ζ
^\<p(x)-f(z)\ + \f(x)-y(xk)\, (6.3.3)
причем \φ(χ) — /(ж)|, \f(x) — у(хк)\ — соответственно
погрешности аналитического и опытного представления модели. Если
теоретический аналог у = f(x) a priori известен, то в принципе
\f(x) — у{х%)\ —> 0 и оценка (6.3.3) принимает вид
Hx)-y(xi)\<Mx)-f(x)\, (6.3.4)
которая требует наилучшего в смысле (6.3.4) аналитического
приближения φ(χ)>
Анализ данных сеточной функции
Уже обсуждали невозможность построения приближений
Ньютона и Лагранжа по данным эксперимента, которые не
удовлетворяют условию однозначности узлов сеточной функции f(xi)
(6.2.3).

105
Если данные опыта удовлетворяют определению
однозначности узлов сеточной функции, то для анализа (с целью проверки
ошибок) данных сеточной функции изначально используем
таблицу разностей. Вычислительная идея такого анализа связана
со стабилизацией разностей высших порядков для неискаженной
ошибками сеточной функции. Поясним сказанное, используя
ПРИМЕР 6.3.1. Функция у = f(x) = (1 + ж)(1 + х2) вместе с ее
разностями задана таблицей примера 6.1.1, в которую вкрались
ошибки:
Хк
0
1
2
3
4
5
6
7
Ук
1
4
15
L4il
85
156
259
400
Ay
3
11
1.26]
[44j
71
103
141
Δ22/
8
115 J
[18|
|_27J
32
38
aV|
\jj 1
iii 1
|n| 1
ш 1
6
A*y
[-4
IaI
|-6
Ш
РЕШЕНИЕ. Из таблиц разностей настоящего и примера 6.1.1
следует, что ошибка при χ = 3 распространяется от места
ошибки у = 41 сектором, причем для уточнения места ошибки
необходим дополнительный эксперимент, увеличивающий числа узлов
сеточной функции у(хк) до 8-10, чтобы добиться в
неповрежденной части таблицы стабилизации высоких разностей. В данной
таблице наблюдается стабилизация А3у разностей: с этого
столбца начинаем восстановление таблицы. После исправления ошибок
таблица принимает вид

106
Xk
0
1
2
3
4
5
6
7
8
1 9
Ук
1
4
15
40
85
156
259
400
585
820
Ay
3
11
25
45
71
103
141
185
245
Δ22/
8
14
20
26
32
38
44
60
A3y
6
6
6
6
6
6
6
A4y
0
0
0
о
0
0
Восстановленная таблица позволяет сделать некоторые выводы
относительно вида восстанавливаемой функции: с учетом \/А4у =
= 0 можно утверждать, что наивысший порядок
восстанавливаемого полинома равен трем. ■
Нарушение стабильности при ошибках опыта в значениях у(хк)
происходит и в таблице разностей для неравномерной сетки. Для
уточнения возможного места ошибки требуется серьезное
увеличение числа узлов сеточной функции (в теории вероятностей в
этом случае говорят об увеличении размера выборки \(у £ Y)\)·
Возможны иные приемы проверки адекватности
(достоверности) результатов эксперимента представленных сеточной
функцией у{хк)ч например, используя среднеквадратичные приближения
или приближения сплайнами.

107
ЛЕКЦИЯ 16
Приближение функций (сплайн-интерполяция)
6.4. Приближение сплайнами
Для приближения по Ньютону и Лагранжу характерным
является сохранение структуры многочлена на всем промежутке
интерполяции (говорят интерполирующий многочлен глобален).
Иногда погрешность приближения \φ{χ) — f{x)\ улучшается за
счет выбора разных многочленов на разных шагах
интерполирования. В этом случае говорят об представлении табличной
функции кусочными (локальными) интерполяционными
многочленами (сплайн-приближение). В основе сплайн-приближения φ{χ)
функции /(ж), заданной сеточной функцией ук{хк), лежит
понятие сплайна (spline — планка, англ.): если упругую линейку,
поставленную на ребро, зафиксировать в точках сеточной функции
yi(xi), то полученная кривая является соответствующим сплайн-
приближением, причем из механики известно, что изгиб упругого
стержня, описываемый в равновесии дифференциальным уравне-
нием 4-го порядка φ [χ) = 0, имеет решение:
φ(χ) = с0 + cix + с2х2 + с3ж3, (6.4.1)
которое не испытывает разрывы не только самой функции, но и
ее первой и второй производных.
Итак, сравнительно с приближением Ньютона или Лагранжа
для приближения сплайном характерно:
— гладкость интерполирующей кривой,
— замена глобального интерполирующего полинома
кусочными (локальными) интерполирующими полиномами,
— совпадение кривизны локальных интерполирующих
полиномов в точках таблицы Уг{хг)·
С учетом механической аналогии (6.4.1) удобно для
последующих вычислений переписать локальный интерполирующий на

108
χ € [χ^,χ^4-ι] многочлен'to виде
4>k{x) = о>к +Ък(х-хк) + ск{х - Xk)2 + dk{x - Хк)3 (к = 0,n- 1),
(6.4.2)
постоянные которого определяются:
— условиями гладкости локальных сплайнов
<ff{xk - 0) = if/(zk + 0) (fc = МП),
<//'(** - 0) = у>"(я* + 0) (fc = МГ=Т),
т. е. значения первой и второй производных локальных
многочленов слева и справа от внутренних точек сетки хк €]хо,хп[
совпадают;
— условием свободных концов локальных сплайнов
φ"(χο) = φ"{χη) = 0. (6.4.4)
Совокупность локальных многочленов (6.4.2) с условиями (6.4,3),
(6.4.4) называют кубическим сплайном сеточной функции Ух(хг).
Выясним условия существования и структуру сплайнов при
последовательном увеличении числа точек fc сеточной функции
/(zfe:fc = 0,l,2,...).
Случай к = 0 Л
Сеточную функцию изначально приближаем одним сплайном
φο{χ) = а0 + Ь0(х - х0) + с0(х - хо)2 + d0(x - хо)3 (к = 0,п- 1),
(6.4.5)
поэтому условия гладкости (6.4.3) отсутствуют (условие
гладкости существует на границе двух сплайнов), а постоянные сплайна
формально определяются значениями сеточной функции
¥>о(яо) = ао = 2/о,
¥>ο(αι) = «о + Ьо(ж1 ~ жо) + со (я ι - z0)2 + do (я? ι - хо)3 = 2/1

109
и условием свободных концов (6.4.4), которое для к = 0,1
принимает вид
φ'1{χ) = 2ск + 6dk{x -хк) = 0(к = 0,1) -> с0 = d0 = 0,
что позволяет переписать (6.4.5) в виде
¥>о(яо) = ао = 2/о, ¥>o(a>i) = ао + b0{xi - хо) = Уъ
пригодном для вычисления постоянных линейного сплайна
φο{χ) = «о + Ьо{х - хо) {к = 0,η - 1).
Случай fc = 0,1,2
Такая сеточная функция приближается двумя сплайнами
φο{χ) =а0 + Ь0(х - х0) + с0{х - хо)2 + d0{x - хо)3,
ψι{χ) = αϊ + b\(x - χι) + ci(a; - χι)2 + d\{x - xif,
для существования которых формально необходимо вычислить
восемь параметров: ао)Ьо,с<ь^о>аъЬъсь^1· С учетом заданной
сеточной функции f(xk : fc = 0,1,2) выпишем восемь
уравнений для вычисления восьми неизвестных параметров локальных
сплайнов:
выражения для локальных сплайнов:
1)φ0{χ = χο) =ао = г/0,
2) ¥>0(ж = Βι) = «о + Ьо(«1 - я?о) + со{х\ - хо)2 + d0{xi - хо)3 =
= Уь
3)^ι(χ = χι) = αι =2/ι,
4)</?ι(χ = ж2) =«ι + ЫЖ2 -»i) + ci(x2 -ζχ)2 + di(x2 ~^i)3 =
= 2/25

no
условия гладкости локальных сплайнов:
5) φ'Ό(χ = χχ) = b0 + 2c0(xi - xo) + 3d0{xi - x0)2 = φ[(χ = χι) =
6) <Pq{x - x\) = 2c0(xi - x0) + 6d0(xi - xq) = Ψι(χ = x\) = 2ci\
условия свободных концов локальных сплайнов:
7)φ'ν(χ = хо) = 2с0 = О,
8) ψι(χ = х2) — 2с\ + 6di(x2 - χι) = 0.
Полученные восемь уравнений рассматриваем как СЛАУ
относительно неизвестных параметров ao,bojCojdo,ai,bi,ci.,c?i, решение
которых существует, если СЛАУ совместна (сопуствующая
проблема однородности или неоднородности СЛАУ зависит в каждом
конкретном случае от заданной сеточной функции f(xk), что
затрудняет в общем случае анализ совместности СЛАУ).
Необходимо уменьшить число уравнений СЛАУ относительно
неизвестных параметров локальных сплайнов, так как условия 1),
4) и 7) тривиальны и позволяют вычислить
1)ао = Уо, 3) ai =2/i, 7) с0 = 0, (6.4.6)
а оставшаяся СЛАУ состоит из пяти уравнений:
2) φ0(χ = χι) = а0 + b0(xi - х0) + c0(xi - х0)2 + d0(xi - х0)3 =
= Уи
4) ψι(χ = х2) = ai+ bi(x2 - χι) + ci(x2 - χι)2 + di(x2 - Χι)3 =
= 2/2,
5) φ'0(χ = χι) = b0 + 2c0{xi - xo) + 3d0{xi - xo)2 = ψ\{χ = x\) =
6) φ'ό{χ = χι) = 2c0{xi - xo) + 6d0{xi - xo) = φ"{χ = xi) = 2ci,
8) φι(χ = x2) = 2ci + 6di(x2 - χι) = 0,
(6.4.7)

Ill
решение которых с учетом (6.4.6) позволяет построить
минимальный кубический сплайн для сеточной функции f(xk : к = 0,1,2).
ПРИМЕР 6.4.1. Построить минимальный кубический сплайн
для сеточной функции y(xi), заданной таблицей и вычислить
р(2.5).
Хк
Ук
?,
-20
0
6
2
0
РЕШЕНИЕ. Для заданной f(xk) вычисляем по (6.4.6)
1) \а0 = -20 , 3) αχ = 6 , 7) \с0 = 0
и последовательно преобразуем систему (6.4.7) в общем виде:
2) - 20 + Ьо(0 + 2) + d0(0 + 2)3 = 6 -> | Ь0 + 4d0 = 1з],
4) 6 + 6ι(2 - 0) + ci(2 - О)2 4· di(2 - О)3 = 0 -> 13 + h + 2сг + 4di = 0
5) Ь0 + 3d0(0 + 2)2 = <^Ί(х = χι) = bi, -> | b0 + Ш0 =~bT
6)6d0(0 4- 2) = 2ci -> I 6d0 = ci
8)2ci + 6di(2 - 0) = 0 -> 1 ci + 6di = 0
ci = 6do = -6di -H do = -d\
Тогда из 5) с учетом 2) следует
bi = b0 + Ш0 = 13 - 4d0 + Ш0 = 13 + 8d0 = h
что позволяет по 4) вычислить 3 + 13 + 8d0 + 4d0 + 4d0 = 0 -»
, следовательно по 2), 6) и 5) соответственно
-»
do
-1 = -άΛ
вычисляем
Ьо = 13 - 4d0 = 17 , ci = -6 , bi = 17 + 12(-1) = 5
что завершает вычисление постоянных:
ао = -20, Ьо = 17, со = 0, do = -1, αϊ = 6, 6χ = 5, ci = -6, di = 1

112
локальных сплайнов
φ0{χ£ [-2,0]) = ж3-6ж2 + 5ж + 6, (рг(х € [0,2]) =
= -ж3 - 6ж2 + Ъх + 6. ■
Представляет интерес хотя бы на примерах сравнить
геометрию сплайнов при наращивании числа точек сеточной функции
или упрощении структуры локальных сплайнов.
ПРИМЕР 6.4.2. Построить кубический сплайн для сеточной
функции y(xk : к = 0,1,2,3), включающей сеточную функцию
примера 6.4.1:
Хк
Ук
-2
-20
0
6
2
0
4
10
РЕШЕНИЕ. Вычисляем параметры локальных сплайнов
<Ро(х £ [хо, χι]) = °>о + Ьо{х - хо) + с0(х - хо)2 + d0{x - ж0)3,
ψι(χ е [a?i,a?2]) = αϊ + bi(x — xi) + ci(x - χι)2 + di(x -χι)3,
ψ2{Χ £ [X2, хъ\) = 0>2 + Мж - X2) + C2{x ~ X2)2 + d>2(x - Ж2)3
с учетом f(xk : к = 0,1,2,3):
α0 = -20
αϊ
α2 = 0 , со = 0
φ0(χ = Χι = 0) = -20 + bo(0 + 2) + d0(0 + 2)3 = уг
= 6->|
b0 + 4d0 = 13
φι(χ = Ж2 = 2) = 6 + 6ι(2 - 0) + ci(2 - О)2 + di(2 - О)3 = y2 =
= 0
3 + δι + 2ci + 4di = 0
^2(ж = χ3 = 4) = 0 + Ь2(4 ~ 2) + с2(4 - 2)2 + d2(4 - 2)3 = 10
3 — ΊΠ —
2/3
b2 + 2c2 + 4d2 = 5
φ'0(χ = «ι = 0) = δ0 + 3d0(0 + 2)2 = ^(ж = хг = 0) =
= fti -Н Ьо + 12do = bi

113
φ[(χ = x2=z2) = bl + 2ci(2 - 0) + 3d0{2 - О)2 = φ'2(χ = x2 = 2)
Ь2 -> bi + 4ci + 12di = δ2
<р{((я = Ж1 = 0) = 6d0(0 + 2) = <ρί (ж = x2 = 2) =
= 2ci -> I
6cfo = c\
φ'{(χ = χ2 = 2) = 2ci + 6di(2 - 0) = ^{x = x2 = 2)
= 2C2 -H cl + 3^1 = c2
^(ж = ж3 = 4) = c2(4 - 2) + 6d2(4 - 2) = 0 -* \c2 + 6d2 = 0
Тогда 6χ = Ьо + 12do = 13 -4d0 + 12do = 13 + 8d0 = h , что
позволяет по 4) вычислить 3 + 13 + 8do + 4cfo + 4cfo = 0 ->
do
di L следовательно по 2), 6) и 5) соответственно
вычисляем
b0 = 13 - 4d0 = 17 , ci = -6 L bi = 17 + 12(-1) = 5
что завершает вычисление постоянных:
α0 = —20, bo = 17, со = 0, do = —1, αϊ = 6, &i = 5, c\ = —6, di
локальных сплайнов
φ0{χ e [-2,0]) = χ3 - 6ж2 + 5ж + 6, ¥?ι(ж € [0,2]) =
= -ж3 - 6ж2 + Ъх + 6. ■
6.5. Линейный сплайн
Хотя физически оправданы локально-кубические сплайны, для
аппроксимации применяют также линейные сплайны с локальным
линейным многочленом вида
4>k{x) = flfe + bk(x - Xk) (fe = 0, η - 1).
(6.5.1)

114
При приближении линейным сплайном y(xi), х% 6 [а, Ъ] заменяют
ломаной (кусочно-линейной), т. е. производят линеаризацию
сеточной функции y(xi), поэтому отсутствуют условия гладкости
(6.4.3) (условие гладкости существует на границе двух сплайнов)
и свободных концов (6.4.4), причем постоянные сплайна
определяются значениями сеточной функции на концах отрезка
интерполирования [Xk,%k+l]:
Vk{xk) =а>к= Ук, Ч>к{хк+\) = а>к+ bo(xk+i ~ хк) = Vk+ъ
ПРИМЕР 6.5.1. Построить линейный сплайн для сеточной
функции у{хк : к = 0,1,2,3) примера 6.4.2:
Хк
Ук
-2
-20
0
6
2
0
4
10
РЕШЕНИЕ. Вычисляем параметры локальных линейных
сплайнов
<Ро(в £ fro, χι\) = αο + b0(x - (-2)) -*
ίφ0{χ = χο) = αο = 2/ο = -20]
φ0(χ = χι = 0) = -20 + 260 = 2/ι = 6 -> b0 = 13 Ι
-»
¥>о(я € [a?0,a?i]) = 6+ 13ж
ψι(χ e [a?i,a?2]) = αι +bi(a; - 0) -*
f ¥>i(a? = ^i) =αι =yi = l
->
-»
i(/?i(x = ^i) =ai =yi = 6]
φι(χ = ж2 = 2) = 6 + 2bi = 2/2 = 0 -* bi = -3j
¥>i(a e [a?i,a?2]) =6 — За: ;

115
Ψ2{χ € [»2,жз]) = а2 + b2{x - 2) -*
<р2(я = ж2) =а2 = у2 = 0]
хз = 4) = 0 + 2Ъ2 = уз = 10 -* Ь2 = 51
->
\</>2(ж
-»
V?2(# £ [#2,#з]) = -10 + Ъх
Для линеаризации нелинейных соотношений используют
аппроксимацию линейным сплайном. ■

116
ЛЕКЦИЯ 17
Точечное квадратичное приближение функций
Математическая статистика в результате массового
эксперимента объемом \W*\ и исключения повторяющихся элементов
W* = {w* : \w*\ ^ 1} наблюдаемой случайной величины w
получает генеральную совокупность W.
Определение. Генеральной совокупностью случайной
величины w называют множество всех ее различных наблюдаемых
значений (пространство элементарных событий) W — {щ : г =
= 1,п}, позволяющее по результатам массового опыта
восстановить ее истинные (теоретические) закономерности.
При очень большом \W\ переходят к заданию выборки на
отрезке (в статистике говорят на интервале), например: {vb\ — u^),
(w2 - w3),..., (щ - щ+ι),...,(%- t&fe+i), причем точки,
лежащие на границе отрезка, учитываются в обоих граничных
отрезках. Границы интервала могут не совпадать с точками W (на
практике большие выборки делят на 15-20 отрезков).
Полигоном частот дискретной выборки называют
построенную в системе координат 0г&,-|1и?| ломанную, соединяющую точки:
(ώχ, \wl|),..., (u)k, |wjl|). Гистограммой частот непрерывной
выборки, построенную в координатах 0u)i\w*\/(u)i — ΐϊ^+ι), называют
множество прямоугольников со сторонами: (w\ — t&2), \wl\/(wi —
-ώ2);... {wk - Wfc+i), \w%\/(wk ~ tSjb+i)·
Увеличение объема выборки не единственный способ
восстановить теоретические закономерности и характеристики
наблюдаемой в опыте w. Обработку ограниченной выборки W: \W\ « \W\
для восстановления закономерностей и характеристик
наблюдаемой в массовом эксперименте объемом \W\ случайной величины
w называют выравниванием (сглаживанием) статистического
ряда W, что составляет содержание основных задач
математической статистики: статистической оценки (параметров) и
проверки правдоподобия (статистических гипотез). В задаче
статистической оценки вычисляют по ограниченной выборке W статис-

117
тические характеристики параметров случайной величины
(точечная оценка) или их доверительные интервалы (интервальная
оценка). Оценки должны удовлетворять условиям:
состоятельности, pi = pi = ρ*; несмещенности, (w) = Mm; эффективности,
{w — (w))2 = Thv —l· ruin. Выборку W называют
представительной, если статистические закономерности и характеристики W и
W совпадают (или не превышают погрешности опыта).
6.6. Метод наименьших квадратов (МНК)
Точечное квадратичное оценивание. МНК
Наблюдаемая в опыте величина yk — j(pk) + £к отличается от
истинного (теоретического) значения f(wk) на величину
случайной ошибки (флуктуации) £, которая связана с ошибками
измерения и др. Как правило, ошибки ξ = у — f(w) массового опыта
распределены по нормальному закону с параметрами: М£ = Ό =>
=» Щ = М(£-М£)2 = Щ2-(М£)2 = М£2 -* МЮ>£ = ММ£2 = Щ2,
причем значение Μξ = 0 в массовом опыте объясняется
отсутствием систематической ошибки измерения случайной величины.
Если в опыте присутствует систематическая ошибка, М£ φ 0, то
ее устраняют и восстанавливают М£ = 0.
Ограниченность объема выборки W может исказить
(например за счет ошибок измерения) статистическое поведение
генеральной совокупности W", что нарушит условие несмещенности,
(w) = Мг5, и возникнет задача сглаживания случайностей W
приближением tp(w,..., щ а,..., d), где w,..., и — переменные,
а,..., d — параметры.
Если параметры а,..., d приближения y?(iu,..., щ а,..., d)
оценивают по выборке W сеточной функцией уk = y(wk, · · ·, Щ · {к =
= 1, п)), а в качестве критерия согласия принято условие min
суммы квадратов невязок (ошибок, EIRror)
η
™£ = Σ(ψ(™><, .··,«*)- Ук? -»· min, (6.6.1)
fc=l

118
то такое приближение называют точечным среднеквадратичным.
Точечное квадратичное оценивание рассмотрим на примере
сглаживания выборки yk = y{wk) приближением φ(ιν) = а + bw:
Ιξ -> min <
= 2^2(a + bwk - ун) = 0,
flERf ftA,
—— = 2 2^(a + b^fc - i/jb)wjb = 0.
(6.6.2)
k=l
С учетом определения статистического среднего
Σ к уг^п sr-лп о
η
η
η
η
(6.6.3)
перепишем (6.6.2) для дальнейших преобразований в виде:
а<ю> + Ь(«;2> = (toy) }*\Ы №) )\Ь)~~\ <«*> У ' 1 '
VF^'W = I , . , 9> I — симметрична, положительно определен-
V И <«° > /
ная матрица. Решаем методом Гаусса:
а - {у} ~ b{w)
(6.6.4)
{y)(w} - b(w)2 + b{w2) = {tot/) =Ф-
6= -s
(wy) - (y){iv)
(w2) - (w)*
„ /,λ (g)(w)-(wy)/...\ _,
a- (y)(w2) ~ (w)(wy)
(го2) - (to)2
(6.6.5)

119
Критерии согласия приближения с опытом.
В качестве критерия согласия (меры адекватности)
приближения с опытом (выборкой), следовательно упорядоченности разных
приближений между собой, служит
— число обусловленности cond(WTW) (6.6.4);
— ERf квадратов ошибок £к = Ук — f{wk) (так называемый
анализ ошибок). Для f(w) = а + bw с учетом (6.6.3) ER£ (6.6.1)
принимает вид
Ш£ = ^{а + Ъюк-ук)2 =
k=l
п(а2 + b2(w2) + {у2) + 2ab(w) - 2а(у) - 2Ъ{шу)), (6.6.6)
причем с учетом М(£) = 0 и обработки эксперимента [15] для
теоретической кривой очевидно
(ЕК£)/п = (М(0)2 = Ю>е
ПРИМЕР 6.6.1. Сгладить выборку, используя приближения
φ(χ) — а + Ьх, φ(χ) = с + d/x, φ(χ) = е + f/x2, и оценить их
качество (упорядоченность):
Хк
УЫ) = Ук
1
105
2
30
4
11,25
5
9
10
6
РЕШЕНИЕ. Вычисляем коэффициенты а,Ь:\ φ(χ) = a + bx
<ar >= (1 + 2 + 4 + 5 + 10)/5 = 4,4;
<у >= (105 + 30 + 11,25 + 9 + 6)/5 = 32,25;
<у2 >= (1052 + 302 + 11,252 + 92 + 62)/5 = 2433,7125;
<х2 >= (1 + 4 + 16 + 25 + Ю0)/5 = 146/5 = 29,2;
<ху >= (1 · 105 + 2 · 30 + 4 · И, 25 + 5 · 9 + Ю · 6) /5 = 63;

120
32,25-29,2-4,4-63 „„ „nr
a= 29,2-4,4* ~67'5305;
63-4,4-32,25^,
b- 29,2-4,42 8'°183·
Вычисляем с,d:\ ψ(χ) = с + d/ж —> Χ = 1/χ : φ(Χ) = с + dX
<Χ>=< 1/χ >= (1 + 1/2 + 1/4 + 1/5 + 1/10)/5 = 0,41,
< Χ2 > =< 1/χ2 >= (1 + 1/4 + 1/16 + 1/25 + 1/100)/5 = 0,2725,
< Xy >= (1 · 105 + 1/2 · 30 + 1/4 · 11,25 + 1/5 · 9 + 1/10 · 6)/5 =
= 25,0425;
32,25-0,2725-0,41-25,0425
C = 0,2725-0,412 * ~14'1695;
_ 25,0425-0,41-32,25 „
d~ 0,2725-0,412 ~Ш>2184·
Вычисляем е, /: φ(χ) = е + f/x2 -> Χ — 1/χ2 : φ(Χ) = e + fX
<Χ> = 0,2725;
< Xy > = (1 · 105 + 1/4 · 30 + 1/16 · 11,25 4-1/25 · 9+
+ 1/100 ·6)/5 = 22,7251;
< Χ2 > = (1 + 1/16 + 1/256 + 1/625 + 1/10000)/5 = 0,2136;
32,25 · 0,2136 - 0,2725 · 22,7251
0,2136-0,27252
= 22,7251-0,2725-32,25 _
7 0,2136-0,27252
5,
Оцениваем по Щ (6.6.1) качество (упорядоченность)
приближений:
ψ{χ) = 67,5305 - 8,0183а; чЩиО, 0001 -> ΈΆξ » 4505 -►
-» cond{WTW) « 82;

121
φ(χ) = -14,1695 + 113,2184/ж -* Μξ π 0,0002 -* ER£ « 225 -*
->cond(VTTW)«13,4;
у = f(x) = ь + Ш-+М£ = 0-*Щ = 0,cond{WTW) « 8,5 -*
—> теоретическое приближение. ■
ЗАМЕЧАНИЕ. Если порядок η аппроксимирующего
многочлена φ(χ) = X^_0Cj#J и число η сеточной функции у(хк)· (к =
0, п) при ajj+i — ajj = Axi = hi = h = const совпадают, то точечное
квадратичное приближение совпадает с интерполянтой Ньютона.
При очень большой выборке \W\ переходят к заданию
выборки на отрезке (в статистике говорят на интервале), например,
{wi - ги2), {w2 - ги3),..., {w{ - wi+i), ...,(wk- wjh-i), причем
точки, лежащие на границе отрезка, учитываются в обоих граничных
отрезках. Границы интервала могут не совпадать с точками W
(на практике большие выборки делят на 15-20 отрезков).
ПРИМЕР 6.6.2. Вычислить выборочные характеристики (М£,
Щ) опыта, представленного таблицей:
W»,Wt+l
Wi
ы*)
\Wi\/W
Wi
15,25
20
20
0,04
-40
25,35
30
70
0,14
-30
35,45
40
60
0,12
-20
45,55
50
40
0,08
-10
55,65
60
50
ОД
0
65,75
70
60
0,12
10
75,85
80
70
0,14
20
85,95
90
50
ОД
30
95,105
100
80
0,16
40
*)|ΐϋί| — число точек на отрезке [wi,iUt+i], & = 9 — число отрезков
статистического ряда.
РЕШЕНИЕ. Для упрощения вычислений статистический ряд
(упорядоченная выборка) целесообразно преобразовать:
— вычисляем среднее значение интервала w% = (wi + Wi+1)/2;
— вычисляем W — Σί~ |u>t| = 500;

122
— учитывая постоянный шаг выборки, сдвигаем среднее гйг-,
Wi = Wi + const = W{ — 60. Оцениваем статистические
характеристики выборки:
MtS = (-40) · 0,04 + (-30) · 0,14 + (-20) · 0,12 + (-10) · 0,08+
+ 10 · 0,12 + 20 · 0,14 + 30 · 0,1 + 40 · 0,16 =
= 4,4 =* MIU = 64,4,
Ш; = 402 · 0,04 + 302 · 0,14 + 202 · 0,12 + 102 · 0,08 + 102 · 0,12+
+ 202 · 0,14 + 302 · 0,1 + 402 · 0,16 = 660,
Ш = 660-4,42 = 640,64.
Далее производим точечное сглаживание средних интервальных
оценок МНК.

123
ЛЕКЦИЯ 18
Интегральное квадратичное приближение
6.7. Интегральное квадратичное приближение. МНК
Рассматривается задача о приближении функции f(x(£ [α, δ]))
обобщенным многочленом
j=0
В качестве меры такого приближения используют интегральное
среднеквадратичное отклонение (аналог (6.6.1)):
1= [ [φ(χ) - f(x)}2dX> (6.7.1)
J a
причем наилучший обобщенный полином тот, для которого
I=»min. (6.7.2)
Условие (6.7.1) позволяет записать СЛАУ относительно йеизвест-
ных (f)j : j = 0, η
= / Ш ~ /(*)№<** = 0, (j = 0, η), (6.7.3)
J a
причем решение СЛАУ (6.7.3) единственно, если
<М#о) Φΐ (χθ) ··· Фп{щ)
Φο{χι) Φι{χι) .·· Φπ{χι)
Φθ{Χη) Φΐ{Χη)
Φη{Χη)

124
ПРИМЕР 6.7.1. Аппроксимировать квадратично функцию
f(x) = Зж на χ £ [—1,1] полиномом первой степени.
РЕШЕНИЕ. Ищем решение в форме φ(χ) = со + с\х, причем
Фо{х) = 1,01 (я) = ж на [-1,1]:
0о(#о) 01 (яо)
0ο(#ι) 0ι(^ι)
^0.
Параметры полинома cq,ci наилучшего приближения определяем
из соотношений (6.7.1), (6.7.2):
Π = / [^(ж) - /(z)]2cte -* min -* I =
./α
= / (co + ciE - 3*)2ώ -> min.
Минимизируя I = J_1{co + c\x — 3x)2dx по параметрам со и с\\
ОТГ /»1
(i) —- = 2 / (с0 + С1Ж-Зж)^ = 0,
cte0 J-i
(ii) ^ = 2 / (c0 + ci« - Зж)^ж = О,
vci У-i
что позволяет выразить параметры а, Ь приближения
, ι
(i)
0со
СоЖ +
:2 3*
1пЗ
J -1
3 1
2со~ьз + зьз
= 0ч
со =
31пЗ
1,214
(ϋ)
5ci
сож2 ciar31
2ci
3
£
ι1 /-1
J_i J-i
3
3* recto = О -)■
ci
3 Z"1
:2У-1
3xardx

125
Последний интеграл интегрируется по частям
/ x3xdx = I ^- ■ 3х\n3dx = т^- f xd(3x) =
J-ι J^lnx 1пЗУ_! V '
=άΗί'-/,1Η=[("8')1·-Εΐ]
1
1η3
s + i
3 1
+
3 1η3 31η3
10
ЗЬЗ 31η2 3
2 о'
что позволяет вычислить
с\
1пЗ 1п23
4,551-3,314 = 1,237
и записать интегральное квадратичное приближение функции
f(x) = 3х на [—1,1] в виде
<р(х) = с0 + cix = 1,214 + 1,237ж (ж е [-1,1])
P. S. Для сравнения приведем интегральное квадратичное
приближение функции у(х) = 3х на [—1,1] полиномом второй степени
φ(χ) = с0 + ах + с2ж2 = 0,8013 + 1,5344ж + 0,9806ж2:
0,8013,
1 4
С0~~шз + з1^з
5 4
С2 =
1пЗ Ь23
65 170
91n3 91n23 31nd3
40
+ ^r^rz ~ 0,9806.
Возможно упростить нахождение многочлена наилучшего
приближения, если система функций (f>j(x), (j — 0,п) ортогональна.
Определение 6.7.1. Систему функций
Φί{χ), (j = 0,n)
(6.7.4)

126
называют ортогональной на [а,Ь], если
(Φπι,Φι) = / φηι{χ)φι{χ)άχ = О, тф1,
J a
причем
\\фт{х)\\ = [ Ф\х)ах (6.7.5)
J a
называют нормой функции фт(х) на [а,Ь].
Очевидно У<^/(аг), у которых ||<^(ж)|| φ 0, допускают
нормировку:
Ш*)1Г
Определение 6.7.2. Если
Vie0^:||^(x)|| = l,
то систему функций 0j(#), (j = 0,га) называют ортонормирован-
ной, причем для ортономированной системы функций
выполняется условие
гЪ
φπι{χ)Φΐ{χ)άΧ = 5m/,
/
J a
где 5mi — постоянная Кронекера
1, т = /,
тф1.
^,={о!
Целесообразно отметить, что в примере 6.7.1 система функций
фо = 1,0ι = х линейного приближения ψ(χ) = со + с\х на [—1,1]
ортогональна
(фо = 1,01 = ж) = / 0o(x)0i(x)dx = / 1 · жйж = 0,
что исключает возможное упрощение нахождения полинома
наилучшего приближения примера.

127
ЛЕКЦИЯ 19
Численное дифференцирование
6.8. Численное дифференцирование
Задача численного дифференцирования возникает при
обработке данных эксперимента в форме сеточной функции {#&, уь : к —
= О, п} или вычислениях, например, для непрерывных функций,
непосредственное дифференцирование которых по тем или иным
причинам затруднено.
Аппроксимация сеточных функций
Целесообразно отметить, что для непрерывной функции одного
аргумента у = у(х Ε [а,Ь]) определение производной (с учетом
существования производной)
у'{х)= lim ^ (6.8.1)
у v ; Δ*->ο Αχ ν }
приводит к использованию в численных расчетах:
— правых разностных производных и правых разностей
у'(х + 0) = lim yi+1 "Ш = lim ^; (6.8.2)
Axi->0 Xi+i — Xi Axi-ϊΟ АХ{
— левых разностных производных и левых разностей
у'(х-0)= lim ViZVizL= lim ^»Ц (6.8.3)
Аж<-1->0 Xf— Xi-i Аж<-1->0 AXi-i
— центральных разностных производных и центральных
разностей
у'(х = 0)=1- lim (Wl* + W^Wziy (б.8.4)
2 Δζί,Δζί-ι-ίΟ \Χί+ι ~ Xi Xi—Xi-iJ
которые для равномерной сетки Xi+i — х{=.х{ — χ{_λ = h = const
соответственно принимают вид

128
— для правых разностей
левых разностей
г/Ог + О)**^; (6.8.5)
ν'(χ-ο)*"ψ; (6.8.6)
— центральных разностей
y(ar = 0)Ri2V Λ J-2T· (6·8·?)
Далее для численного интегрирования можно использовать в
качестве аппроксимирующего приближения полиномы Ньютона или
Лагранжа, что позволяет для полинома Ньютона
/ ч и л и(и — 1) А о
у(х) =у0 + ^Ауо + + ν 2! ;A2g/o + . ·.
t ιχ(ιχ - 1)... (и - η + 1) Λ п χ -х0
+ — -Δη</0, где и = ——-
η! Αι
записать его производную в форме
dy dy du
dx du dx
1/1 , 2u-l 2 ( 3u2-6u + 2 3
Λ Vl!A2/0 + ~2Г У0 + 3! W0 + '"" J ' (6'8'8)
Аналогично
„( . _ d t( d __ dy du __ ( d dy\ du dy d2u
dx dx dudx \dxduj dx dudx2
(d dy\ du __ (du\ d2y
dxdujdx \dx J du2'

129
а, следовательно,
d2y(x)
dx2
ν"(χ) = Ίλ (д22/о + {и - 1)Δ32/0 + ...)· (6-8.9)
h2
Для нахождения производных функции у(х) для табличных
значений аргумента χ служат более простые формулы, следующие
из основных при χ = хо, следовательно, при и = 0:
dy(x)
dx
Аналогично
h \1\АУо~2\&2Уо+у&3Уо +
(6.8.10)
У"{х) = У(х) = ^ = U (Δ2Ι* - Δ»* + ...). (6.8.11)
ώι
dx h2
ПРИМЕР 6.8.1. Используя приближение сеточной функции
y(xi), заданной таблицей (пример 6.6.1), вычислить у'(х = 2,5)
для xq = 0, хо = 1, хо — 2 и 2/" (ж = -0,5) для xq = 3
Zfc
J/fc
0
1
1
4
2
15
3
40
4
85
РЕШЕНИЕ. Строим таблицу разностей (правых) для
заданной y(xk):
Хк
0
1
2
з
4
Ук
1
4
15
40
85
Ау
3
11
25
45
А2у
8
14
20
А3у
6
6
6
Д4у
0
0
0

130
Сетка равномерная: Лж& = h — 1.
С учетом h — 1 для жо — О вычисляем и = (х — xq)/\ — χ.
Далее по таблице разностей вычисляем у'(2,5) при xq
(7.1.8), причем и = {х~ x0)/h - (2,5 - 0)/1 = х = 2,5:
<//(ж = 2,5)1^=0 =
1 2-2,5-Г
ii3+-ir-8+
3-6,25-6-2,5 + 2,;
+ 3! 6
= 3 + 4 · 4 + 18,75 - 15 + 2 = 24,75.
Для ψ'(χ = 2,5)|Я0Я1 ->u=(i- а*)//» = (2,5 - 1)/1 = 1,5:
^'(^ = 2,5)1^=! = -
U, 2-1,5-1,,
ΪΤ11 + ^! 14+
3-2,25-6-1,5 + 2 '
+ si 6
= 11 + 14 + 6,75 - 9 + 2 = 24,75.
Для вычисления у'(х — 2,5) для xq = 2 следует достроить таблицу
разностей до А3у включительно.
Вычисляем аналогично для ψ'(χ =. 2,5)1x0=3 -> « = (x—xo)/h =
= (2,5 — 3)/1 = —0,5:
¥>'(* = 2,5)|ΙΟ=3 = 15 + ξ^25 + °'5(0^~1)20+
| 0,5(0,5-1)(-0,5-2)6
3!
113
= 15 + 12,5 - 2,5 + - - - = 25 + 0,375 = 25,375.
Δ Ζ Ζ

131
Для ψ(χ = -0,5) \х„^2 -+ и = (χ - x0)/h = (-0, б - 2)/l = -2,5:
(-2,5)(-3,5)(-4,5)
+ 3! 6 =
К 7 о Qi t;
= 15-62,5+ 87,5---- = 40- — =0,625,
что соответствует полученным ранее. ■
Аппроксимация функций на отрезке
Практический интерес представляет задача вычисления
производных непрерывной функции, например, /(ж), однако
аналитическое вычисление производных (по определению производной)
затруднено (например, вследствие сложности ее структуры). В
принципе всегда возможно
отказаться от аналитического вычисления производных
заданной функции, например /ж), и
— ограничиться численным вычислением производной
функции, например /(ж), χ £ [α,δ], построив для f(x) на [а, Ь] таблицу
Xi,Vi- i = 0,n.
Мерой качества аппроксимации φ(χ), τ. е.
|/(*)-¥>(*)!=> о
служит устойчивость значений производных при изменении
шага таблицы:
ПРИМЕР 6.8.2. Сравнить значение производной функции у =
sin(x = 0,5) с ее значением, вычисленным по определению
производной (правой) для различных значений Ах.
РЕШЕНИЕ. Итак у' = (вт(ж)); = совж^о.б = 0,8775826
(напомним sin 0,5 « 0,479425538). В таблице приведены значения
производной, вычисленные по определению для разных Лж:

132
Δχ
0,01
ΙΟ"3
10~4
ю-5
ю-6
ю-11
ю-13
ю-16
sin(a;+Aa;)—sin x
Αχ
0,875170888
0,877343307
0,877558589
0,8775801647
0,8775823
0,87758263
0,8775857
0,8744267
Итак, даже для такой гладкой функции, какой является у =
sin(x = 0,5), устойчивость значений ее производной при Δ —» 0
нарушается, что связано с погрешностью математической модели
для sin ж.
В разложении sin(x + Ах) = sin χ + Ах cos χ — (Ax)2 cos χ/2 +...
(Αχ)η sinfx + ηπ/2)
• · · -\ : указанной в таблице неустойчивости
η!
производной не наблюдается. ■

133
ЛЕКЦИЯ 20
Численное интегрирование (метод прямоугольников)
ГЛАВА 7. Численное интегрирование
7.1. Первообразная и определенный интеграл
В анализе рассмотрение операции интегрирования для
простоты начинают с вычисления определенного интеграла f f(x)dx.
Толкование f f(x)dx как площади под кривой f{x(E [α, b]))
связано с вычислением первообразной функции F(x)\^:
F{b)-F{a)= I f{x)dx. (7.1.1)
J a
Действительно, по теореме о среднем можно записать
рЪ /»Ь
/ f{x)dx =< f{x) > dx=< /(я?) > (6 - а), (7.1.2)
J a J a
что в прямоугольной декартовой системе координат Оху
соответствует площади прямоугольника со сторонами < f(x) > и (Ь — а).
Поэтому при численном вычислении интегралов используют
понятия: квадратура, квадратурные формулы и другие, причем слово
"квадратура" обозначает площадь.
Задача численного интегрирования возникает для
подынтегральных функций, аналитическое интегрирование которых
нецелесообразно или затруднено при обработке табличных данных
{xfe, у: к = 0,п} эксперимента.
Рассмотрим некоторые численные методы вычисления
определенных собственных интегралов Ja f(x)dx с непрерывной
подынтегральной функцией у = f(x) одного аргумента без особенностей
на [а, Ь] (несобственные интегралы рассмотрены в параграфе 7.5).

134
7,2. Метод прямоугольников
Итак, существует непрерывная на [а, Ь] подынтегральная
функция у = f(x) (рис. 7.2.1), причем интеграл fa f(x)dx
необходимо вычислить численно (хотя в принципе может существовать и
аналитическое решение). Иное представление подинтегральной
функции у = f(x) в принципе не меняет полученные
квадратурные формулы в рассматриваемых методах.
♦ /(*)
■►я
Xq Х\ Х\ Χι Χι Χ} Хз
Рис.7.2.1. Метод прямоугольников
Разобьем [а, Ь] узлами сетки {χο,χχ, * * · ,^п} на η равных
частей, (Ь — а)/п — Η, что позволяет записать квадратурные
формулы:
левых прямоугольников
S= t f(x)dx = h[f(x0) + f{Xl) + f(x2) + .. · + /(*„„!)] + Rh
J a
Ri < -^h-max\f'(a^x < b)|, (7.2.1)

135
правых прямоугольников
S= f f{x)dx = h[f{Xl) + f{x2) + f{x3) + -.. + f(xn)] + Дг,
J a
Яг < ~h-max\f'(a^x^b)\, (7.2.2)
центральных прямоугольников
S = J f{x)dx = /*[/(< #0 >) + /(< xi >) + ...
••• + f(<Xn-l >)]+Яс>
лс< ^F^2·max\f"(a<x^b)l (7·2·3)
где Ri,RT,Rc соответственно погрешности приведенных
квадратур, полученные при условии существования соответствующих
производных на [а,Ь], где
lim Ri,r,c = 0, (7.2.4)
причем
(5 - Ri) < (£ - Rc) ^(S- Rr) ^Ri^Rc> Rr.
Как показывают расчеты, условие (7.2.4) выполняется раньше,
поэтому среди "кондовых вычислителей" существует прием: при
вычислениях на ЭВМ не использовать формулы погрешностей
квадратурных формул, а проводить расчеты последовательно при
Л, h/2, /г/4,... вплоть до получения устойчивых результатов
счета прямо на квадратурах.
Замечание 7.2,1. Если пределы интегрирования совпадают с
узлами сетки, xq = а и хп = Ь, то квадратурные формулы
называют формулами замкнутого тина, а иначе — открытого типа.

136
ПРИМЕР 7.2.1. Вычислить интеграл
r2 xdx
+ хг
по формуле квадратур центральных прямоугольников с
точностью ε < 0,01.
РЕШЕНИЕ. Напомним квадратуру центральных
прямоугольников
S == / f(x)dx = h[f(< хо >) + /(< χι >) + /(< χ2 >) + ...
J a
••· + /(<Χη_! >)]+Дс,
где
<xi>=Xi+£i+1 (i = 0,l,2,...,n-l), /(х)=ж/(1 + ж2),
Rc < ^/ι2 · max \f"(a < χ < 6)|.
Вычисляем последовательно max \f"(a ^ χ ^b)\ подинтегральной
функции:
w/ \ _ 1 ж(-1)2ж _ 1-ж2
; (X) ~ ΪΤί2 + (1 + x2)2 ~ (1 + ж2)2'
f"f ϊ - ~2x (1-ж2)(-2)2ж _ 2ж3 - 6ж
; W ~ (1 + ж2)2 + (1 + x2)3 "(l + a:2)3'
.,„, , _ 2 · 3s2 - 6 (2ж3 - 6а)(-2)2я _ -6ж4 + 36ж2 - 6 _
{Х)~ (1 + χψ + {1 + χψ ~ (1 + х2)4
= 0 -> (1 + ж2) φ 0 -> ж4 - 6ж2 + 1 = 0 -»· ж2_2 =
= 3±л/8 = 3±2,83.

137
Итак χχ,2 = ±-\/5,83 > Ь = 2 не удовлетворяет условию задачи,
жз,4 = ±л/0> 17 —> жз = +\/0) 17 — корень уравнения. Вычисляем
/"(я?) « Л/0Д7 =
2у/0Д7(0,17 - 3)
1 + 0,17)3
1,457
~*Яс ^ ^г λ2'тах |/,,(:Е = г/5^7)| =
2~°-/»2.1,457<0,01->Л2< *' ™
24
0,08236 -» /г < 0,286,
100-1,457
что позволяет выбрать в качестве h = 0,25 для квадратуры
центральных прямоугольников и вычислить последовательно
/(< х0 >= 1/8) =
/(< X! >= 3/8) =
/(< Ж2 >= 5/8) =
/(< Хз >= 7/8) =
1 _8_
8(1 + 1/64) ~ 65'
3 _ 24
8(1 + 9/64) ~ 73'
5 _ 40
8(1 + 25/64) ~~ 89'
7 56
8(1 + 49/64) 113'
/(< хА >= 9/8) = 72/145; /(< х5 >= 11/8) = 88/185;
/(< Хв >= 13/8) = 104/233; /(< х7 >= 15/8) = 120/289,
и, следовательно,
с2
[2 xdx 1/8 24 40 56 72 88 104 120 \
Jo 1 + χ2 4 \65 73 89 113 145 185 233 289/
ι 0,25(0,12308 + 0,32877 + 0,44944 + 0,49558 + 0,47568 + 0,44635 + 0,41522) ι
«0,25-3,23067 «0,80767.

138
P. S. Интеграл, предложенный для численного решения, имеет
аналитическое решение
/о
Xdx i ln(l + x2)jg = 0,5 In5 ^ 0,80472,
/о 1 + *2 2
подтвержающее точность численного решения.

139
ЛЕКЦИЯ 21
Численное интегрирование (метод трапеций)
7.3. Метод трапеций
Аналогично в методе трапеций подинтегральная функция f(x)
в пределах интегрирования [а, Ь] узлами сетки {жо? #ъ - · · > #п}
разбивается на η равных частей, (Ь — а)/п = h (рис. 7.3.1).
*-#
Χι Хг X} Хя
Рис. 7.3.1. Метод трапеций
Тогда с учетом обозначений рис. 7.3.1 можно записать:
/ /(s)ite = 2 /ι< yi > +Rt = h^ < Vt > +Rt =
= Λ(< J/i > + < J/2 >+···+< У» >+···+< Уп >) +
+ Д« = Л
/(^Ο) , t, Ч , ,/ χ , , /fat»)
+ Rt,
(7.3.1)

140
где
/Ы + /(si)
/(xi) + /Ы
<yi>= <у2>=
/(a?i_i) +/(а;<) /(»n-i) + f(xn)
■ ■ ■, < Vi >= ^ , · · ·, < Уп >= ^
Rt^b-^h2.max\f"(a^x4b)\.
ПРИМЕР 7.З.1. Вычислить интеграл (пример 7.2.1)
r xdx
Ι-
+ x2
по формуле квадратур трапеций с точностью ε ^ 0,01.
РЕШЕНИЕ. Выпишем квадратуру прямоугольной трапеции
S
где
/ f(x)d
J a
\х = h
^ + /(*1) + /Ы + -+^
+ Rt,
<Уг>=т+?+1 (« = 0,1,2,...,η-1), /(х)=х/(1 + х2),
-Rt < l^j^2 · max |/"(о < χ ^ 6)|.
Приведем необходимые результаты расчетов примера 7.2.1:
1-х2
f'(x) =
Я*) =
(1 + х2)2'
2х3 - 6ж
(1 + х2)3'
/"'(*) = "6^3χψ 6 - 0 ^ (1 + χ2) ^ 0 -ч χ4 - 6х2 + 1
= 0 -» х?,2 = 3 ± >/§ = 3 ± 2,83.

141
Итак, х\$ = ±V5,83 > b = 2 не удовлетворяет условию задачи,
жз,4 = ±л/0,17 —> жз = +\/0,17 — корень уравнения.
f"{x « v/0^7 =
2^0^7(0,17 - 3)
1 + 0,17)3
1,457.
Проводим с учетом заданной погрешности ε ^ 0,01 оценку сверху
h:
b — a.
2-0,
Rt ^ -j^h2 ■ max \f"{x = ^/0^7)\ = —/ι2 · 1,457 < 0,01 -*
^h2<
1-6
0,04118->/i< 0,2029;
100-1,457
что позволяет выбрать в качестве h = 0,2 для квадратуры
трапеций и вычислить последовательно
f(x0 = 0) = 0;
/(Жз = 0,6) =
/(я4 = 0,8)
0,6
1 + 0,36
0,8
= 0,34483
= 0,44118
= 0,48780
1 + 0,64
/(х5 = 1) = 0,5;
/(ж6 = 1,2) = 1,2/2,44 » 0,49180
f(x7 = 1,4) = 1,4/2,96 « 0,47297:
f(x8 = 1,6) = 1,6/3,56 « 0,44944
f(Xg = 1,8) = 1,8/4,24 w 0,42453
/(жю = 2) = 2/5 = 0,4

142
и, следовательно,
L
xdx
- 5- = °> 2(0/2+0,19231 + 0,34483+0,44118+0,48780+0,5+
о 1 + я2
+ 0,49180 + 0,47297 + 0,44944 + 0,42453 + 0,4/2) » 0,80097,
причем с учетом точного значения / = 0,80472 абсолютная
погрешность численного вычисления 0,80472 — 0,80097 = 0,00376 <
< ε = 0,01.
7.4. Погрешность квадратур прямоугольников и
трапеций
Точное вычисление собственного интеграла J f{x)dx
возможно, если известна первообразная подинтегральной функции
F(b) - F{a) = / f{x)dx =< f{x) > (b - α), (7.4.1)
J a
причем
</w>.MzfW (7,.2)
или квадратура тождественна F(b) — F(a).
Практический интерес представляет вычисление < f(x) > без
использования первообразной F(x)\^ и
— при известной f(x) или
— при наличии табличных данных (сеточной функции) {ж^, у^:
к = 0, п} эксперимента.
В случае известной подинтегральной функции f(x £ [а,Ь])
вычисление < f(x) > без использования первообразной по сути
сводится к использованию приближенных квадратурных формул
(прямоугольников, трапеций и др.)
М(/(я)) =< f{x) >= - Σ /(< * >)> (7·4·4)
η7ΞΌ

143
абсолютные погрешности которых, используя статистические
характеристики (в частности дисперсию)
И>№)) = - Σ(/(< хг>~< № >)2, (7.4.5)
можно оценить методом итерации на примере квадратурной
формулы центральных прямоугольников.
ПРИМЕР 7.4.1. Оценить, используя дисперсию (7.4.5),
абсолютную погрешность вычисления интеграла примера 7.2.1.
РЕШЕНИЕ. Напомним вычисления примера 7.2.1:
Г2 xdx _ 1 / 8 24 40 Ъ6_ _72_ j*8_ 104 120 \
/о 1 + х2 ~" 4 \6Ь + 73 + 89 + ИЗ + 145 + 185 + 233 + 289/ ~
и 0,25 · 3,23067 и 0,80767,
что позволяет вычислить < f(x) >(7.4.4) и D(< f(x) >) (7.4.4):
< f(x) >= 0,40383; D(< f{x) >) = 0,125[(0,1230 - 0,40383)2+
+(0,32877-0,40383)2+(0,44944-0,40383)4(0,49558-0,40383)2+
+(0,49655-0,40383)2 + (0,47568-0,40383)2+(0,44635-0,40383)2+
+(0,41522 - 0,40383)2] « 0,01383
и сравнить с абсолютной погрешностью точно вычисленного
интеграла
L
Х(лХ χ , /- 2\ι2
5- = ^1п(1 +"* )1о = 0,51n5 « 0,80472,
1 + хг 2
'0
подтвержающее справедливость статистических оценок.

144
ЛЕКЦИЯ 22
Численное интегрирование (метод парабол)
7.5. Метод парабол (Симпсона)
Квадратурная формула Симпсона по структуре сложнее
квадратурных формул прямоугольников и трапеций, однако при
одинаковом η = (b — a)/h обеспечивают более высокую точность
вычислений, так как парабола у = ах2 + Ьх + с сравнительно с
прямоугольниками и трапецией в принципе лучше приближает по-
динтегральную функцию /(ж). Предварительна полезна
Лемма 7.5.1. Площадь S криволинейной трапеции,
ограниченной параболой у = ах2 + Ьх + с, равна
£ = %1+4у2 + уз)/3. (7.5.1)
Доказательство.
/h ph ph ph
[αχ2 + Ьх + c)dx = / ax2dx + / bxdx + с I dx =
-h J-h J-h J-h
ax3 Ьх2 \\ h /e% , ο „ ν
— + — + cx\\ = ~{2ah2 + 6c).
С другой стороны,
Vl\x=-h +4y2U=0 +Vs\x=h =
= {ax2 +bx + c)\x^h+A {ax2 +bx + c)\x==0+ {ax2 +bx + c)\x=h =
= {ah2 -bh + c + ^c + ah2 +bh + c) = 2ah2 + 6c ->
h h
~* S = 3(yi\x=~h + 4^2Ι^=° + Уз\х=н) = τ{2ah2 + 6c),
что доказывает справедливость утверждения леммы. ■

145
Ло Х\ X·) Х"к Х^
Рис.7.5.1. Метод парабол (Симпсопа)
Согласно лемме 7.5.1 отрезок [а, Ь] интегрирования узлами
сетки {#о,#ъ · · · >χη} разбивают на η равных частей, (Ь — а)/п = h
(рис. 7.5.1), причем η = 0(mod2), что позволяет записать
формулу квадратур парабол (Симпсона):
J a
f{x)dx = -[(t/o + 4t/i +t/2) + (j/2 + 4y3 + J/4) + l·
+ · · · + (tfn-2 + 4yn-l + Vn)] + #S = з [(2/0 + Уп) +
+ *{У1 + Уг + · · · + Vn-i)+
+ 2(У2 + У4 + --- + Уп-2)} + Яз,
предельная абсолютная погрешность которой имеет вид
(7.5.2)
Rs<
(Ь-а)*
/ (а < χ < 6)
/г4 (6 — а) »»
180n4 J v~ " ~ " "J 180
Вывод Rs аналогичен методу прямоугольников.
/ (а^х^Ь). (7.5.3)

146
ПРИМЕР 7.5.1. Вычислить интеграл примера 7.2.1
L
2 xdx
по формуле квадратур парабол с точностью ε ^ 0,01.
РЕШЕНИЕ. Напомним квадратуру криволинейной трапеции
ή
S
= / f{x)dx =
J a
= з [(2/0+42/1+2/2) + (2/2+4г/34-2/4) + - · · + (2/η-2+42/η-ι+2/η)]+#5 =
= д[(Уо+Уп) + 4(у1 + 2/з + ••• + tfn-i) +2(2/2 + 2/4 + · ·- + 2/n-2)] + #s,
где
f(x) = x/(l + x% Ис^Н4{Ьта)Г(а^х^Ь).
Приведем необходимые результаты расчетов примера 7.2.1
.,, , 1 — X 2χΆ-6χ
f (Ж) " (1 + Х2)2' / И - (1 +0,2)3'
-6ж4 + 36ж2 - 6
/'"(*) =
(1 + я2)4
/ν _ -24s3 + 72ж (-6ж4 + 36ж2-6)(-4)2ж _
; {Х)~ (1 + я2)4 + {1 + χψ
- 24s5 - 240s3 + 120s
(1 + s2)5
у, . _ 120s4 - 720s2 + 120s (24s5 - 240s3 + 120а?)(-5)2аг
} {Х}~ (Ι + s2)5 + (1 + s2)5
_ 120(1 - s2)(l - 14s2 + s4)
(l + s2)<* ~0,_>
-> (1 + х2) φ 0, (1 - х2) = 0 -> sli2 = ±1;
-> ж4 - 14s2 + 1 = 0 -> х\л = 7 ± л/48 = 7 ± 6,93.

147
Итак, #1,2 = ±1 -> Ει = 1 — корень уравнения; я?з,4 = 7 + \/6,93 >
> 6 = 2 не удовлетворяет; #5,6 = +л/0,07 — корень уравнения.
/"(з« ν/θ^07)«19,5.
Проводим с учетом заданной погрешности ε ^ 0,01 оценку сверху
h:
Rt < ^£hA · max |/""(я = λ/(Μ)7)| < 0,01 -* fc2 ->
loll
-*■ /ι4 < 0,046154 -> Λ < 0,4635,
что в принципе позволяет выбрать в качестве
h = 0,4 -> η = Ц^ = ^-^ = 5 φ 0(mod2).
Выбираем /i = 0,25-+n = 8 = 0(mod2). Однако рационально
выбрать /г = 0,2—>n = 10 = 0(mod2), что совпадает с
узлами примера 7.3.1 и позволяет сравнить результаты расчетов по
различным формулам квадратур:
f(x0 = 0) = 0;
/<*ι = 0,2)-ϊ^=0,1Μ31;
f(x2 = 0,4) = Т-М^= 0,34483;
/(*3 = 0'6) = ОТ = М4118;
/(Χ4 = 0'8) = Π^64)=0'48780;
/(aj5 = l) = 0I5;

148
f{x6 = 1,2) = 1,2/2,44 и 0,49180
/(ж7 = 1,4) = 1,4/2,96 и 0,47297
/(я?в = 1,6) = 1,6/3,56 « 0,44944
f(x9 = 1,8) = 1,8/4,24 и 0,42453
/(ало = 2) = 2/5 = 0,4
и, следовательно,
ί
xdx h
1 + х2 3
[(г/о + Уп) + 4(yi + Уз + · · · + j/„-i)+
+ 2(у2 + У4 + · · · +1/„_2)] = -у ((0 + 0,4)+
+ 4(0,19231 + 0,44118 4- 0,5 + 0,47297 + 0,42453)+
+ 2(0,34483 + 0,48780 + 0,49180 + 0,44944)) =
= ~(0,4 + 4 · 2,03099 + 2-1,77387 =
О
= ^12,0717 = 0,80478
подтвердить
— точность численного решения квадратур парабол, S =
= 0,80767;
— большую точность численного решения квадратур
параболы сравнительно с квадратурами центральных прямоугольников
(0,80478, но при h = 0,25) и трапеций (0,80097 при h = 0,2).
Методы прямоугольников, трапеций и парабол (Симпсона)
являются каноническими, если используют при вычислении
интегралов на равномерной сетке какой-либо один метод с
соответствующей квадратурной формулой. Однако графический анализ по-
динтегральной функции позволяет использовать на каждом шаге
сетки тот метод, который из геометрических соображений
обеспечивают наименьшую абсолютную погрешность вычислений.

149
ЛЕКЦИЯ 23
Численное вычисление несобственных интегралов
7.6. Численное вычисление несобственных интегралов
В предыдущих разделах рассматривались численные методы
вычисления собственных интегралов вида f f(x)dx, причем f(x)
— непрерывная на [а, Ь] функция без особенностей.
Несобственными интегралами называют интегралы с
особенностями
— границ отрезка интегрирования, например, ] — ос, Ь], ]а, ос[,
] — ос, оо[ или
— подинтегральной функции, например, f(x) не существует
или разрывна на отрезке интегрирования.
Существуют аналитические методы вычисления
несобственных интегралов, однако приведенные методы вичисления связаны
с различными приемами использования полученных ранее
квадратурных формул, позволяющими обойти особенности
несобственного интеграла.
Особенности подинтегральной функции f(x).
Если подинтегральная функция f(x) на концах [а, Ь]
обращается в |оо|, то следует воспользоваться квадратурой
центрального прямоугольника (погрешность численного вычисления
несобственного интеграла можно устранить последовательным
уменьшением шага интегрирования, Λ, h/2, h/A и т. д).
ПРИМЕР 7.6.1. Вычислить интеграл
Г1 dx
Ail-*4
с точностью ε ^ 0,01.
РЕШЕНИЕ. Подинтегральная функция f(x) = 1/(1 — χ4) на
концах отрезка интегрирования [—1,1] обращается в бесконеч-

150
ность. Поэтому воспользуемся квадратурой центральных
прямоугольников
S = / f(x)dx = h[f{< х0 >) + /(< χι >) + /(< Х2 >) + · · ·
J a
··· + /(< Яп-1 >)]+Дс.
где
^+Ж;+1
<Λί >=-I-yi±i (i = 0,1,2,... ,η-1), /(я)=я/(1 + *2),
Rc < -^/ι2 -maxI/"(α ^ ж < fc)|,
которая позволяет обойти численно указанные трудности.
Вычисляем последовательно производные подинтегральной функции:
fix) = [(1 - χ4)"1]' = (-1)(1 - х4)~2(-4х3) = (4х3)(1- х4)-2;
f"{x) = [4х3(1 - χ4)"2]' = (3 · 4х2)(1 - х4)~2+
„3,, „44-8/ ,ч/ AJS, 12х2 + 20х6
+ 4х3(1-х4)~3(-2)(-4х3)
(1-х4)3 '
.,„. ч 24х+120х5 (12х2 + 20ж6)(-3)(-4х3)
(1-х4)3 (1-х4)4
= 0:(1-ι4)/0^^ ±1;
120х9 + 240х5 + 24 Λ „ 4
(1 - ж4)4
120х9 + 240х5 + 24х = (120х8 + 240х4 + 24)х =
= (5х8 + 10.x4 + 1)х = 0 -* ая = 0;
5х8 + Юх4 + 1 = 0^ х4,3>4,5 = ~5±^05~4·5
-5±\/5 _ -5 ± 2,236
10 10

151
Итак, #2> #з, ж4> #5 — не удовлетворяет условию задачи, χ ι = 0 —
корень уравнения, но f"(x = 0) = 0, что исключает вычисление
шага интегрирования h.
Прерываем решение примера, чтобы сформулировать
процедуру Канторовича, связанную с переходом от "плохой" (с
особенностями) подинтегральной функции с другой:
— аддитивной, f(x) = ф(х) + д(х),
— мультипликативной, f(x) = ф(х) · д(х), причем д(х)
сохраняет особенность интеграла, но интегрируется проще, а ф(х) на
отрезке интегрирования f(x) образует собственный интеграл.
ПРИМЕР 7.6.1 (продолжение решения).
Представляем подинтегральную функцию в виде
1 0,5 0,5
1-ж4 1-х2 1 + ж2'
причем д(х) = 0,5/(1 — х2) сохраняет особенности
подинтегральной функции f(x) = 1/(1 — ж4), а ф(х) = 0,5/(1 + х2) — функция
без особенностей.
Вычисляем производные д(х) = 0,5/(1 — х2):
f{x) = [(1 - χ2)"1]' = (-1)(1 - х*Г\-2х) = (2х)(1 - χ2)-2;
f"(x) = [2x(l-x2)-2}' =
= (2)(1 - χ2)'2 + 2х(1 - х2Г*(-2)(-2х) = ^±^;
f'M 6'2 , (2 + 6х2)(-3)(-2х)
(1-х2)3 ' (1-х2)4
= 0:(1-х4)^0-»х^±1;
24ж + 24х3 п 1Л 4
(1-х2)4
24х3 + 24х = 24(х2 + 24)х = 0 -> χι = 0;
χ2 + 24 = 0 ->· х2,з = sT^iA.

152
Итак Χ2,χά не удовлетворяет условию задачи, х\ = 0 — корень
уравнения, но f"(x — 0) φ 0. Вычисляем
f"(x = 0) = 2 ^ йс < h-^h2 · max \f"(x = 0)| =
2_0-λ2·2<ο,οι-»λ2^ Ь12
24 " ' ^100-2
«0,06->/ι<0,245,
что позволяет выбрать в качестве h — 0,2 для квадратуры
центральных прямоугольников и вычислить последовательно
/(< х0 >= -0,9) = 1_о81 = 5,2632;
/(< X! >= -0,7) = ^—^ = 1,9608;
/(<x2>=_0,5) = I-^g = l,(3);
/(< х3 >= -0,3) = 1_\ω = 1,0989;
/(< Ха >= -0,1) = /(< хь >= 0,1) = 1, (01);
/(< а* >= 0,3) = 1,0989;
/(<Ж7>=0,5) = 1(3);
/(< я?8 >= 0,7) = 1,9608;
f{<x9 >= 0,9) =5,2632
и, следовательно,
\ ί Г^Ч = °»Х(5'2632 + Х'9608 + 1>3333 + *>0989+
2 J_ι 1-х
+ 1,0101)2 «2,1333.
Интеграл с подинтегральной функцией ^>(ж) = 0,5/(1 + ж2) без
особенностей
1 Г1 dx
2j_1l + x*'

153
предложенный для численного решения, выражается в
элементарных функциях (более того, является табличным интегралом)
/
dx
= arctan χ + const
1 + x
и не представляет интереса в разделе несобственных интегралов
(хотя интересен для численного решения).
Если точка разрыва подинтегральной функции лежит внутри
отрезка интегрирования, то целесообразно
— выбрать шаг h в формулах квадратур, чтобы узлы не
попадали в точки разрыва подинтегральной функции или
— разбить отрезок интегрирования (по числу точек разрыва),
чтобы точки разрыва попали на границу подотрезков
интегрирования, а далее воспользоваться квадратурами центрального
прямоугольника.
Особенности пределов интегрирования
Подинтегральная функция f(x) непрерывна, однако
существуют особенности границ отрезка интегрирования, например,
] — оо, Ь], ]а, оо[, ] — оо, оо[.
Используют прием, связанный с выбором новой переменной
интегрирования, которая устраняет особенность предела
интегрирования. В случае несобственного интеграла f f(x)dx переходят к
новой переменной, которая позволяет перейти к новым пределам
интегрирования без особенностей, например
α α αατ Г х — а -* t — Ο,Ί
α=ζ :->t = l >dx = ~- -^ -» < >
(1 — t)2, у x = оо —> t — 1J
a adt
г и — > dX = — "Τ" Г7Г
1-t χ (I-*)2
adt
и вычислить интеграл методом центральных прямоугольников.

154
ЛЕКЦИЯ 24
Задача Коти ОДУ. Численные методы Пикара и
Эйлера
ГЛАВА 8.Численное решение обыкновенных
дифференциальных уравнений (ОДУ)
8.1. Определения. Задача Коши
Пусть у = у{х) — непрерывная функция одного аргумента,
определенная в открытой (возможно замкнутой) области V(x,y).
Определение 8.1.1. Уравнение F(x,y,yV\... ,y(n)) = 0,
связывающее функцию г/(ж), ее производные (дифференциалы) и
единственную независимую переменную ж, называют
обыкновенным дифференциальным уравнением (ОДУ), причем наивысший
порядок производной (дифференциала) у(х) называют порядком
ОДУ. ОДУ называют линейным, если F(a?,y,j/^\ ... ,у^) = О
линейно относительно у(х),у(т),..., у^.
Замечание 8.1.1. Наряду с F{x^y^y^\ ...,у^) = О
используют обозначения L[y] = /(ж, у), где L оператор, действующий на
у(*)·
Определение 8.1.2. Всякую дифференцируемую η раз
функцию у = у(х), обращающую дифференциальное уравнение η-го
порядка в тождество, называют решением (интегралом,
интегральной кривой) дифференциального уравнения (нахождение решения
дифференциального уравнения называют его интегрированием).
Начальной задачей (Коши) называют нахождение частного
решения дифференциального уравнения η-го порядка,
удовлетворяющее η начальным условиям. Так начальной задачей (Коши)
ОДУ является поиск решения F(x,y(x £ [α,b]),yW) = 0, (у(х =
= α) = consta), которыми и ограничимся далее. Напомним
следующую теорему

155
Теорема (Коши) 8.1.2. Если f(x,y) непрерывна в
замкнутой области Т>(х, у) плоскости Оху и имеет ограниченную
частную производную | $'у' 1 < оо, то для\/{х,у) £ V(x,y)
существует единственное решение у(х) обыкновенного
дифференциального уравнения η = 1-го порядка.
Численные решения ОДУ связаны с дискретизацией f(x,y) €
£ Х>(ж,у), следовательно, у (χ) —» y(xk) = 2/fe, поэтому обязательна
проверка полученного дискретного численного решения
lim yk = у. (8.1.3)
k-+oo
Ограничимся рассмотрением численно-аналитического метода
последовательных приближений (метод Пикара) и собственно
численных методов (метод Эйлера и его модификации) решения ОДУ.
8.2. Метод последовательных приближений (Пикара)
ОДУ
Рассмотрим численно-аналитическое решение
дифференциального уравнения 1-го порядка с начальным условием
—у(х) = у'{х) = /(ж, у), (у{х0) = 2/(0) = (уо) -»
-> L[y) = /(χ,у), Г[у(х = а)] = а. (8.2.1)
Предварительно, напомним, что утверждает
Лемма 8.2.1. Задача Коши для (8.2.1) тождественна
итеративному решению интегрального уравнения
Ук+ι = Уо + f{x}yk)dx. (8.2.2)
Jxo
Доказательство. Для непрерывной на χ Ε [α, 6] функции F(x)
определение первообразной, записанное эквивалентно в
дифференциальной форме
F'(x) = f(x) -* dF(s) = /(z)cte,

156
позволяет у = ук(к == 0,1,2,...) рассматривать в качестве
первообразной F(x) функции f{x,yk), т. е. у'к = F(k = f{x,yk(k =
= 0,1,2,...), что позволяет записать (8.2.1) рекурсивно
рХ ПХ
dy(x) — F(x)dx -» / dy(x) = / F(x)dx -»
J Xq J Xq
pX
I f(x,yk)dx,
Jxa
pX
-* ук+i = 2/o +
что соответсвует утверждению леммы. ■
Суть метода Пикара (последовательных приближений) состоит
в использовании интегрального уравнения (8.2.2) для построения
итеративного алгоритма решения задачи Коши (8.2.1): заданное
начальное условие у(хо) = Уо позволяет переписать (8.2.2)
в первом приближении
/ /(^ Уо)аа
Jxn
2/ι = 2/ο +
и в принципе вычислить j/χ, далее вычислить
второе приближение
Г
2/2 = 2/0 +
'хо
/ f{x,yi)dx,
J X(\
третье приближение
2/3 = 2/0+ / f{x,y2)dx
Jxo
и т. д., что позволяет по индукции сформировать итерационный
алгоритм
Ук =2/0+ / f{x,yk-i)dx,
J xq

157
причем Мшк-^ооУк = У (8.1.3), и решить задачу Коши для (8.2.1).
ПРИМЕР 8.2.1. Решить численно-аналитически задачу Коши
dy/dx = у-х, {у{х = 0) = 2/о = 1)·
РЕШЕНИЕ. Условия теоремы Коши выполнены:
— /(ж, у) = у — χ непрерывна, причем /(ж, у) = 0, если у = ж,
— частная производная ограничена, df(x,y)/dy = 1.
Решаем по схеме Пикара с учетом /(у, х) = у — х:
первое приближение
[х Г1 х2
2/1 = 1+ / /(ж, 2/0)cte = 1 + / (1 - z)cte = 1 + χ - —,
«/ж0 ^0 ^
второе приближение
з
fx ί1 χ2
= 1+/ f(x,y!)dx = l+ {l+x- — -x)dx = l + x
Jxo J0 *
Χ
2/2
третье приближение
fx f1 χ" x^
2/3 = 2/0+/ f(x,V2)dx = l+ {l + x-7r-^-x)dx = l + x~—,
Jxo Jo
2—-x)dx = l + x~ —
что позволяет сформировать по индукции yk+i = 1 + х — η^ и
вычислить
I X \ Χ
у(х) = lim Ук = Ит 1 + ж - —- ) = 1 + х ~ lim —- = 1 + χ.
fc-»oo fc-»oo у k\ J fc-»oo к\
с учетом, что
h
rn"> rn ηη rn rn
lim — = lim — lim - lim — ... lim — = 0.
k—юо h\ fc—юо 1 k—юо 2 fc—юо 3 k—>oo ft

158
Проверка: у(х) = 1 + ж—>t/ = l + # — χ = у — χ.
8.3. Численное решение ОДУ (методы Эйлера)
Численные методы не позволяют получить общее решение
дифференциального уравнения, а приводят лишь к частному
численному решению. Универсальным численным методом решения
дифференциальных уравнений является метод конечных
разностей, сущностью которого является:
— замена области непрерывного изменения аргумента
(аргументов) дискретным множеством точек (говорят сеткой или
сеточной областью) с постоянным (переменным) шагом Л,
например, χ £ [а, Ь] -» Х{ (Ξ [а, 6], где Xi+\ = Х{ + Ах = х\ + Л,
— аппроксимация дифференциального уравнения L[y] = /(ж, у)
дискретным (разностным) аналогом Lh[yh] = /h, которое
является комбинацией значений сеточной функции у(Ш) на шаблоне
сетки (т. е. подмножестве узлов сетки относительно некоторого
ее фиксированного узла) Ш, причем lim/^o \у — УН\ = 0.
Из численных методов решения ОДУ у'(х) — f(x,y), {у{%о) =
= Уо) (8.2.1) рассмотрен метод Эйлера (и его модификации),
основанный на разностном представлении производной, которая для
правой разности, постоянном шаге h = Xk+ι — %k принимает вид
л/ ч // ч ι. Ук+1 ~ Ук Ук+1 ~Ук /о о 1\
/(ж*,Ук) = У (х) = Л hm — ~ τ · 8.3.1)
Окончательно в методе Эйлера для (8.2.1)
Ук+ι = Ук + hy'{x) =ук + hf{xk,yk), (у0 = у{хо)), (8.3.2)
а абсолютная погрешность вычисления следует из разложения
h h2
у{х + К) = у(х) + jy'{x) + — max \y"(a < я <Ь)|, (8.3.3)

159
что позволяет в принципе вычислить шаг Л, если h не задан
условием задачи.
ПРИМЕР 8.3.1. Решить методом Эйлера дифференциальное
уравнение
dy/dx = у - х, (у(х = 0) = г/о = 1), (0 < χ < 0.5).
Вычисления провести при h = 0,1.
РЕШЕНИЕ. Проверка условий теоремы Коши в примере 8.2.1.
Решение методом Эйлера (8.3.2) с учетом f(y,x) = 2/ — #, (у(х —
= 0) =2/о = 1 > 5), (0 ^ χ ^ 0,5) последовательно принимает вид
хо = 0 -» г/1 = г/о + hf(x0,y0) = 1,5 + 0,1(1,5 - 0) =
-1,5 + 0,15 = 1,65;
xi = 0,1 :2/2 = 2/1 + hf{xuyi) = 1,65 + 0,1(1,65 - 0,1) =
= 1,65 + 0,155 = 1,805;
*2 = 0,2 :2/з = 2/2 + hf{x2,2/2) = 1,805 + 0,1(1,805 - 0,2) =
= 1,805 + 0,1605 = 1,9655;
хз = 0,3 :2/4 = 2/з + Л/(я?з, Уз) = 1,9655 + 0,1(1,9655 - 0,3) =
= 2,13205;
ж4 = 0,4 -> 2/5 = 2/4 + Л/(я?4,У4) = 2,13205 + 0,1(2,13205 - 0,4) =
= 2,305255;
ж5 = 0,5 -* г/б = 2/5 + Л/(ж5,2/δ) =
= 2,305255 + 0,1(2,305255 - 0,5) = 2,48578.
Проводим вычисления абсолютной погрешности приведенных
расчетов: \таху"(х)\ = \у'(х) — χ — 1\ = 0,805255 реализуется
при χ = 0,5; Тогда
^-|maxy'(s)| = ^--0,805255 « 0,004026275. ■

160
ПРИМЕР 8.3.2. Решить методом Эйлера дифференциальное
уравнение
dy/dx -у-χ, {у{х = 0) = 2/о = 1,5), (0 < χ < 0,5).
Вычисления провести с абсолютной погрешностью ε = 0,01.
РЕШЕНИЕ. Проверка условий теоремы Коши в примере 8.2.1.
В качестве предварительной оценки h вычисляют у" (х = 0):
-г-у'{х) = -7-0/ ~ х) = 2/'0*0 -1 = 2/-ж-1-> у"(х = 0) =
= 1,5-0-1 = 0,5,
откуда с учетом (8.3.3) следует
у"(0) 0,5
Проводим предварительный расчет при /ι ^ 0,2. Например,
можно использовать данные расчетов предыдущей задачи, где h =
= 0,1 ^ 0,2. Тогда последовательно получаем
хо = 0,2/о = 1,5 -)· у"(жо) = 1,5-0- 1 = 0,5;
Ж1 = 0, l;yi = 1,65 -» у"(жх) = 1,65 -0,1 - 1 = 0,55;
Ж2 = 0,2; у2 = 1,805 -► у"(ж2) = 1,805-0,2- 1 = 0,605;
хг = 0,3; уз = 1,9655 -> у"(ж3) = 1,9655 - 0,3 - 1 = 0,6655;
х4 = 0,4; j/4 = 2,13205 -» у"(ж4) = 2,13205 - 0,4 - 1 = 0,73205;
Х4 = 0,5; № = 2,305255 -» у"(ж5) = 2,305255 - 0,5 - 1 =
= 0,8053255.
Пересчитываем h:
"2«^ = δβί0·02484"ΛΐΜ57β·
Итак, расчет можно проводить при h = 0,1.
Следовательно при решении задачи Коши примера 8.3.2 можно
воспользоваться данными расчетов примера 8.3.1. ■

161
ЛЕКЦИЯ 25
Численное решение краевой задачи ОДУ
8.4. Определения. Краевая задача ОДУ
Определение 8.4.1 Обыкновенное дифференциальное
уравнение
F(x>y(x),y{I)Jn\---,y(n))=0 (8-4.1)
с η ^ 2 краевыми условиями называют ОДУ с краевыми
условиями, причем построение частного решения ОДУ,
удовлетворяющего краевым условиям, называют краевой задачей ОДУ.
Численное решение ОДУ 2-го порядка рассмотрим на примере
F(x,y(xe[a,b}),y{n,y{n)) = 0,
(у(х = а) = consta,y(x = b) = constb) —» (8.4.2)
L[y] = /(*, у), Г[у(х = α)] = α, Г[у(х = 6)] = /J,
которое является краевой задачей ОДУ, причем в качестве
граничных значений могут выступать ±оо.
8.5. Численное решение краевой задачи ОДУ
При решении краевой задачи линейного ОДУ второго порядка
(8.4.2) с учетом, что справедливо
Допущение 8.5.1. О Зу(х) краевой задаче линейного ОДУ
(8.4.2). Ограничимся вычислительными методами наименьших
квадратов (интегральный вариант) и разностной факторизацией,
рассмотренные в настоящих лекциях.
Метод наименьших квадратов
Решение краевой задачи линейного ОДУ
L[y] = у" + р(х)у' + q(x)y = f(x) (8.5.1)

162
допускает представление
η
(8.5.2)
г=1
где фг — базисные (линейно независимые) функции и
а, г = О,
Г[фг(х = а)]
О, i > О,
Г[Л(х = Ь)]| Д
краевые условия. С учетом представления невязки (остатков)
/3, < = О,
г> О
Я(*, d, С2) = Цфо] " /(*) + ^ СМх)
г-1
условие минимальности интегрального среднеквадратичного
приближения (невязки)
ι = Ι*(ΐ[φο]-ί(χ)+Σαφ{(χ)
-> mm,
(8.5.3)
позволяет определить постоянные С\,Сч из решения СЛАУ (с
симметричной матрицей коэффициентов):
( д!{СъС2)
dI(CltC2)
дСо
= о,
= 0
(8.5.4)
и записать частное решение краевой задачи ОДУ.

163
Замечание 8.5.1. Если Г[у(х = а)] — а — О, Г[у(х = Ь)] =
= β — О, т. е. краевые условия однородны, то можно положить
0о(*)=О.
Замечание 8.5.2. Можно использовать точечный вариант
метода наименьших квадратов, если заданную область
интегрирования краевой задачи ОДУ заменить сеточной функцией (в первом
приближении равномерной и плотной).
Метод разностной факторизации (метод прогонки)
Чтобы упростить пояснения, метод разностной факторизации
рассмотрим на примере решения линейного ОДУ 2-го порядка
Щ = у" + q(x)y = /(*), (Г[у(х = а)} = ув> Г[у(х = Ь)} = до),
(8.5.5)
которое перепишем в конечных разностях на трехточечном Ш=
= {xi+i,Xi,Xi-i}:
Lh[y } = J~2 + ЯгУг = fi ->
-> y<+i - 2у< + Уг-1 + ЯгУгЬ2 = /*Л2. (8.5.6)
Последнее выражение (8.5.6) с учетом
й_1 + [^ - % + у<+1 = h2fu {г = 1,2,... ,п - 1) (8.5.7)
допускает представление в виде СЛАУ Aj/ = Ь, которое можно
решить методами линейной алгебры, например методом
прогонки с учетом Г[у(х = а)] — ?/а, Г[у(х = Ь)] = до) или используя
матричную факторизацию (1/[/-разложение).
ПРИМЕР 8.5.1. Решить методом наименьших квадратов
(интегральный вариант) обыкновенное дифференциальное уравнение
2-го порядка с краевыми условиями
у" + (1 + х2)у + 1 = 0, {у{х = -1) = у{х = 1) = 0).
РЕШЕНИЕ. В качестве базисных функций выбираем
полиномы
фп(х) = х2п~2{1 - χ2), η = 1,2,....

164
С учетом порядка решаемого дифференциального уравнения
достаточно ограничиться:
η = 1 —» φι(χ) = 1 — #2,
п = 2-*ф2{х) = х2{1-х2),
что позволяет записать общее решение дифференциального
уравнения в виде
у{х) = С1{1-х2) + С2х2{1-х2),
которое удовлетворяет граничным условиям краевой задачи ОДУ.
С учетом невязки (остатков)
η
R(x, CUC2) = Цф0] - f(x) + Σ Цф{(х) =
t=l
= 1 - (1 + x4)d + (2 - Их2 - х6)С2
записываем интегральный вариант метода наименьших
квадратов (8.5.3):
/ = / [1 - (1 + x4)Ci + (2 - Их2 - x6)C2]2dx -> min,
минимизируя которое, получаем СЛАУ относительно
неизвестных С\,С2:
д!
dd
-2 ί [1 - (1 + x4)Ci + (2 - llx2 - x6)C2](l = x4)dx = 0,
£r=2 [1 - (1 + x4)Ci + (2 - Их2 - x6)C2](2 - llx2-
C02 J-i
- x6)dx = 0,

165
которое после интегрирования можно переписать в виде
( 68 ,3548 5
I 45Сг + Ш5С2 = 1>
] 3548^ 63404 38
I. 1155 1 4095 21
Решая последнюю СЛАУ (с симметричной матрицей
коэффициентов), получаем С\ = 0,985, Сч = —0,078, что позволяет записать
частное решение
у(х) = 0,985(1 - ж2) - 0,078ж2(1 - ж2)
краевой задачи ОДУ. ■

166
ЛЕКЦИЯ 26
Численное решение уравнений в частных производных
ГЛАВА 9.Численное решение уравнений в частных
производных
9.1. Определения. Классификация
Если у = у(х) — непрерывная функция более чем одного
аргумента, определенная в замкнутой области Х>, то ее дифференциал
имеет вид
dy=mdXl+?Mdx^...+mdXitt ρ,,.,,
ΟΧ χ ΟΧ2 ОХ п
соответственно, и дифференциальное уравнение (называемое
дифференциальным уравнением в частных производных) в принципе
учитывает частные производные у = у(х) по тем и (или) иным
переменным и поэтому приобретает более сложный вид
сравнительно с ОДУ. Теория дифференциальных уравнений в частных
производных достаточно сложна и не будет нашим предметом. Более
частными (но полезными для аналитиков) являются уравнения
процессов и систем математической физики, с которыми
связывают параболические, элиптические и гиперболические
дифференциальные уравнения в частных производных:
μ α— = f(t^x) — уравнение переноса,
at ox
dU{t,x,y) u dU2{t,x,y) . , dU2{t,x,y)
ot = kx ox* + ку~~д? УР^нение
теплопроводности,
dU2(x,y) dU2{x,y) ,, ч „ ,Λι1Λ4
^V^ + hV"^ = Яж> У) —Уравнение Пуассона, (9.1.2)
охг иуг

167
параболическая, эллиптическая, гиперболическая "природа"
которых проистекает из структуры их характеристических
уравнений. Предметом последующего рассмотрения будут простейшие
дифференциальные уравнения в частных производных
параболического типа (с которыми традиционно связывают задачи
стационарной и нестационарной теплопроводности) и их методы
численного решения.
9.2. Численное решение уравнений в частных
производных
Численное решение дифференциальных уравнений в частных
производных рассмотрено на примере двумерного уравнения
нестационарной теплопроводности
dT{t,x,y) _ и dT*(t,x,y) , L dT\t,x,y) fc*= const*,
dt ~x дх2 +у ду1 ' fc^constj,, ν·ΖΛ)
которое решаем традиционной численной процедурой —
разностной факторизацией, предварительно упростив структуру (9.2.1).
Стационарный случай, двумерная среда
При
«•('■«•»)-0 (9.2.1)
at
принимает вид дифференциального уравнения в частных
производных 2-го порядка для стационарной теплопроводности
двумерной термоизотропной среды
к972^ + кдТ^У) =0,(кх = ку = к = const φ 0), (9.2.2)
которое с учетом определений частных производных в разностной
форме (правая разность)
дТ(х,у) _АТ = Ti>j+l-Titj^ дТ(х,у) „AT = Ti+1j-Tij
дх Ах hx ' ду Ay hy

168
дТ*(х,у) _ Δ (ϋ) = A(Titj+1 - Titj) = TiJ+l - 2Titj + Tid.x
дх2 ~ Ax hi hi
дТЦх,у) ^АШ
A(Ti+ij —Tjj) __ Tj+ij - 2Tjj + Tj-ij
ду2 ~ Ay hi Щ
(9.2.3)
принимает вид дискретного аналога дифференциального
уравнения в частных производных 2-го порядка с граничными
(краевыми) условиями
k _ +k _ _o,
T(0,y) = Το,Γ(χ,Ο) = TuT(xn,y) = T2,T(x,yn) = T3, (9.2.4)
которое можно представить, например, на двумерной
равномерной (hx = h) сетке (0 ^ i, j ^ п) в виде в СЛАУ и решить
относительно VTij методами линейной алгебры.
Стационарный случай, одномерная среда: Τ = Τ (χ)
dT(x) = Щх± йДТ = Tj+i-Tj
dT2 = ЭГ2(Ж) ^ A (H) = A(TJ+i - Tj) = Ti+1 - 27) + T,·-!
cte2 дх2 ~ Ах hi hi
(9.2.5)
принимает вид дискретного аналога ОДУ 2-го порядка с
граничными (краевыми) условиями
fcri+1-2T-,+Ti_1 = 0 {τ{χ = 0) = Γθ)Τ(χη) = Τι) (9 2 6)

169
Нестационарный случай, одномерная термоизотропная
среда
(9.2.1) при d2T(t,x)/dy2 = 0 принимает вид
дифференциального уравнения в частных производных 2-го порядка для
нестационарной теплопроводности одномерной термоизотропной среды:
д~^т^ = кх^^ {к* = constф0)) (9-2J)
которое с учетом определений частных производных в разностной
форме (правая разность)
dT(t,x) _ ΔΤ _ TT|i+i -TrJ
дх Ах hx
dT2{t,x) _ Δ|^ = A(TT|j+i -TTj) = Trj+i - TTrj + Trj-i
dx2 ~ Ax h\ h2x
(9.2.8)
9T(t,x) AT _ TT+hx - TT,X . .
~~dt Δί" ~ Ы (9 9)
принимает вид дискретного аналога дифференциального
уравнения в частных производных 2-го порядка с граничными
(краевыми) условиями
Ττ+ι,χ — ΤΤιΧ _ TTj+i — 2Trj + TT)J-i
ы hi
T(t = 0, ж = 0) = T0,T(t = 0, ж = 0 = Τι. (9.2.10)
Уравнение (9.2.10) с учетом граничных условий можно
представить, например, на одномерной равномерной (hx = К) сетке (0 ^
^ 3 ^ п) в виде в СЛАУ и решить относительно VTT); методами
линейной алгебры. Для решения дифференциальных уравнений
в частных производных используем термодинамическое развитие

170
метода конечных элементов (МКЭ) — метод конечных
термодинамически равновесных элементов (МКТРЭ).
9.3. Численное решение задач теплопроводности.
МКТРЭ
Метод конечных элементов (МКЭ) широко распространен
среди вычислительных методов в силу его универсальности, легкости
физической интепретации и относительной простоты.
Основная идея МКЭ состоит в том, что решение краевой
задачи ищется в виде набора функций (обычно это полиномы),
определенных на некоторых сеточных областях (сеточных конечных
элементах), что позволяет перейти от дифференциальных
уравнений к системе линейных алгебраических уравнений (СЛАУ).
Для плоской задачи наиболее часто применяют КЭ в виде
треугольников и четырехугольников, причем наиболее простые
(линейные) КЭ имеют узлы только в вершинах и прямые стороны. В
одномерном случае линейный элемент — отрезок прямой с узлами
на концах. В качестве функции КЭ обычно применяют полином,
порядок которого зависит от числа узлов элемента. В двумерном
случае линейный элемент имеет форму треугольника. В
трехмерном случае линейный элемент имеет форму тетраэдра. Иначе
говоря, в линейном случае число постоянных полинома равно числу
узлов КЭ, а узлы совпадают с границами или вершинами
элемента. Для сравнения оптимальности вариантов разбиения тела
используют значения математического ожидания Μ и дисперсии
D: вариант с меньшими МиВ соответствует более оптимальному
разбиению сеточной области и потому предпочтителен.
Термодинамическая "некорректность" дискретных уравнений
теплопроводности (где Tij неизвестная температура в ij-м узле
(!) сетки) с учетом улучшения МКЭ его термодинамическим
аналогом — методом конечных термодинамически равновесных
элементов (МКТРЭ) позволяет переписать (9.2.4)(9.2.6) в виде

171
fc— -2 - 0,
(T(0,y) = T0,T(a?,0) = Γι,Τ(χη,ι/) = T2>T{x,yn) = T3), (9.3.1)
fcTj+1-2Tj + TJ-1 = 0 (т(ж = o) = ^ r(jBn) = Γι))
где Tij — температура не в ij~ом узле сетки, а в ij-й замкнутой
(безузельной!) области: для одномерной задачи — отрезок
стержня, для двумерной задачи — квадрат(треугольник), для
трехмерной задачи — куб (тетраэдр). Идеологию МКТРЭ уточняют
замечания
Замечание 9.3.1. Неизвестные Г^(или Tj) КТРЭ п-мерных
задач (9.3.1) — это статистические (арифметические) средние
соприкасающихся (по границе, а не в точке) соседних конечных
элементов (так называемых "ближайших соседей"):
(г) η = 1(стержень) : Г,- = (Tj+i + Τ,-0/2 (j = 1,2,...),
f Δ: Tij = (Γ^_ι + Ti+ltj + Γ<ιί+ι)/3, (г J = 1,2,...)
(η) η = 2 <
Ι Π : Τ*,, = №,,·_! + Ti-ij + Τ*j+i + Т<+и)/4 и др.,
(9.3.2)
что термодинамически обосновано (в термодинамике говорят не о
температуре в точке, а температуре области). Проверка
устойчивости полученных результатов осуществляется статистическими
средними (Μ(·),Β(·)) на статистике КТРЭ при их
последовательном измельчении (вплоть до получения устойчивых результатов).
Замечание 9.3.2. МКТРЭ не ограничен уравнениями
теплопроводности и распространяется на любые дифференциальные
уравнения, для неизвестных которых физически естественна
макроскопическая природа, например давление, энергия, энтропия и
др.
Программа (Turbo Pascal 6.0) решения методом Гаусса
уравнения (9.3.2, (i)) приведена в приложении, в котором для
наглядности вложен отладочный пример для двумерной термоизотропной

172
двумерной области (пластины) с граничными условиям,
разбитой КТРЭ (рис. 9.3.1). Такой подход можно использовать и для
решения нестационарных задач, рассматривая область конечных
термодинамически равновесных элементов задачи на каждом
шаге t(τ) как страницу книги, каждая страница которой
пронумерована временем t(r + 1) = t(r) + At, τ = 1,2,3,
*(/)
0
To
y(0 \
τ
I -mi
Τ
121
Τ
131
τ
τ,
τ12
τ
Χ22
Τ
132
Τ3
τ13
τ„
Тзз
Χ
Τ
Рис. 9.З.1. Область безузельных квадратных
элементов уравнения двумерной стационарной изотропной
теплопроводности (МКТРЭ)] (9.3.2)
Матрица коэффициентов (А) СЛАУ разностных уравнений при
большом количестве конечных элементов становится
разреженной, что приводит к проблеме рационального хранения
разреженных матриц в памяти ЭВМ. Это не единственная проблема
разностных уравнений: матрица коэффициентов СЛАУ разностных
уравнений не только разрежена, но рыхло заполнена
коэффициентами, с чем связана проблема рациональной перенумерации узлов
сетки для формирования компактной ленточной матрицы [4].

173
ПРИЛОЖЕНИЕ
МЕТОДЫ ПОСТРОЕНИЯ
АНАЛИТИЧЕСКИХ МОДЕЛЕЙ
В приложении изложены некоторые методы вычислений,
полезные при построении аналитических (физических и физических
аппроксимационных) моделей прикладной математики.
П.1. ПРИМЕР ФИЗИЧЕСКОЙ МОДЕЛИ
Модель плотного газа (аксиоматический подход) [16]
Необходимость усовершенствования модели газа Ван-дер-
Ваальса связана с ее несовершенством, вызванным
ограниченностью модели твердых сфер с притяжением. Например, структура
свободной энергии газа Ван-дер-Ваальса
i/>(T,V) = -NTIn [T3/2eV(l - Nb/V)/(NA3)j - N2a/V,
где a — постоянная притяжения, b = 2π£>3/3, — не обеспечивает
правильный ход термодинамических величин: изохоры газа Ван-
дер-Ваальса — прямые линии, а теплоемкость Cv совпадает с Cv
идеального газа. Однако более существенным недостатком газа

174
Ван-дер-Ваальса, как впрочем и других сред, описываемых двух-
параметрическими (например вириальными) уравнениями
состояния в том, что они удовлетворяют закону соответственных
состояний, следовательно, описывают термодинамически подобные
газы. Так, с учетом соотношений {дР/дУ)ткр = (дР2/д2У)ткр = О
в критической точке (индекс "кр") приведенное уравнение
состояния газа Ван-дер-Ваальса принимает вид:
Р/Ркр = 8(T/TKp)/[S(V/VKp) - 1] - 3(VKP/V)2,
где а — 3P7cp(Vr?Cp/iV)2, Ъ = VKp/(3N), причем для реальных
веществ Ζκρ = ΡκρνκρΙ(ΝΤκρ) = 3/8 не согласуется с опытом.
Несовпадение с опытом вызвано тем, что "...широкий класс газов
ведет себя почти в соответствии с потенциалом Lennard-Jones "12"
"6", который рассматривает атомы как
сферически-симметричные поля, а не как твердые сферы". Запишем кластерную модель
частиц плотного газа.
АКСИОМА (статики, сплошности модели) 1. Атомы газа
распределены в объеме V равномерно с плотностью 1/v:
V = Nva + Nvvv = Νυ, ναφνν, (Π.1.1)
где να = 7τ£>3/6, D — диаметр атома, а смысл vv выясняет
АКСИОМА (динамики модели) 2. Атомы модели
перемещаются поступательно в объеме V и сталкиваются друг с
другом, причем линейный размер vv соответствует средней длине
свободного пробега атомов модели.
Тогда с учетом замечаний, что взаимодействие между
атомами газа не очень мало только при столкновениях, т. е. когда
расстояние между ними соответствует диаметру D атома, верна
АКСИОМА (взаимодействия ближайших соседей) 3.
Потенциальная энергия Ής({χ}) плотногазовой системы с неупругими
столкновениями определена в виде суммы бинарных — 7iq(xi,Xj),
триарных — Wq(xi,Xj,xi),..., п-арных — %q(xiiXj,...,Xk)
взаимодействий (столкновений) атомов модели с его ближайшими
соседями:

175
Kq{{x}) =^q(xuXj) +Hq{Xi,Xj,Xl) Η + Uq(xi, Xj, . . . , £fc),
(Π.1.2)
причем, например,
1 W
и,(ЗД) = 2 Σ HRiiWRii ~ D),
F(R-) = ί €°lm(R™/Rij)n -n{RmIRij)m}l{n-m), при Ду = D
[ ij) lO при Ri:i > D.
(П.1.3)
АКСИОМА (столкновений) 4. Лгаолш п^ш столкновении
образуют кластер, включающий центральный атом и его
ближайших соседей.
Итак, столкновения атомов модели плотного газа не
обязательно парные, что отразится на величине решетчатых сумм Ак:
Ак = J2i X^i==1 ni(Ri/Ri)K, где I = 1,2,3,... — число кластеров
модели в каждый момент времени t.
С помощью аксиом 2,3 можно записать Ής({χ}) в виде
ад*>) = Σ ^|^ед, - я) =
гфЗ
Ne0
2(п - т)
Апт
Атп
(П.1.4)
когда Tiq({x}) не зависит от переменных интегрирования:
NlJv'Jv
ехр
Апт
Nep/2
η — т
γΝ
Ж6ХР" 2(n-m)
- Атп
dVi...dVN =
Ne0
Апт
- Атп

176
что позволяет записать выражение для свободной энергии ч/>(Т, и,
να), следовательно, выразить термодинамические функции
плотного газа в неравновесном состоянии:
Ε(Τ,υ,υα)
Ne0/2
η — πι
Anm
- АтП
m/3
+ 3NT/2
Ρ(Τ,υ,να) = -(дф/dV) = NT/V,
S{T,v,va) = iVln[T3/2Ve5/2/(^A3)].
(П.1.5)
В равновесном состоянии в принципе можно пользоваться ψ(Τ,ν,
να), однако теперь не все Τ,ν,να являются независимыми
переменными, а связаны условием минимума свободной энергии дф(Т,у,
va)/dva = 0:
θψ(Τ,υ,υα)
0va
NT Ne0 nm
ν 6να η ~ τη
Vm
n/3
-A
m/3"
= o,
l/v =
eo nm
6Tva η — πι
An[^
n/3
«/»'
(Π.1.6)
Условие минимума (П.1.6) позволяет рассматривать ψ(Τ,ν,να) в
качестве свободной энергии Гельмгольца модели плотного газа,
следовательно, выразить его уравнения состояния:
Е =
Nen
2(п — т)
Апт
Λγηϊϊ
+
\««/
SNT
2
+
3υ
(dva\ 1
Кдт)у\

177
dVJT ν [ \θυ
T3/2ye5/2 NT /dVa
ТЛ
iVA3 ν
где (Ουα/θΤ)ν = -(v/T)(dva/dv)T и
2
дТ
(П.1.7)
(dv/dva)T = -
nmco
vaJ 6T(n — m)
Полученные функции состояния с учетом постоянных п,га,бо,^т
могут быть использованы для описания поведения плотных газов
в окрестности критической точки:
Ζκρ = 1 - (θνα/Θν)Τκρ = -v(d2va/dv2)Tl[p, (d3vJdv3)TKp = О
на границе испарения и при описании процессов различной
природы.
П.2. РАЗНОСТНЫЕ УРАВНЕНИЯ
Полезный результат следует для одномерного стационарного
процесса волнового сжатия вещества (т. е. г = j ^ \/σ^· φ
Φ 0,Vcy = 0 и ε33 φ 0; i φ j ■$ Va<j,6y = 0, где ση = σ22 =
= —Ρ + σ/2, азз = — {Ρ + о) для твердого тела и оц = —£>ijP для
жидкости), причем для идеальной жидкости законы сохранения
принимают вид
Vvs = Vfcfae -г;р),
Ρ - Ρ0 = «aVp/Vb,
£-E0-(P + P0)(Vo-^)/2,

178
где vp — массовая скорость вещества за фронтом vs ударной
волны (скорость деформации вещества), а индексом "О" помечено
начальное состояние вещества, Ρ = Pq + ΔΡ, V = Vo + AV и др.
При
At =» 0 : d£(P, V) = (V& - K)dP/2 ~ {Ρ + P0)dV/2 = -iW,
что характерно для изэнтропических процессов изолированной
системы, а за интервал времени At φ 0 изменение энергии
тела вдоль адиабаты Гюгонио имеет вид
АЕ = -Ρ(ΔΚ) + (ΔΡ)(Δ7)/2,
что характерно для необратимых процессов адиабатически
изолированной системы, причем
-APAV/2 = vJ/2.
П.З. ЗАМЕНА ПЕРЕМЕННЫХ ФУНКЦИЙ
В вычислительной математике представляет интерес задача,
связанная с изменением множества (подмножества) переменных
известной непрерывной функции, например, f(x,y) -> f(u,v), где
u(x,y),v(x,y).
Такая задача может привести к более наглядному
вычислительному эксперименту и более простым дальнейшим
преобразованиям над исследуемой функцией с иными переменными.
Например, в классической механике совершается переход в
гамильтониане
?ί({»ι,α:2,··.,»η},{έι,±2,···,*η}) ->
->^({ж1,Ж2,...,Жп},{Р1)Р2,...,Рп}),

179
что позволяет, рассматривая новые переменные как независимые,
разделить их. Сказанное проиллюстрируем примером
рационального преобразования функции z{x,у; α,6, с), где х,у —
независимые переменные, а, Ь, с — параметры:
г(ж,!/;а,6,с) = хау / ехр ( )ydy+—ln exp(--)ydy,
Jo v я/ У Jo V x'
(П.3.1)
вводя новую переменную
ξ(χ,ν) = / ежр (--) ydy, (П.3.2)
упрощающую, например преобразования (П3.1),
xh
z{x,y^;a,b) = хаУ ξ+ -1ηξ, (Π.3.3)
хотя в данном случае трудностей с интегрированием (П.3.2) нет,
например
«*; с) = jf1 ехр (-£) ydy = f^Z^M, (П.3.4)
причем иное выражение для ξ получим из условия
9z{x,y,£) ay>xb η ν C( ΙΛ Хб te1"^
Οξ у£ ухаУ у
(П.3.5)
Подстановка (П.3.5) в (П.3.3) позволяет записать окончательно
хЬ хЪ( Ьх1~ау\ /тто.ч
2? х,у;а,6) = + — In , (П.3.6)
У У \ У )
что с учетом (П.3.4) тождественно 2г(ж,у;а,Ь, с). По сути
введение новой переменной ξ связано с исключением параметра с:
ζ(χ9 у; а, δ, с) -» ζ(χ, у, ζ; а, Ь).

180
Π.4. ΠΡΟΓΡΑΜΜΑ РЕШЕНИЯ МКТРЭ
ДВУМЕРНОГО СТАЦИОНАРНОГО
УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ
Замечание. Для наглядности в программу вложен отладочный
пример для двумерной термоизотропной сеточной области
размером 3x3 КТРЭ с граничными условиями (рис. 9.3.1).
program gauss;
const а[1,1]=4;а[1,2]=-1;а[1,3]=0;а[1,4]=-1;а[1,5]=0;
,а[1,6]=0;а[1,7]=0;а[1,8]=0;а[1,9]=0;а[1Д0]=Го + Τι=55;
а[2,1]=-1;а[2,2]=4;а[2,3]=-1;а[2,4]=0;а[2,5]=-1;
а[2,6]=0;а[2,7]=0;а[2,8]=0;а[2,9]=0;а[2,10]=Г1=30;
а[3,1]=0;а[3,2]=-1;а[3,3]=4;а[3,4]=0;а[3,5]=0;
а[3,6]=-1;а[3,7]=0;а[3,8]=0;а[3,9]=0;а[3,10]=Т1 + Т2=70;
а[4,1]=-1;а[4,2]=0;а[4,3]=0;а[4,4]=4;а[4,5]=-1;
а[4,6]=0;а[4,7]=-1;а[4,8]=0;а[4,9]=0;а[4,10]=Т0=25;
а[5,1]=0;а[5,2]=-1;а[5,3]=0;а[5,4]=-1;а[5,5]=4;
а[5,6]=-1;а[5,7]=0;а[5,8]=-1;а[5,9]=0;а[5,10]=0
а[6Д]=0;а[6,2]=0;а[6,3]=-1;а[6,4]=0;а[6,5]=-1;
а[6,6]=4;а[6,7]=0;а[6,8]=0;а[6,9]=-1;а[6,10]=Г2=40;
а[7,1]=0;а[7,2]=0;а[7,3]=0;а[7,4]=-1;а[7,5]=0;
а[7,6]=0;а[7,7]=4;а[7,8]=-1;а[7,9]=0;а[7,10]=Г0 + Г3=70;
а[8Д]=0;а[8,2]=0;а[8,3]=0;а[8,4]=0;а[8,5]=-1;
а[8,6]=0;а[8,7]=-1;а[8,8]=4;а[8,9]=-1;а[8,10]=Т3=45;
а[9,1]=0;а[9,2]=0;а[9,3]=0;а[9,4]=0;а[9,5]=0;
а[9,6]=-1;а[9,7]=0;а[9,8]=-1;а[9,9]=4;а[9,10]=Т2 + Г3=85;
type mat=array [l..20,1..21] of real;
vec=array [1..20] of real;
var a:mat; x:vec; i,N :integer; s:real;
procedure matr(N:integer; var a:mat); var ij-.integer;
begin for i:=l to N do for j:=l to N+l do
begin write('a[',i:2,Vj:2,']='); read(a[ij]) end end;
procedure gauss(N:integer;var a:mat;var x:vec;var s:real);

181
var ij,k,l,kl,Nl:integer; r:real; begin
N1:=N+1; for k:=l to N do begin kl:=k+l; s:=a[k,k]; j:=k;
for i:=kl to N do begin r:=a[i,k]; if abs(r)»abs(s) then begin s:=r;
j:=i end end;
if s=0.0 then exit; if j«»k then for i:=k to N1 do begin r:=a[k,i];
a[k,i]:==a[j,i]; a[j,i]:=r end;
for j:=kl to N1 do a[kj]:=a[kj]/s;
for i:=kl to N do begin r:=a[i,k]; for j:=kl to N1 do a[i j]:=a[i j]-
a[kj]*r end end;
if s«»0.0 then for i:=N downto 1 do begin s:=a[i,Nl]; for j:=i+l to
N do s:=s-a[ij]*x£j]; x[i]:=s end end;
begin repeat
\νΓΗβ('Ν='); readln(N); matr(N,a); gauss(N,a,x,s);
if s«»0.0 then
for i:=l to N do ιντίΐβΙηί'χ',ίΛ,^',χΗ) else writeln('(iet=0,)
until false
end.

182
ЛИТЕРАТУРА
1. Бахвалов Я. С., Жидков Я. Я, Кобельков Г. М.
Численные методы. — М.: Наука, 1987.
2. Волков Е. Л. Численные методы. — М.: Наука, 1987.
3. Демидович Б. Я., Марон Я. А. Основы вычислительной
математики. — М.: Наука, 1970.
4. Джордж, А,^ Лю Дж. Численное решение больших
разреженных систем уравнений. — М.: Мир, 1984.
5. Калиткин Я. Я. Численные методы. — М.: Наука, 1978.
6. Косарев Я. И. 12 лекций по вычислительной математике.
— М.: Изд-во. МФТИ, 1995.
7. Тихонов Л. Я, Костомаров Д. Я Вводные лекции по
прикладной математике. — М.: Наука, 1984.
8. Ильин Ε. Α., ПознякЭ. Г. Линейная алгебра. —М.: Наука,
1984.
9. Ильин Я. Я. Методы неполной факторизации для решения
алгебраических систем. — М.: Наука, 1995.
10. Самарский А. А. Введение в численные методы. — М.:
Наука, 1987.
11. Севастьянов Б. А. Курс теории вероятностей и
математической статистики. — М.: Наука, 1982.
12. Фадеев Д. К,, Фадеева Я. Я. Вычислительные методы
линейной алгебры. — М.: Физматгиз, 1963.
13. Форсайт Дж., Малькольм М., Моулер К. Машинные
методы математических вычислений. —М.: Мир, 1980.
14. Хорн Р., Джонсон Ч. Матричный анализ. — М.: Мир,
1989.

183
15. Протасов И. Д. Уравнения термодинамических состояний
биологических сред в кн.: Современные проблемы биомеханики.
Вып. 6. — Рига: Зинатне, 1989.
16. Протасов Я. Д. Уравнение состояния с параметрами
некоторых сплошных сред// Вестник МГУ — Сер.1 Математика.
Механика. — 1997. — 2. С. 35-39.
17. Протасов И. Д. Теория игр и исследование операций. —
М.: Гелиос АРВ, 2003.
Сборники задач и программ, практикумы
18. Воробьева Г. Н., Данилова А. Я. Практикум по
вычислительной математике. — М.: Высшая школа, 1990.
19. Дробышевич В. И., Дымников В. Я, Ривин Г. С. Задачи
по вычислительной математике. — М.: Наука, 1980.
20. Лабораторный практикум по курсу "Основы
вычислительной математики". — М.: МЗ-Пресс, 2001.
21. Плис А. И., Сливина Я. А. Лабораторный практикум по
высшей математике. — М.: Высшая школа, 1994.
22. Ракитин В. И., Первушин В. Е. Практическое
руководство по методам вычислений с приложением программ для
персональных компьютеров. — М.: Высшая школа, 1998.
23. Сборник задач по методам вычислений / Под редакцией
Я. И. Монастырского — М.: Наука, 1994.

184
ОГЛАВЛЕНИЕ
Лекция 1. Вычислительный эксперимент 3
Лекция 2. Погрешности чисел. Значащие цифры. Округление ... 9
Лекция 3. Верные значащие цифры. Погрешности действий .... 15
Лекция 4. Матрица: типы, действия, характеристики 19
Лекция 5. Матрицы: разбиение, разложение (факторизация) .... 27
Лекция 6. Обратная матрица. Обращение матриц 32
Лекция 7. Характеристики матрицы СЛАУ 39
Лекция 8. Прямые методы решения СЛАУ 47
Лекция 9. Решение СЛАУ (методы итерации) 57
Лекция 10. Решение специальных СЛАУ 64
Лекция 11. Решение нелинейных уравнений (метод секущих) .... 73
Лекция 12. Решение нелинейных уравнений (метод касательных) . .81
Лекция 13. Разности. Разностные уравнения 87
Лекция 14. Приближение функций (полином Ньютона) 93
Лекция 15. Приближение функций (полином Лагранжа) 100
Лекция 16. Приближение функций (сплайн-интерполяция) .... 107
Лекция 17. Точечное квадратичное приближение функций .... 116
Лекция 18. Интегральное квадратичное приближение 123
Лекция 19. Численное дифференцирование 127
Лекция 20. Численное интегрирование (метод прямоугольников) . 133
Лекция 21. Численное интегрирование (метод трапеций) 139
Лекция 22. Численное интегрирование (метод парабол) 144
Лекция 23. Численное вычисление несобственных интегралов ... 149
Лекция 24. Задача Коши ОДУ. Численные методы Пикара и Эйлера 154
Лекция 25. Численное решение краевой задачи ОДУ 161
Лекции 26. Численное решение уравнений в частных производных 166
Приложение. Методы построения аналитических моделей .... 173
П.1. Пример физической модели 173
П.2. Разностные уравнения 177
П.З. Замена переменных функции 178
П.4. Программа решения МКТРЭ двумерного стационарного
уравнения теплопроводности 180
Литература 182

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