Л.И. ТУРЧАК, П.В. ПЛОТНИКОВ
основы
ЧИСЛЕННЫХ
МЕТОДОВ
ИЗДАНИЕ ВТОРОЕ, ПЕРЕРАБОТАННОЕ И ДОПОЛНЕННОЕ
Допущено Министерством образования
Российской Федерации в качестве учебного пособия
для студентов высших учебных заведений
МОСКВА
ФИЗМАТЛИТ
2003

УДК 519.6
ББК 22.19
Т89
Турчак Л. И., Плотников П. В. Основы численных методов: Учебное
пособие. — 2-е изд., перераб. и доп. — М.: ФИЗМАТЛИТ, 2003. — 304 с. — ISBN
5-9221-0153-6.
Содержит основные сведения о численных методах, необходимые для перво-
первоначального знакомства с предметом. Излагаются основы численных методов для
систем линейных и нелинейных уравнений, а также дифференциальных и инте-
интегральных уравнений. Имеется много задач, примеров и алгоритмов для облегчения
понимания логической структуры рассматриваемых методов и их использования в
расчетах на компьютерах.
Первое издание — 1987 г.
Для студентов вузов.
Табл. 21. Ил. 83. Библиогр. 66 назв.
© ФИЗМАТЛИТ, 2002, 2003
ISBN 5-9221-0153-6	© Л.И. Турчак, П.В. Плотников, 2002, 2003

ОГЛАВЛЕНИЕ
Предисловие 	 7
Введение	 10
1 Этапы решения задачи на компьютере A0). 2 Математические модели A2).
3 Численные методы A3).
ГЛАВА 1
ТОЧНОСТЬ ВЫЧИСЛИТЕЛЬНОГО ЭКСПЕРИМЕНТА
§ 1. Приближенные числа 	 15
1 Числа с плавающей точкой A5). 2 Понятие погрешности A7). 3 Действия над
приближенными числами A8).
§ 2. Погрешности вычислений 	 20
1 Источники погрешностей B0). 2 Уменьшение погрешностей B1). 3 О реше-
решении квадратного уравнения B3).
§ 3. Устойчивость. Корректность. Сходимость 	 26
1 Устойчивость B6). 2 Корректность B7). 3 Неустойчивость методов B8).
4	Понятие сходимости B9). Упражнения B9).
ГЛАВА 2
АППРОКСИМАЦИЯ ФУНКЦИЙ
§ 1. Понятие о приближении функций 	 31
1 Постановка задачи C1). 2 Точечная аппроксимация C2). 3 Непрерывная ап-
аппроксимация. Равномерное приближение C4). 4 Вычисление многочленов C5).
§ 2. Использование рядов	 36
1 Элементарные функции C6). 2 Многочлены Чебышева C9). 3 Рациональные
приближения D4).
§ 3. Интерполирование	 47
1 Линейная и квадратичная интерполяции D7). 2 Многочлен Лагранжа D9).
3 Многочлен Ньютона E0). 4 Точность интерполяции E4). 5 Сплайны E5).
6 О других формулах интерполяции E8). 7 Функции двух переменных E9).
§ 4. Подбор эмпирических формул	 60
1 Характер опытных данных F0). 2 Эмпирические формулы F1). 3 Определение
параметров эмпирической зависимости F3). 4 Метод наименьших квадратов F6).
5	Локальное сглаживание данных F9). Упражнения G0).

Оглавление
ГЛАВА 3
ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ
\ 1. Численное дифференцирование	 72
1 Аппроксимация производных G2). 2 Погрешность численного дифференциро-
дифференцирования G3). 3 Использование интерполяционных формул G5). 4 Метод неопре-
неопределенных коэффициентов G9). 5 Улучшение аппроксимации (81). 6 Частные
производные (82).
\ 2. Численное интегрирование	 85
1 Вводные замечания (85). 2 Методы прямоугольников и трапеций (88). 3 Ме-
Метод Симпсона (91). 4 Использование сплайнов (93). 5 Погрешность численного
интегрирования (94). 6 Адаптивные алгоритмы (97). 7 О других методах. Осо-
Особые случаи A00). 8 Кратные интегралы A02). 9 Метод Монте-Карло A04).
Упражнения A06).
ГЛАВА 4
СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ
\ 1. Основные понятия	 107
1 Линейные системы A07). 2 О методах решения линейных систем (ПО). 3 Дру-
Другие задачи линейной алгебры A11).
\ 2. Прямые методы	 113
1 Вводные замечания A13). 2 Метод Гаусса A14). 3 Определитель и обратная
матрица A20). 4 Метод прогонки A21). 5 О других прямых методах A23).
\ 3. Итерационные методы	 124
1 Уточнение решения A24). 2 Метод простой итерации A26). 3 Метод Гаусса-
Зейделя A27).
\ 4. Задачи на собственные значения 	 131
1 Основные понятия A31). 2 Метод вращений A35). 3 Трехдиагональные мат-
матрицы A39). 4 Частичная проблема собственных значений A41). Упражнения A43).
ГЛАВА 5
НЕЛИНЕЙНЫЕ УРАВНЕНИЯ
\ 1. Уравнения с одним неизвестным 	 145
1 Вводные замечания A45). 2 Метод деления отрезка пополам (метод бисек-
ции). A46). 3 Метод хорд A48). 4 Метод Ньютона (метод касательных). A49).
5 Метод простой итерации A51).
\2. О решении алгебраических уравнений 	 152
1 Действительные корни A52). 2 Комплексные корни A53).
\ 3. Системы уравнений 	 154
1 Вводные замечания A54). 2 Метод простой итерации и метод Зейделя A55).
3 Метод Ньютона A55). Упражнения A58).
ГЛАВА 6
МЕТОДЫ ОПТИМИЗАЦИИ
\ 1. Основные понятия 	 160
1 Определения A60). 2 Задачи оптимизации A61). 3 Пример постановки зада-
задачи A62).

Оглавление
§2. Одномерная оптимизация 	 162
1 Задачи на экстремум A62). 2 Методы поиска A64). 3 Метод золотого сече-
сечения A66). 4 Метод Ньютона A70).
§ 3. Многомерные задачи оптимизации 	 172
1 Минимум функции нескольких переменных A72). 2 Метод покоординатного
спуска A74). 3 Метод градиентного спуска A76).
§ 4. Задачи с ограничениями	 178
1 Метод штрафных функций. A78). 2 Линейное программирование A80). 3 Ге-
Геометрический метод A83). 4 Симплекс-метод A86). 5 Задача о ресурсах A89).
Упражнения A92).
ГЛАВА 7
ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ
§ 1. Основные понятия 	 194
1 Постановка задач A94). 2 О методах решения A96). 3 Разностные мето-
методы A98).
§2. Задача Коши 	 201
1 Общие сведения B01). 2 Метод Эйлера B02). 3 Модификации метода Эйле-
Эйлера B05). 4 Методы Рунге-Кутта B07). 5 Многошаговые методы B10). 6 По-
Повышение точности результатов B12).
§ 3. Краевые задачи	 214
1 Предварительные замечания B14). 2 Метод стрельбы B16). 3 Методы конеч-
конечных разностей B18). Упражнения B22).
ГЛАВА 8
УРАВНЕНИЯ С ЧАСТНЫМИ ПРОИЗВОДНЫМИ
§ 1. Элементы теории разностных схем	 224
1 Вводные замечания B24). 2 О построении разностных схем B26). 3 Сходи-
Сходимость. Аппроксимация. Устойчивость B30).
§ 2. Уравнения первого порядка 	 236
1 Линейное уравнение переноса B36). 2 Квазилинейное уравнение. Разрывные
решения B44). 3 Консервативные схемы B50). 4 Системы уравнений. Характе-
Характеристики B51).
§ 3. Уравнения второго порядка 	 254
1 Волновое уравнение B54). 2 Уравнение теплопроводности B58). 3 Понятие
о схемах расщепления B62). 4 Уравнение Лапласа B65). Упражнения B69).
ГЛАВА 9
ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ
§ 4. Постановка задач	 271
1 Вводные замечания B71). 2 Виды интегральных уравнений B72).
§ 5. Методы решения	 273
1 Методы последовательных приближений B73). 2 Численные методы B75).

6	Оглавление
§ 6. Сингулярные уравнения		278
1 Сингулярные интегралы B78). 2 Численное решение сингулярных интеграль-
интегральных уравнений B82). Упражнения B85).
Приложение А. Структурограммы 		286
Приложение Б. Многочлены Чебышева 		288
Литература		290
Предметный указатель 		293

ПРЕДИСЛОВИЕ
Внедрение компьютеров во все сферы человеческой деятельности тре-
требует от специалистов разного профиля овладения навыками использования
вычислительной техники. Повышается уровень подготовки студентов ву-
вузов, которые уже с первых курсов приобщаются к использованию компью-
компьютеров и простейших численных методов, не говоря уже о том, что при
выполнении курсовых и дипломных работ применение вычислительной
техники стало нормой.
Применение компьютеров приобрело сейчас массовый характер. Они
используются не только при научных и инженерных расчетах, но и для
хранения и обработки информации, при решении ряда других задач и даже
в быту. Тем не менее использование компьютера для проведения математи-
математических вычислений не потеряло своей актуальности. Причем оно распро-
распространилось не только на точные, технические и экономические науки, но
и на такие традиционно нематематические специальности, как медицина,
лингвистика, психология и др. Возникла многочисленная категория спе-
специалистов — пользователей компьютеров, использующих их в качестве
вычислительного инструмента и поэтому нуждающихся в литературе по
соответствующим дисциплинам.
Одной из основных дисциплин является вычислительная математика.
Она изучает методы построения и исследования численных методов реше-
решения математических задач, которые моделируют различные процессы.
Численные методы разрабатывают и исследуют, как правило, высоко-
высококвалифицированные специалисты - математики. Что касается подавляющей
части студентов нематематических специальностей и инженерно - техни-
технических работников, то для них главным является понимание основных идей,
методов, особенностей и областей их применения. Следует также иметь в
виду, что указанная категория читателей не обладает достаточными мате-
математическими знаниями для подробного исследования численных методов.
К тому же в этом нет особой необходимости специалисту - нематематику,
использующему численные методы как готовый инструмент в своей прак-
практической работе.
В предлагаемом учебном пособии в сжатом виде приводятся основные
необходимые сведения о численных методах решения различных приклад-
прикладных задач. Изложение проводится на доступном для студентов втуза уровне.
При необходимости напоминаются основные сведения из курса высшей

Предисловие
математики. Для многих рассматриваемых методов приводятся алгоритмы,
а также примеры решения задач, способствующие лучшему пониманию
материала. Книга написана с учетом особенностей применения численных
методов при решении задач с использованием компьютеров.
Поскольку данное учебное пособие не ориентировано на студентов
конкретной специальности, то приведенные в нем задачи носят общий
характер. Большой выбор интересных задач содержится в книгах при-
прикладного характера, включенных в список литературы. Они, в частно-
частности, могут быть использованы при выполнении курсовых и дипломных
работ, а также в научно-исследовательской работе студентов. В список
литературы включены также некоторые пособия по численным методам,
которые авторы использовали в работе над данным пособием. Читатель
может найти в них более подробные сведения по интересующим его
разделам курса. Разумеется, это далеко не полный перечень литературы
по численным методам и их приложениям.
При изложении материала сказался стиль чтения авторами курсов лек-
лекций для студентов нематематических специальностей втузов и слушателей
факультета повышения квалификации. Книга будет полезна студентам и
специалистам при первоначальном знакомстве с предметом. Она может
служить кратким справочным пособием, которое студенты могут исполь-
использовать при выполнении расчетных заданий.
Книга содержит девять глав. После каждой главы приведены упраж-
упражнения для самостоятельного решения. Их выполнение читателем будет
способствовать лучшему усвоению материала. Упражнения повышенной
трудности отмечены звездочкой.
В главе 1 излагаются основные понятия, связанные с погрешностями
вычислений. Рассматриваются источники погрешностей при расчетах на
компьютерах.
Глава 2 посвящена различным способам аппроксимации (приближения)
функций. При рассмотрении интерполирования дано понятие сплайнов,
которые получили широкое распространение в вычислительной практике.
Выписаны некоторые формулы, которые могут быть полезными при само-
самостоятельной работе.
Вопросы численного дифференцирования и численного интегрирования
изложены в главе 3. Здесь же приведены выражения для аппроксимаций
производных, которые могут быть использованы при построении разност-
разностных схем для решения дифференциальных уравнений. Среди методов чис-
численного интегрирования упомянуто использование сплайнов. Приведено
также понятие адаптивных алгоритмов, которые сейчас широко исполь-
используются и при решении других задач.
Глава 4 содержит основные сведения по численному решению задач
линейной алгебры. При первоначальном знакомстве можно опустить § 4,
в котором излагаются некоторые задачи на собственные значения, по-
поскольку эта тема носит специальный характер.
В главе 5 изложены основные методы решения нелинейных уравнений
(алгебраических и трансцендентных) и их систем.

Предисловие
Глава 6, посвященная методам решения задач оптимизации, содержит
также элементы линейного программирования.
Методы решения задач Коши и краевых задач для обыкновенных диф-
дифференциальных уравнений излагаются в главе 7.
В главе 8 излагаются численные методы решения уравнений с частными
производными и приводятся некоторые элементы теории разностных схем.
Глава 9, посвященная интегральным уравнениям, носит ознакомитель-
ознакомительный характер, и при первом чтении может быть опущена. Вместе с тем
следует отметить, что решение интегральных уравнений, в том числе и
сингулярных, необходимо во многих областях науки (механике, физике и др,).
При подготовке ко второму изданию материал книги был перерабо-
переработан, дополнен изложением некоторых новых вопросов. Были исправлены
замеченные опечатки, а также добавлена часть упражнений.
Авторы искренне признательны академику О. М. Белоцерковскому
за ценные замечания по рукописи. Полезные предложения по улучше-
улучшению содержания высказали 3. С. Волк, Ю. Г. Евтушенко, И. К. Лифанов,
В. Б. Миносцев, Г. П. Тиняков, Э. Г. Шифрин и другие товарищи, прочитав-
прочитавшие рукопись или отдельные ее части. Большую помощь в работе над кни-
книгой оказал В. В. Щенников. Всем им авторы выражают свою глубокую бла-
благодарность.

ВВЕДЕНИЕ
1. Этапы решения задачи на компьютере. Вычислительная техни-
техника нашла эффективное применение при проведении трудоемких расчетов
в научных исследованиях. Действительно, современные компьютеры за
одну секунду выполняют такой объем вычислений, на который человеку
не хватит всей жизни.
При решении задачи на компьютере основная роль все-таки принад-
принадлежит человеку. Машина лишь выполняет его задания по разработанной
программе. Роль человека и машины легко уяснить, если процесс решения
задачи разбить на следующие этапы.
Постановка задачи. Этот этап заключается в содержательной
(физической) постановке задачи и определении конечных целей решения.
Построение математической модели (математическая
формулировка задачи). Модель должна правильно (адекватно) описывать
основные законы физического процесса. Построение или выбор математи-
математической модели из существующих требует глубокого понимания проблемы
и знания соответствующих разделов математики.
Разработка численного метода. Поскольку компьютер
может выполнять лишь простейшие операции, он «не понимает» постанов-
постановки задачи даже в математической формулировке. Для ее решения должен
быть найден численный метод, позволяющий свести задачу к некоторому
вычислительному алгоритму. Разработкой численных методов занимаются
специалисты в области вычислительной математики. Специалисту - при-
прикладнику для решения задачи, как правило, необходимо из имеющегося
арсенала методов выбрать тот, который наиболее пригоден в данном кон-
конкретном случае.
Разработка алгоритма. Процесс решения задачи (вычисли-
(вычислительный процесс) записывается в виде последовательности элементарных
арифметических и логических операций, приводящей к конечному резуль-
результату и называемой алгоритмом решения задачи. Алгоритм можно наглядно
изобразить в виде блок-схемы, структурограммы и т. п. Опытный вычисли-
вычислитель зачастую может и не прибегать к такому наглядному представлению
алгоритма, непосредственно переходя к следующему этапу.
Программирование. Алгоритм решения задачи записывается на
понятном машине языке в виде точно определенной последовательности

Введение	11
операций — программы для компьютера. Составление программы (про-
(программирование) обычно производится с помощью некоторого промежу-
промежуточного (алгоритмического) языка, а ее трансляция (перевод на машинный
язык) осуществляется самой вычислительной системой.
Отладка программы. Составленная программа содержит разного
рода ошибки, неточности, описки. Отладка программы на машине включает
контроль программы, диагностику (поиск и определение содержания) оши-
ошибок, их исправление. Программа испытывается на решении контрольных
(тестовых) задач для получения уверенности в достоверности результатов.
Проведение расчетов. На этом этапе готовятся исходные
данные для расчетов и проводится счет по отлаженной программе. При
этом для уменьшения ручного труда по обработке результатов желательно
использовать удобные формы выдачи результатов, особенно их графическое
представление {визуализацию).
Анализ результатов. Результаты расчетов анализируются,
оформляется научно-техническая документация.
Если при решении конкретной задачи возможно использование уже
имеющихся прикладных программных средств, то некоторые из перечис-
перечисленных этапов могут быть опущены. Так, для решения многих (хотя и
достаточно узких) классов задач созданы программные продукты, суще-
существенно облегчающие труд вычислителя. Речь может идти, например, о том,
что вычислитель сообщает программе только математическую модель (или
даже постановку задачи) и исходные данные, а выбор метода, проведение
расчетов, выдачу результатов программа берет на себя. Но даже в этом слу-
случае нельзя забывать о том, что полученное решение обычно является лишь
приближенным, что каждая модель и каждый метод имеют свои области
применимости. Следовательно, специалисту, использующему компьютер
для решения прикладных задач, необходимо иметь представление об осно-
основах математического моделирования, численных методов, о возможностях
компьютеров, уметь анализировать полученные результаты с точки зрения
их достоверности.
Следует отметить еще один важный момент в процессе решения задачи
с помощью компьютера. Это — экономичность выбранного способа ре-
решения задачи, численного метода, модели компьютера. В частности, если
задача допускает простое аналитическое решение или измерение, то вряд
ли целесообразно делать вычисления на машине. Иногда решение зада-
задачи производят с помощью большого вычислительного комплекса, хотя это
можно было осуществить с использованием персонального компьютера.
Не умаляя значения физического эксперимента, нужно все-таки отме-
отметить неуклонно возрастающую долю вычислений на компьютере в общем
объеме решения научно-технических задач. В связи с этим наряду с увели-
увеличением парка вычислительных машин и повышением их «интеллектуаль-
«интеллектуальных» возможностей возрастает интерес к математическому моделированию
и разработке численных методов.

12	Введение
2. Математические модели. Основное требование, предъявляемое
к математической модели, — адекватность рассматриваемому явлению,
т. е. она должна достаточно точно (в рамках допустимых погрешностей)
отражать характерные черты явления. Вместе с тем она должна обладать
сравнительной простотой и доступностью исследования.
Приведем примеры некоторых математических моделей, оказавших
огромное влияние на развитие различных отраслей науки и техники.
При построении математических моделей получают некоторые математи-
математические соотношения (как правило, уравнения).
Пример. Пусть в начальный момент времени t = 0 тело, нахо-
находящееся на высоте /го, начинает двигаться вертикально вниз с началь-
начальной скоростью vq. Требуется найти закон движения тела, т. е. построить
математическую модель, которая позволила бы математически описать
данную задачу и определить параметры движения в любой момент вре-
времени.
Прежде чем строить указанную модель, нужно принять некоторые
допущения, если они не заданы. В частности, предположим, что данное
тело обладает средней плотностью, значительно превышающей плот-
плотность воздуха, а его форма близка к шару. В этом случае можно пре-
пренебречь сопротивлением воздуха и рассматривать свободное падение те-
тела с учетом ускорения д. Соответствующие соотношения для высоты h
и скорости v в любой момент времени t хорошо известны из школьного
курса физики. Они имеют вид
h = ho-vot	—, v = vQ+gt.	@.1)
Эти формулы являются искомой математической моделью свободного
падения тела. Область применения данной модели ограничена случаями,
в которых можно пренебречь сопротивлением воздуха. Во многих зада-
задачах о движении тел в атмосфере планеты модель @.1) не может быть
использована, поскольку при ее применении мы получили бы неверный
результат. К таким задачам относятся движение капли, вход в атмосферу
тел малой плотности, спуск на парашюте и др. Здесь необходимо постро-
построить более точную математическую модель, учитывающую сопротивление
воздуха. Если обозначить через F(t) силу сопротивления, действующую
на тело массой га, то его движение можно описать с помощью уравнений
^	^ = ~v.	@.2)
К этой системе уравнений необходимо добавить начальные условия при
t = 0:
v = г>о, h = Hq.	@.3)
Соотношения @.2) и @.3) являются математической моделью для
задачи движения тела в атмосфере. Существуют и другие, более сложные

Введение	13
модели подобных задач (например, задача о движении планера). Заметим
также, что модель @.1) легко получается из @.2), @.3) при F = 0.
Известно большое число математических моделей различных про-
процессов или явлений. Укажем некоторые из них, широко используемые
в механике. Модель абсолютно твердого тела позволила получить урав-
уравнения движения тел в динамике полета. Модель идеального газа приве-
привела к системе уравнений Эйлера, описывающей невязкие потоки газов.
В гидродинамике широко известна модель на основе уравнений Навье -
Стокса, в кинетической теории газов — уравнения Больцмана. В меха-
механике деформируемого твердого тела известны математические модели,
описывающие различные среды (упругую, упруго - пластичную и др.).
Имеются математические модели и для описания задач экономики,
социологии, медицины, лингвистики и др.
Адекватность и сравнительная простота модели не исчерпывают
предъявляемых к ней требований. Обратим еще внимание на необходи-
необходимость правильной оценки области применимости математической моде-
модели. Например, модель свободно падающего тела, в которой пренебрегают
сопротивлением воздуха, весьма эффективна для твердых тел с большой
средней плотностью и формой поверхности, близкой к сферической. Вме-
Вместе с тем в ряде других случаев (движения капельки жидкости, парашют-
парашютного устройства и др.) для решения задачи уже недостаточно известных
из курса физики простейших формул. Здесь необходимы более сложные
математические модели, учитывающие сопротивление воздуха и другие
факторы.
Отметим, что успех решения задачи в значительной степени опреде-
определяется выбором математической модели: здесь в первую очередь нужны
глубокие знания в той области, к которой принадлежит поставленная
задача. Кроме того, необходимы знания соответствующих разделов мате-
математики и возможностей компьютеров.
3. Численные методы. С помощью математического моделирова-
моделирования решение научно-технической задачи сводится к решению матема-
математической задачи, являющейся ее моделью. Для решения математических
задач используются следующие основные группы методов: аналитиче-
аналитические, графические и численные.
При использовании аналитических методов решение задачи удается
выразить с помощью формул. В частности, если математическая задача
состоит в решении простейших алгебраических или трансцендентных
уравнений, дифференциальных уравнений и т. п., то использование из-
известных из курса математики приемов сразу приводит к цели. К сожале-
сожалению, на практике это бывает достаточно редко.
Графические методы позволяют в ряде случаев оценить порядок ис-
искомой величины. Основная идея этих методов состоит в том, что реше-
решение находится путем геометрических построений. Например, для нахож-
нахождения корней уравнения f(x) = 0 строится график функции у = f(x),
точки пересечения которого с осью абсцисс и будут искомыми корнями.

14	Введение
Графические методы могут применяться для получения начального
приближения к решению, которое затем уточняется с помощью числен-
численных методов.
Основным инструментом для решения сложных математических задач
в настоящее время являются численные методы, позволяющие свести ре-
решение задачи к выполнению конечного числа арифметических действий
над числами; при этом результаты получаются в виде числовых значений.
Подчеркнем важные отличия численных методов от аналитических.
Во-первых, численные методы позволяют получить лишь приближенное
решение задачи. Во-вторых, они обычно позволяют получить лишь реше-
решение задачи с конкретными значениями параметров и исходных данных.
Поясним второе отличие на примере. По формулам @.1) (по аналити-
аналитическому решению) можно проанализировать, как изменится закон дви-
движения при изменении параметра g и начальных значений vq, ho. Если
в модели @.2), @.3) выражение для F(t) имеет простой вид (напри-
(например, F = const), то можно получить аналитическое решение, аналогич-
аналогичное @.1). Это решение тоже легко исследовать на предмет зависимости
от изменения параметров и начальных условий. Если же выражение для
F(t) достаточно сложно, то задачу @.2), @.3) можно решить только чис-
численно. При этом вместо общей формулы решения в результате расчета
будут получены значения v и h для некоторого набора моментов време-
времени t при конкретных значениях д, га, vq, ho. Для получения решения
при других значениях параметров и (или) других начальных условиях
необходимо провести новый расчет. Для анализа зависимости решения
от параметров и начальных условий необходима большая серия расчетов.
Несмотря на эти недостатки, численные методы незаменимы в слож-
сложных задачах, которые не допускают аналитического решения.
Многие численные методы разработаны давно, однако при вычисле-
вычислениях вручную они могли использоваться лишь для решения не слиш-
слишком трудоемких задач. С появлением компьютеров начался период бур-
бурного развития численных методов и их внедрения в практику. Только
вычислительной машине под силу выполнить за короткое время объем
вычислений в миллиарды, триллионы и более операций, необходимых
для решения многих современных задач.
Численный метод наряду с возможностью получения результата за
приемлемое время должен обладать и еще одним важным качеством —
не вносить в вычислительный процесс значительных погрешностей.

ГЛАВА 1
ТОЧНОСТЬ ВЫЧИСЛИТЕЛЬНОГО
ЭКСПЕРИМЕНТА
§ 1. Приближенные числа
1. Числа с плавающей точкой. Числа могут быть представлены в
памяти компьютеров различными способами. Современные компьютеры
(процессоры), как правило, позволяют обрабатывать целые числа, а также
дробные числа в форме с плавающей точкой 1\
Как известно, множество целых чисел бесконечно. Однако процессор из-
за ограниченности его разрядной сетки может оперировать лишь с некото-
некоторым конечным подмножеством этого множества. В современных компью-
компьютерах для хранения целого числа обычно отводится 4 байта памяти 2\
что позволяет представлять целые числа, находящиеся примерно в диапа-
диапазоне от -2 • 109 до 2-109.
При решении научно-технических задач в основном используются дей-
действительные (вещественные) числа. В компьютерах они представляются
в форме с плавающей точкой. Десятичное число D в этой форме записи
имеет вид D = ±га • 10п, где тип — соответственно мантисса числа и
его порядок. Например, число -273.9 можно записать в виде: -2739 • 1СГ1,
—2.739 • 102, —0.2739 • 103. Последняя запись — нормализованная форма
числа с плавающей точкой. Таким образом, если представить мантиссу
числа в виде т = 0. &\&i... d&, то при d\ ф 0 получим нормализованную
форму числа с плавающей точкой. В дальнейшем, говоря о числах с пла-
плавающей точкой, будем иметь в виду именно эту форму. Обычная же запись
числа в виде —273.9 называется формой записи с фиксированной точкой.
В настоящее время такое представление используется в компьютерах, как
правило, только на этапе ввода и вывода чисел.
Все сказанное выше о числах с плавающей точкой распространяется и
на числа, записанные в других системах счисления. Число А в системе счис-
счисления с основанием а можно представить в виде А = ±0. а\п2 ... ctk • ап,
где ai, а2,..., ctk — целые числа из диапазона 0,..., а — 1. Из этой записи
следует, что подмножество действительных чисел, с которым оперирует
1 ' В англоязычных странах целую и дробную части при десятичной записи числа
разделяют точкой, а не запятой. Часто аналогично поступают и в специализирован-
специализированной литературе на русском языке.
2^ Здесь и далее при упоминании характеристик компьютеров имеются в ви-
виду прежде всего широко распространенные персональные компьютеры на базе
процессоров фирмы Intel.

16
Гл. 1. Точность вычислительного эксперимента
конкретный компьютер, не является бесконечным: оно конечно и опреде-
определяется разрядностью &, а также границами порядка nb n2 (ni ^ п ^ п2).
Можно показать, что это подмножество содержит
N = 2(а - 1)(п2 - ni +
A.1)
чисел, наименьшим и наибольшим по модулю являются соответственно
числа
и
= A-а
A.2)
называемые машинным нулем и машинной бесконечностью.
Границы порядка rii, n2 определяют ограниченность действительных
чисел по величине, а разрядность к — дискретность распределения их
на отрезке числовой оси. Например, в случае десятичных чисел при че-
четырехразрядном представлении все значения, находящиеся в промежутке
@.28505,0.28515), представляются числом 0.2851 (при выполнении округ-
округления). Если к этому числу 0.2851 прибавить число, меньшее по модулю
половины единицы последнего разряда (т. е. меньшее по модулю 0.00005),
в результате получится то же самое число 0.2851.
В настоящее время большинство производителей процессоров в ос-
основном придерживаются стандарта IEEE 754 ^ для арифметических опе-
операций над двоичными числами с плавающей точкой. Данный стандарт пре-
предусматривает наличие, в частности, двух двоичных (а = 2) форматов:
с одинарной точностью и с двойной точностью. Приведем для этих форма-
форматов размер отводимой памяти, значения к, nb n2 и приближенные значения
Мо и Мое. Заметим, что стандарт IEEE 754 предусматривает обработку
чисел, меньших по модулю Мо, но не меньших Mq*, правда, с меньшей
разрядностью к.
Точность
Одинарная
Двойная
Байты
4
8
к
24
53
пг
-125
-1021
П2
128
1024
Мо
1.2-108
2.2-Ю-308
м0*
1.4-Ю-45
4.9-Ю-324
Мое
3.4 -1038
1.8-10308
Поскольку для человека более удобной является десятичная система
счисления, возникает вопрос о том, скольким десятичным разрядам соот-
соответствует указанная двоичная разрядность к. Можно считать, что к соот-
соответствует 6-9 десятичным разрядам при одинарной и 15-17 разрядам при
двойной точности.
В современных языках программирования предусмотрены типы дан-
данных для представления вещественных чисел с одинарной и двойной
точностью. Например, в языке Си это типы float и double, в языке Па-
Паскаль — single и double, в языке Фортран — real и double precision. Обычно
эти представления соответствуют стандарту IEEE 754.
IEEE — институт инженеров по электротехнике и электронике (США).

§ 1. Приближенные числа	17
Таким образом, компьютер оперирует с приближенными значениями
действительных чисел. Мерой точности приближенных чисел является по-
погрешность.
2. Понятие погрешности. Различают два вида погрешностей —
абсолютную и относительную. Абсолютная погрешность некоторого числа
равна разности между его истинным значением и приближенным значени-
значением, полученным в результате вычисления или измерения. Относительная
погрешность — это отношение абсолютной погрешности к приближенно-
приближенному значению числа.
Таким образом, если а — приближенное значение числа х, то выра-
выражения для абсолютной и относительной погрешностей запишутся соответ-
соответственно в виде
Ах = х - а, 5х = Ах/а.
К сожалению, истинное значение величины х обычно неизвестно. По-
Поэтому приведенные выражения для погрешностей практически не могут
быть использованы. Имеется лишь приближенное значение а и нужно най-
найти его предельную погрешность Аа, являющуюся верхней оценкой модуля
абсолютной погрешности, т. е. |Дж| ^ Аа. В дальнейшем значение Аа
принимается в качестве абсолютной погрешности приближенного числа а.
В этом случае истинное значение х находится в интервале (а — Аа, а + Аа).
Для приближенного числа, полученного в результате округления,
абсолютная погрешность Аа принимается равной половине единицы по-
последнего разряда числа. Например, значение а = 0.734 могло быть полу-
получено округлением чисел 0.73441, 0.73353 и др. При этом |Дж| ^ 0.0005,
и полагаем Аа = 0.0005. Если при вычислениях на компьютере округление
не производится, а цифры, выходящие за разрядную сетку машины, отбра-
отбрасываются, то максимально возможная погрешность результата выполнения
операции в два раза больше по сравнению со случаем округления.
Приведем примеры оценки абсолютной погрешности при некоторых
значениях приближенной величины а:
а
Аа
51.7
0.05
-0.0031
0.00005
16
0.5
16.00
0.005
Предельное значение относительной погрешности — отношение пре-
предельной абсолютной погрешности к абсолютной величине приближенного
числа:
да = Да/|а|.
Например, ?(-2.3) = 0.05/2.3 и 0.022 B.2%). Заметим, что погрешность
округляется всегда в сторону увеличения. В данном случае ?(—2.3) и 0.03.
Приведенные оценки погрешностей приближенных чисел справедливы,
если в записи этих чисел все значащие цифры верные. Напомним, что зна-
значащими цифрами считаются все цифры данного числа, начиная с первой
ненулевой цифры. Например, в числе 0.037 две значащие цифры: 3 и 7,
2 Л.И. Турчак, П.В. Плотников

18	Гл. 1. Точность вычислительного эксперимента
а в числе 14.80 все четыре цифры значащие. Кроме того, при измене-
изменении формы записи числа (например, при записи в форме с плавающей
точкой) число значащих цифр не должно меняться, т. е. нужно соблюдать
равносильность преобразований. Например, записи 7500 = 0.7500 • 104
и 0.110 • 102 = 11.0 равносильные, а записи 7500 = 0.75 • 104 и 0.110 •
• 102 = 11 неравносильные.
3. Действия над приближенными числами. Сформулируем пра-
правила оценки предельных погрешностей при выполнении операций над
приближенными числами.
При сложении или вычитании чисел их абсолютные погрешности скла-
складываются. При умножении или делении чисел друг на друга их относи-
относительные погрешности складываются. При возведении в степень прибли-
приближенного числа его относительная погрешность умножается на показатель
степени.
Для случая двух приближенных чисел а и b эти правила можно запи-
записать в виде формул
А(а ±Ь) = Аа + ДЬ, 5(а • Ь) = 6а + 5Ъ,
5(а/Ъ) = 5а + 5b, 5(ak) = кба.
Относительная погрешность суммы положительных слагаемых заклю-
заключена между наибольшим и наименьшим значениями относительных по-
погрешностей этих слагаемых. Действительно, пусть а > 0, b > 0, га =
= min(Ea, 5b), M = maxEa, 5b). Тогда
., _, А(а + Ъ) Аа + АЪ а5а + Ъ6Ъ аМ + ЪМ Ъ/Г
5(а + Ъ) = —	-^ =	— =	— ^	— = М.
а + о	а + о	а + о	а + о
Аналогично, 6(а + Ь) ^ га. На практике для оценки погрешности прини-
принимается наибольшее значение М.
Пример 1. Найти относительную погрешность функции
а + 6
У =
х3A — х)
Используя формулы A.3), получаем
v	п 2 \ \п + /I	т1 11 - т
Ll^l	I I
Полученная оценка относительной погрешности содержит в знаменате-
знаменателе выражение |1 — х\. Ясно, что при х и 1 можно получить очень
большую погрешность. В связи с этим рассмотрим подробнее случай
вычитания близких чисел.

§ 1. Приближенные числа	19
Запишем выражение для относительной погрешности разности двух
чисел в виде
., _, А(а-Ъ) Аа + АЬ
При а и 6 эта погрешность может быть сколь угодно большой.
П р и м е р 2. Пусть а = 2520, Ь = 2518. В этом случае имеем аб-
абсолютные погрешности исходных данных А а = А 6 = 0.5 и относительные
погрешности да и 56 = 0.5/2518 и 0.0002 @.02 %). Относительная пог-
погрешность разности равна
Следовательно, при малых погрешностях в исходных данных мы по-
получили весьма неточный результат. Нетрудно подсчитать, что даже при
случайных изменениях а и Ь на единицу в последних разрядах их разность
может принимать значения 0,1, 2, 3, 4. Поэтому при организации вычисли-
вычислительных алгоритмов следует избегать вычитания близких чисел; при воз-
возможности алгоритм нужно видоизменить во избежание потери точности
на некотором этапе вычислений.
Из рассмотренных правил следует, что при сложении или вычитании
приближенных чисел желательно, чтобы эти числа обладали одинаковыми
абсолютными погрешностями, т. е. одинаковым числом разрядов после де-
десятичной точки. Например, 38.723 + 4.9 = 43.6; 425.4-0.047 = 425.4.
Учет отброшенных разрядов не повысит точность результатов. При умно-
умножении и делении приближенных чисел количество значащих цифр вырав-
выравнивается по наименьшему из них.
Наряду с приведенными выше оценками погрешностей при выполне-
выполнении некоторых операций над приближенными числами можно записать
аналогичные оценки и для вычисления функций, аргументами которых яв-
являются приближенные числа. Однако более полным оказывается общее
правило, основанное на вычислении приращения (погрешности) функции
при заданных приращениях (погрешностях) аргументов.
Рассмотрим функцию одной переменной у = f(x).. Пусть а — прибли-
приближенное значение аргумента х, Аа — его абсолютная погрешность. Абсо-
Абсолютную погрешность функции можно считать ее приращением, которое
она испытывает при изменении аргумента на Аа. Это приращение мож-
можно заменить дифференциалом: Ау и dy. Тогда для оценки абсолютной
погрешности получим выражение Ay = \f'(a)\Aa.
Аналогичное выражение можно записать для функции нескольких
аргументов. Например, оценка абсолютной погрешности функции и =
= f(x,y,z), приближенные значения аргументов которой соответствен-
соответственно а, 6, с, имеет вид
Au=\ti(x,y,z)\Aa+\f'(x,y,z)\Ab+\f'z(x,y,z)\Ac. A.4)

20	Гл. 1. Точность вычислительного эксперимента
Здесь Да, Д6, Ас — абсолютные погрешности аргументов. Относитель-
Относительная погрешность находится по формуле
Аи
ди =
|/(а,Ь,с)Г
Полученные соотношения можно использовать для вывода оценки пог-
погрешности произвольной функции (таким способом легко получить выра-
выражения A.3)). Например, при с = а — Ъ по формуле A.4) получа-
получаем Ас = \с'а\ Аа + \с'ь\АЪ = Аа + АЪ.
§ 2. Погрешности вычислений
1. Источники погрешностей. На некоторых этапах решения задачи на
компьютере могут возникать погрешности, искажающие результаты вычис-
вычислений. Оценка степени достоверности получаемых результатов является
важнейшим вопросом при организации вычислительных работ. Это осо-
особенно важно при отсутствии опытных или других данных для сравнения,
которое могло бы в некоторой степени показать надежность используемого
численного метода и достоверность получаемых результатов.
Рассмотрим источники погрешностей на отдельных этапах решения за-
задачи.
Математическая модель, принятая для описания данного процесса или
явления, может внести существенные погрешности, если в ней не учтены
какие-либо важные черты рассматриваемой задачи. В частности, математи-
математическая модель может прекрасно работать в одних условиях и быть со-
совершенно неприемлемой в других; поэтому важно правильно учитывать
область ее применимости.
Исходные данные задачи часто являются основным источником пог-
погрешностей. Вместе с погрешностями, вносимыми математической моде-
моделью, их называют неустранимыми погрешностями, поскольку они не могут
быть уменьшены вычислителем ни до начала решения задачи, ни в процес-
процессе ее решения. Проведенный ранее анализ оценки погрешностей при вы-
выполнении арифметических операций показывает, что следует стремиться к
тому, чтобы все исходные данные были примерно одинаковой точнос-
точности. Сильное уточнение одних исходных данных при наличии больших
погрешностей в других, как правило, не приводит к повышению точности
результатов.
Численный метод также является источником погрешностей. Это свя-
связано, например, с заменой интеграла суммой, с усечением рядов при вы-
вычислениях значений функций, с интерполированием табличных данных
и т. п. Как правило, погрешность численного метода регулируема, т. е.
теоретически она может быть уменьшена до любого значения путем из-
изменения некоторого параметра (например, шага интегрирования, числа
членов усеченного ряда и т. п.). Погрешность метода обычно стараются

§ 2. Погрешности вычислений	21
довести до величины, в несколько раз меньшей неустранимой погрешности.
Дальнейшее снижение погрешности не приведет к повышению точности
результатов, а лишь увеличит стоимость расчетов из-за необоснованного
увеличения объема вычислений. Подробнее погрешности методов будем
рассматривать при анализе конкретных численных методов.
При вычислениях с помощью компьютера неизбежны погрешности ок-
округлений, связанные с ограниченностью разрядной сетки машины. При
обычном округлении (которое, как правило, и реализуется в компьютерах)
максимальная относительная погрешность есть
<W = 0.5a1-*,	A.5)
где а — основание системы счисления, к — количество разрядов мантис-
мантиссы числа. При простом отбрасывании лишних разрядов эта погрешность
увеличивается вдвое.
Вычислим по формуле A.5) максимальную погрешность округления
#тах для чисел, представленных в форматах с одинарной и двойной точно-
точностью стандарта IEEE 754. Имеем: а = 2 в обоих случаях, для одинарной
точности к = 24 и 5max ~ 6 • 10~8, для двойной точности к = 53 и
<Wx и Ю-16.
Несмотря на то, что при решении больших задач выполняются милли-
миллиарды и триллионы операций, это вовсе не означает механического умно-
умножения погрешности при одном округлении на число операций, так как при
отдельных действиях погрешности могут компенсировать друг друга (на-
(например, при сложении чисел разных знаков). Вместе с тем иногда погреш-
погрешности округлений в сочетании с плохо организованным алгоритмом могут
сильно исказить результаты. В дальнейшем мы рассмотрим такие случаи.
Перевод чисел из одной системы счисления в другую также может
быть источником погрешности из-за того, что основание одной системы
счисления не является степенью основания другой (например, 10 и 2). Это
может привести к тому, что в новой системе счисления число становится
иррациональным.
Например, число 0.1 при переводе в двоичную систему счисления при-
примет вид 0.000 1100 1100 ... Может оказаться, что с шагом 0.1 нужно при
вычислениях пройти отрезок [0,1] от х = 1 до х = 0; десять шагов не
дадут точного значения х = 0.
2. Уменьшение погрешностей. При рассмотрении погрешностей
результатов арифметических операций отмечалось, что вычитание близ-
близких чисел приводит к увеличению относительной погрешности; поэтому в
алгоритмах следует избегать подобных ситуаций. Рассмотрим также неко-
некоторые другие случаи, когда можно избежать потери точности правильной
организацией вычислений.
Пусть требуется найти сумму пяти четырехразрядных чисел: S =
= 0.2764 + 0.3944 + 1.475 + 26.46 + 1364. Складывая все эти числа, а затем
округляя полученный результат до четырех значащих цифр, получаем S =
= 1393. Однако при вычислении на компьютере округление происходит

22	Гл. 1. Точность вычислительного эксперимента
после каждого сложения. Предполагая условно сетку четырехразрядной,
проследим за вычислением на компьютере суммы чисел от наименьше-
наименьшего к наибольшему, т. е. в порядке их записи: 0.2764 + 0.3944 = 0.6708,
0.6708 + 1.475 = 2.156, 2.156 + 26.46 = 28.62, 28.62 + 1364 = 1393;
получили Si = 1393, т. е. верный результат. Изменим теперь порядок
вычислений и начнем складывать числа последовательно от последнего к
первому: 1364+26.46 = 1390, 1390 + 1.475 = 1391, 1391+0.3944 = 1391,
1391 + 0.2764 = 1391; здесь окончательный результат S2 = 1391, он менее
точный.
Анализ процесса вычислений показывает, что потеря точности здесь
происходит из-за того, что прибавления к большому числу малых чисел не
происходит, поскольку они выходят за рамки разрядной сетки (а+6 = а при
а ^> Ъ). Этих малых чисел может быть очень много, но на результат они все
равно не повлияют, поскольку прибавляются по одному, что и имело место
при вычислении S2- Здесь необходимо придерживаться правила, в соответ-
соответствии с которым сложение чисел нужно проводить по мере их возрастания.
В машинной арифметике из-за погрешности округления существен порядок
выполнения операций, и известные из алгебры законы коммутативности
(и дистрибутивности) здесь не всегда выполняются.
При решении задачи на компьютере нужно использовать подобного рода
маленькие хитрости для улучшения алгоритма и снижения погрешностей
результатов. Например, при вычислении на компьютере значения (а + хJ
величина х может оказаться такой, что результатом сложения а + х полу-
получится а (при ж < а); в этом случае может помочь замена (а + хJ = а2 +
+ 2ах + х2 = а(а + 2х) + х2. Действительно, теперь к а прибавляется
не х, а 2х. Если же при х <С а вычисляется величина (а + хJ — а2, то це-
целесообразно привести ее к виду 2ах + х2, избежав тем самым вычитания
близких величин.
Рассмотрим еще один важный пример — использование рядов для
вычисления значений функций. Запишем, например, разложение функ-
функции sin х по степеням аргумента:
х3 хъ х7
По признаку Лейбница остаток сходящегося знакочередующегося ряда,
т. е. погрешность суммы конечного числа членов, не превышает значе-
значения первого из отброшенных членов (по абсолютной величине).
Вычислим значение функции sin ж при х = 0.5236 C0°). Члены ряда,
меньшие 10~4, не будем учитывать. Вычисления проведем с четырьмя
верными знаками. Получим
sin 0.5236 = 0.5236 - 0.2392 ¦ КГ1 + 0.3279 • ИГ3 = 0.500.
Это отличный результат в рамках принятой точности. Зная из курса высшей
математики, что это разложение синуса справедливо при любом значении
аргумента (—оо < х < +оо), используем его для вычисления функции

§ 2. Погрешности вычислений	23
при х = 6.807 C90°). Опуская вычисления, получаем sin 6.807 и 0.5167.
Относительная погрешность составляет здесь около 3 % (вместо ожида-
ожидаемого значения 0.01 % по признаку Лейбница). Это объясняется погреш-
погрешностями округлений и способом суммирования ряда (слева направо, без
учета величины членов).
Не всегда помогает и повышенная точность вычислений. В част-
частности, для данного ряда при х = 25.6563 ... A470° = 4 • 360° +30°) даже
при учете членов ряда до 10~8 и вычислениях с восемью значащими циф-
цифрами в результате аналогичных вычислений (суммирование слева направо)
получается результат, не имеющий смысла: sin ж и 129.
В программах, использующих степенные ряды для вычисления зна-
значений функций, могут быть приняты различные меры по предотвраще-
предотвращению подобной потери точности. Так, влияние погрешностей округления
существенно уменьшается, если \х\ < 1. Действительно, при вычисле-
вычислении хк допускается абсолютная погрешность
А(хк) = хкд(хк) = хккдх
(см. A.3)), которая при невыполнении неравенства \х\ < 1 может стать
неприемлемо большой. Для тригонометрических функций можно исполь-
использовать формулы приведения, благодаря чему аргумент будет находиться на
отрезке [0, 1]. При вычислении экспоненты аргумент х можно разбить
на сумму целой и дробной частей (ех = еп+а = еп • еа, 0 < а < 1) и
использовать разложение в ряд только для еа, а еп вычислять умно-
умножением. Таким образом, при организации вычислений можно своевре-
своевременно распознать «подводные камни», дающие потерю точности, и по-
попытаться затем исправить положение.
3. О решении квадратного уравнения. Мы убедились в том, что
при численном решении задач на компьютере вычислителя ожидают вся-
всякие «ловушки», которые могут привести к заметной потере точности ре-
результатов или даже к прекращению счета. Хорошей иллюстрацией к этому
является анализ алгоритма решения такой простой задачи, как решение
квадратного уравнения ах2 + Ьх + с = 0. Его корни определяются соот-
соотношениями
-b-VD	-b + VD	2
xi =	, х2 =	, D = b -4ac.	A.6)
la	la
Из анализа этих формул видно, что здесь имеется ряд особенностей вычис-
вычислительного характера, которые необходимо иметь в виду при составлении
алгоритма.
Рассмотрим простейший случай а = 0. Здесь уравнение становится
линейным, и его единственный корень есть х = -с/b, если Ь ф 0.
При а = Ъ = 0 и с ф 0 уравнение не имеет решения, а в случае а =
= Ь = с = 0 его решением будет любое число. Заметим, что в машинной
арифметике редко получаются точно нулевые значения. Поэтому коэффи-
коэффициенты можно сравнивать не с нулем, а с некоторой малой величиной е.

24	Гл. 1. Точность вычислительного эксперимента
Это в свою очередь порождает ряд ситуаций, зависящих от соотношения
между коэффициентами.
Далее необходимо предусмотреть разветвление алгоритма в зависимо-
зависимости от знака дискриминанта D: D > 0 — корни действительные (см. A.6));
D = 0 — корни равные: х\ = х2 = —Ь/Bа); D < 0 — корни комплекс-
комплексные: xlj2 =R±iI, где R = -Ъ/Bа), I = V^D/Ba).
Менее очевидным вопросом является возможность появления погреш-
погрешностей в зависимости от соотношения между коэффициентами уравнения.
Рассмотрим один из важнейших случаев, когда коэффициент Ъ значитель-
значительно превышает по абсолютной величине остальные. При этом Ъ2 ^> 4ас
и возникает опасность вычитания близких чисел в числителе одного из
выражений A.6) из-за того, что \[Ъ и |6|.
Положение можно исправить разными способами. Например, при Ъ > О
формулу для Х2 можно преобразовать следующим образом:
л/D-b VD + b	2с
la
При b < 0 аналогичным способом можно записать формулу для х\.
Более универсальным способом является использование значения sign Ь
(«знак величины Ъ »):
Тогда один из корней может быть вычислен по формуле
x^_b + Signb.VD
2а
Выражение для вычисления значения второго корня можно получить с по-
помощью теоремы Виета. Из соотношения х\Хч = с/а следует, что
х2 = —.	A.9)
ах\
На рис. 1.1 и 1.2 представлены структурограмма (см. приложение А) и
(для сравнения) блок-схема одного из вариантов алгоритма решения квад-
квадратного уравнения с учетом рассмотренных здесь особенностей. При D > О
значения корней вычисляются по формулам A.8), A.9). Заметим, что в при-
приведенном на структурограмме алгоритме предусмотрены еще не все случаи
возможных вычислительных затруднений, которые могут встретиться при
решении квадратных уравнений.
Можно привести некоторые примеры, когда реализация этого алгоритма
на компьютере невозможна. Будем предполагать, что вычисления проводят-
проводятся с двойной точностью.
Пример 1. а = Ю-200, Ъ = -3 • ИГ200, с = 2 • ИГ200. При вы-
вычислении произведений Ь2 и 4ас получается машинный нуль, т. е. D = 0;

§ 2. Погрешности вычислений
25
Да
Да
^Ч. С :
Да \^
Вывод:
х любое
Вывод:
решения
нет
Ввод <
а -
0 s'
с
Вывод х
2,6, С
= 0
^^
Да
Да^<
Вычисле-
Вычисление
Вывод
D = b2-
^^ D *
^^^Нет
Ж1-Ж2-
Ь
~ ~2а
XI, Х2
Нет
- 4ас
' ° /^
Вычисле-
Вычисление
RJ
Вывод
D,R,I
Рис. 1.1. Структурограмма алгоритма решения квадратного уравнения
^Начало J
х = — с/Ъ 	*у Вывод хг—>Г
Конец
ъ
Вычисление
Вывод /_
Рис. 1.2. Блок-схема алгоритма решения квадратного уравнения
решение пойдет по ветви равных корней: х\ = x<i = 1.5. Точные значения
корней, как нетрудно видеть, равны х\ = 1, Х2 = 2.

26	Гл. 1. Точность вычислительного эксперимента
Пример 2. а = 10200, Ь = -3 • 10200, с = 2 • ИГ200. Этот ва-
вариант аналогичен предыдущему случаю с той лишь разницей, что вместо
получения машинного нуля произойдет переполнение и прерывание счета.
Пример 3. а = ИГ200, Ъ = 10200, с = -10200. Это трудный для
реализации на компьютере случай. В практических расчетах встречаются
уравнения с малым коэффициентом при х2. В этом случае Ъ2 ^> 4ас, но
при вычислении Ъ2 произойдет переполнение. Простейшим выходом из
этого положения может быть сведение к случаю а = 0 с обязательной
проверкой других коэффициентов.
Таким образом, анализ даже такой задачи, как решение квадратного
уравнения, показывает, что использование численного алгоритма может
быть сопряжено с некоторыми трудностями.
§ 3. Устойчивость. Корректность. Сходимость
1. Устойчивость. Рассмотрим погрешности исходных данных. Пос-
Поскольку это так называемые неустранимые погрешности и вычислитель не
может с ними бороться, то нужно хотя бы иметь представление об их влия-
влиянии на точность окончательных результатов. Конечно, мы вправе надеяться
на то, что погрешность результатов имеет порядок погрешности исходных
данных. Всегда ли это так? К сожалению, нет. Некоторые задачи весьма
чувствительны к неточностям в исходных данных. Эта чувствительность
характеризуется так называемой устойчивостью.
Пусть в результате решения задачи по исходному значению величи-
величины х находится значение искомой величины у. Если исходная величина
имеет абсолютную погрешность Ах, то решение имеет погрешность Ау.
Задача называется устойчивой по исходному параметру х, если решение у
непрерывно от него зависит, т. е. малое приращение исходной величины Ах
приводит к малому приращению искомой величины Ау. Другими словами,
малые погрешности в исходной величине приводят к малым погрешностям
в решении.
Отсутствие устойчивости означает, что даже незначительные погреш-
погрешности в исходных данных приводят к большим погрешностям в решении
или даже к неверному результату. О неустойчивых задачах также говорят,
что они чувствительны к погрешностям исходных данных.
Приведем пример неустойчивой задачи. Рассмотрим квадратное урав-
уравнение с параметром a
х2 — 2х + sign a = 0.
Функция sign определена в A.7). Решение этого уравнения в зависимости
от значения а таково: х\ = x<i = 1 при а ^ 0; х\^ = 1 ± \[2 при a < 0.
Очевидно, что при a = 0 сколь угодно малая отрицательная погрешность
в задании а приведет к конечной, а не сколь угодно малой погрешности в
решении уравнения.

§ 3. Устойчивость. Корректность. Сходимость	27
Иногда бывает, что теоретически задача устойчива, но тем не менее
чувствительна к погрешностям исходным данных. Приращения исходной
величины, гарантирующие малость приращения искомой величины, оказы-
оказываются в этом случае настолько малыми, что реальные малые приращения,
с которыми имеет дело вычислитель, приводят к большим погрешностям
в решении.
Интересной иллюстрацией такой задачи является так называемый при-
пример Уилкинсона. Рассматривается многочлен
Р[х) = {х- 1)(х - 2)... [х - 20) = х20 - 210?19 + ...
Очевидно, что корнями этого многочлена являются х1 = 1, х2 = 2, ...
... , Х2о = 20. Предположим, что один из коэффициентов многочлена
вычислен с некоторой малой погрешностью. Например, коэффициент —210
при х19 увеличим на 10~7. В результате вычислений с двойной точ-
точностью получим существенно другие значения корней. Приведем для на-
наглядности эти значения, округленные до трех значащих цифр:
хх = 1.00,	х9 = 8.93,
х2 = 2.00,	жю,11 = 10.1 ± 0.60Н,
х3 = 3.00,	^12,13 = 11-8 ± 1.60г,
х4 = 4.00,	?14,15 = 14.0 ± 2.45г,
хъ = 5.00,	?16,17 = 16.7 ± 2.73г,
?6 = 6.00,	?18,19 = 19.5 ± 1.87г,
?7 = 7.00,	?20 = 20.8.
?8 = 8.01,
Таким образом, изменение коэффициента при х19 с -210 до -210 +
+ 10~7 (а это, несомненно, малое изменение в обычной вычислительной
практике) привело к тому, что половина корней стали комплексными. При-
Причина такого явления — чувствительность задачи к погрешностям исход-
исходным данных; вычисления выполнялись достаточно точно, и погрешности
округления не могли привести к таким последствиям. Заметим, что если ко-
коэффициент —210 изменить на значительно меньшее число, чем 10~7, то
изменение значений корней станет малым. Например, при увеличении коэф-
коэффициента —210 на 10~п значения корней, округленные до трех знаков,
совпадут со значениями корней исходного многочлена. Примерно такого
результата и следовало ожидать, поскольку рассматриваемая нами задача
устойчива.
2. Корректность. Задача называется поставленной корректно, если для
любых значений исходных данных из некоторого класса ее решение сущест-
существует, единственно и устойчиво по исходным данным.
Рассмотренная выше неустойчивая задача является некорректно постав-
поставленной. Применять для решения таких задач численные методы, как пра-
правило, нецелесообразно, поскольку возникающие в расчетах погрешности
округления будут сильно возрастать в ходе вычислений, что приведет к
значительному искажению результатов.

28	Гл. 1. Точность вычислительного эксперимента
Вместе с тем отметим, что в настоящее время развиты методы решения
некоторых некорректных задач. Это в основном так называемые методы
регуляризации. Они основываются на замене исходной задачи коррект-
корректно поставленной задачей. Последняя содержит некоторый параметр, при
стремлении которого к нулю решение этой задачи переходит в решение
исходной задачи.
3. Неустойчивость методов. Иногда при решении корректно пос-
поставленной задачи может оказаться неустойчивым метод ее решения. Такие
случаи имели место в § 2. В частности, по этой причине при вычислении
синуса большого аргумента был получен результат, не имеющий смысла.
Рассмотрим еще один пример неустойчивого алгоритма. Построим числен-
численный метод вычисления интеграла
Интегрируя по
1
г
/i = xex~1dx
0
i
In
частям,
= хе-
= хпех-Чх,
J
0
находим
1
л\1- \ex~1dx
0
1
71 = 1,2,...
1
е'
h = x2ex~1dx =
l - 2 xex~1dx = 1 - 2/ь
j
1	1
In = 1хпех-Чх = xnex-1\10 - n [x^e^dx = 1 - n/n_b
о	о
Пользуясь полученным рекуррентным соотношением, вычисляем с
двойной точностью (приводим результат, округленный до трех знача-
значащих цифр)
h = 0.368,	h = 0.146,
h = 0.264,	....
/3 = 0.207,	/it = 0.0558,
h = 0.171,	/i8 = -0.00369.
Значение интеграла /is не может быть отрицательным, посколь-
поскольку подынтегральная функция х18ех~1 на всем отрезке интегрирования
[0,1] неотрицательна. Исследуем источник погрешности. Максимальная
абсолютная погрешность при вычислении 1\ равна 0.5-2~53 и 5-10~17.
Однако на каждом этапе эта погрешность умножается на число, модуль
которого больше единицы (—2, —3,..., —18), что в итоге дает 18! и
и 6.4 • 1015. Это и приводит к результату, не имеющему смысла. Здесь

Упражнения	29
снова причиной накопления погрешностей является алгоритм решения
задачи, который оказался неустойчивым.
Численный алгоритм (метод) называется корректным в случае су-
существования и единственности численного решения при любых значениях
исходных данных, а также в случае устойчивости этого решения относи-
относительно погрешностей исходных данных.
4. Понятие сходимости. При анализе точности вычислительного
процесса одним из важнейших критериев является сходимость числен-
численного метода. Она означает близость получаемого численного решения
задачи к истинному решению. Строгие определения разных оценок бли-
близости могут быть даны лишь с привлечением аппарата функционального
анализа. Здесь мы ограничимся некоторыми понятиями сходимости, необ-
необходимыми для понимания последующего материала.
Рассмотрим понятие сходимости итерационного процесса. Этот про-
процесс состоит в том, что для решения некоторой задачи и нахождения ис-
искомого значения определяемого параметра (например, корня нелинейного
уравнения) строится метод последовательных приближений. В результате
многократного повторения этого процесса (или итераций) получаем по-
последовательность значений xi,X2,...,xn,... Говорят, что эта последова-
последовательность сходится к точному решению х = а, если при неограниченном
возрастании числа итераций предел этой последовательности существует и
равен a: lim xn = а. В этом случае имеем сходящийся численный метод.
п-юо
Другой подход к понятию сходимости используется в методах дис-
дискретизации. Эти методы заключаются в замене задачи с непрерывными
параметрами на задачу, в которой значения функций вычисляются в
фиксированных точках. Это относится, в частности, к численному ин-
интегрированию, решению дифференциальных уравнений и т. п. Здесь под
сходимостью метода понимается стремление значений решения дис-
дискретной модели задачи к соответствующим значениям решения ис-
исходной задачи при стремлении к нулю параметра дискретизации (на-
(например, шага интегрирования).
При рассмотрении сходимости важными понятиями являются ее
вид, порядок и другие характеристики. С общей точки зрения эти по-
понятия рассматривать здесь нецелесообразно; к ним будем обращаться
при изучении конкретных численных методов.
Таким образом, для получения решения задачи с необходимой точ-
точностью ее постановка должна быть корректной, а используемый численный
метод должен обладать устойчивостью (корректностью) и сходимостью.
Упражнения
1.	Представить числа 175.4, —3.169, —0.00874 в нормализованном виде.
2.	Записать в форме с фиксированной точкой числа 0.312 • 103, —0.70 • 101,
0.465-Ю-2.
Т. Получить соотношения A Л) и A.2).

30
Гл. 1. Точность вычислительного эксперимента
4.	Указать максимально возможные абсолютные и относительные погрешно-
погрешности приближенных чисел 27, -14.0, 0.00173, 0.745 • 10~4, -0.245 • 104, -
-0.8960 -102.
5.	Оценить погрешности величин х, у, заданных соотношениями
У =
при а и 32, Ъ и 17, с и 3.7.
6. Найти относительные погрешности при вычислении определителей
7.
8.
9.
0.19 -0.27
1.4 2.3
17.5 10.4
10.4 6.18
Каковы относительные погрешности объема шара и площади поверхности сфе-
сферы, если их радиус известен с точностью до 10 %?
Убедиться, что результат вычисления суммы чисел 103, 0.704, 0.537, 15.2
с округлением до трех значащих цифр после каждого сложения зависит от
порядка суммирования. В каком порядке следует суммировать эти числа?
Что произойдет, если с помощью алгоритма, представленного на рис. 1.1,
попытаться решить квадратное уравнение с коэффициентами а = 10180, Ь =
= —2 • 10180, с = —10~180? Как можно преодолеть возникшие трудности?

ГЛАВА 2
АППРОКСИМАЦИЯ ФУНКЦИЙ
§ 1. Понятие о приближении функций
1. Постановка задачи. Пусть величина у является функцией аргу-
аргумента х. Это означает, что любому значению х из области определения
поставлено в соответствие значение у. Вместе с тем на практике часто
неизвестна явная связь между у и х, т. е. невозможно записать эту связь в
виде некоторой зависимости у = f(x) .Иногда даже известная зависимость
у = f(x) оказывается настолько громоздкой (например, содержит трудно
вычисляемые выражения, сложные интегралы и т. п.), что ее использование
в практических расчетах требует слишком много времени.
Наиболее распространенным и практически важным случаем, когда вид
связи между параметрами х и у неизвестен, является задание этой связи
в виде некоторой таблицы {xi,yi}. Это означает, что дискретному мно-
множеству значений аргумента {xi} поставлено в соответствие множество
значений функции {yi} (i = 0,1,..., п). Эти значения — либо результаты
расчетов, либо экспериментальные данные. На практике нам могут пона-
понадобиться значения величины у и в других точках, отличных от узлов Х{.
Однако получить эти значения можно лишь путем очень сложных расчетов
или проведением дорогостоящих экспериментов.
Таким образом, с точки зрения экономии времени и средств мы при-
приходим к необходимости использования имеющихся табличных данных для
приближенного вычисления искомого параметра у при любом значении (из
некоторой области) определяющего параметра х, поскольку точная связь
у = f(x) не известна (либо нам нецелесообразно ею пользоваться).
Этой цели и служит задача о приближении (аппроксимации) функций:
данную функцию f(x) требуется приближенно заменить (аппроксими-
(аппроксимировать) некоторой функцией ip(x) так, чтобы отклонение (в некотором
смысле) ip(x) от f(x) в заданной области было наименьшим. Функ-
Функция (р(х) при этом называется аппроксимирующей.
Аппроксимация рассмотренного выше типа, при которой приближение
строится на заданном дискретном множестве точек {х^} , называется то-
точечной. К ней относятся интерполирование, среднеквадратичное прибли-
приближение и др. При построении приближения на непрерывном множестве то-
точек (например, на отрезке) аппроксимация называется непрерывной (или
интегральной). К непрерывной аппроксимации относится, например, рав-
равномерное приближение.
Для практики весьма важен случай аппроксимации функции много-
многочленом
ср(х) = Рт(х) = а0 + а\х + а2х2 + ... + ашхш.	B.1)

32	Гл. 2. Аппроксимация функций
В дальнейшем будем чаще всего рассматривать аппроксимацию такого ро-
рода. При этом коэффициенты ctj будут подбираться так, чтобы достичь наи-
наименьшего отклонения многочлена от данной функции.
Что касается самого понятия «малое отклонение», то оно будет уточне-
уточнено в дальнейшем — при рассмотрении конкретных способов аппроксима-
аппроксимации.
2. Точечная аппроксимация. Одним из основных типов точечной ап-
аппроксимации является интерполирование. Оно состоит в следующем: для
данной функции у = f(x) строим интерполирующую функцию <р(х) (на-
(например, многочлен B.1)), принимающую в заданных точках Х{, те же зна-
значения yi, что и функция f(x), т. е.
(f(xi)=yi, i = 0,1,..., п.	B.2)
При этом предполагается, что среди значений xi нет одинаковых, т. е.
xi ф Xk при г ф к. Точки xi называются узлами интерполяции.
Таким образом, близость интерполирующей функции (рис. 2.1, сплош-
сплошная линия) к заданной функции состоит в том, что их значения совпадают
на заданной системе точек.
Интерполирующая функция (р(х) может строиться сразу для всего рас-
рассматриваемого интервала изменения х или отдельно для разных частей
этого интервала. В первом случае гово-
говорят о глобальной интерполяции, во втором
— о кусочной (или локальной) интерполя-
интерполяции.
Как правило, интерполирование ис-
используется для аппроксимации функции в
промежуточных точках между крайними
узлами интерполяции, т. е. при ж о < х <
< хп. Однако иногда оно применяется
и для приближенного вычисления функ-
хп х ции вне рассматриваемого отрезка (х <
Рис. 2.1. Интерполяция и аппрок- < х* > х > х") • Это приближение назы-
симаиия	вают экстраполяцией.
Рассмотрим использование в качестве
функции ip(x) многочлена B.1), называе-
называемого интерполяционным многочленом. При глобальной интерполяции, т. е.
при построении одного многочлена для всего рассматриваемого интервала
изменения х, для нахождения коэффициентов многочлена необходимо ис-
использовать все уравнения системы B.2). Данная система содержит п + 1
уравнение, следовательно, с ее помощью можно определить п + 1 коэф-
коэффициент. Поэтому максимальная степень интерполяционного многочлена
га = п, и многочлен принимает вид
... + апхп.	B.3)

§ 1. Понятие о приближении функций	33
Система уравнений B.2) при использовании в качестве (р(х) многочле-
многочлена B.3) является системой линейных алгебраических уравнений относи-
относительно неизвестных коэффициентов а$ , а\, ... , ап и имеет вид
а0 + агхг + ... + апж? = Уъ
а0 + ai^n + ... + апж™ = уп.
Определитель такой системы в линейной алгебре называется определите-
определителем Вандермонда. Можно показать, что определитель Вандермонда отличен
от нуля, если Х{ ф Xk при г ф к, т. е. если среди узлов интерполяции
нет совпадающих. Следовательно, в этом случае система B.4) имеет един-
единственное решение. Решив систему B.4), можно построить интерполяцион-
интерполяционный многочлен. Такой метод построения интерполяционного многочлена
называется методом неопределенных коэффициентов. Заметим вместе с
тем, что этот метод требует значительного объема вычислений, особенно
при большом числе узлов. Существуют более простые алгоритмы построе-
построения интерполяционных многочленов, которые будут рассмотрены в § 3.
Как видим, при интерполировании основным условием является про-
прохождение графика интерполирующей функции через данные значения
функции в узлах интерполяции. Однако в ряде случаев выполнение этого
условия затруднительно или даже нецелесообразно.
Например, при большом количестве узлов интерполяции в случае
глобальной интерполяции получается высокая степень многочлена B.3).
Кроме того, табличные данные могли быть получены путем измерений и
содержать ошибки. Построение аппроксимирующей функции с услови-
условием обязательного прохождения ее графика через эти экспериментальные
точки означало бы тщательное повторение допущенных при измерениях
ошибок. Выход из этого положения может быть найден выбором такой
функции, график которой проходит близко от данных точек (см. рис. 2.1,
штриховая линия). Понятие «близко» уточняется при рассмотрении раз-
разных видов приближения.
Одним из таких видов является среднеквадратичное приближение.
Если при этом используется многочлен B.1), то га ^ п; случай га = п
соответствует глобальной интерполяции. На практике стараются подо-
подобрать аппроксимирующую функцию как можно более простого вида,
например, многочлен степени га = 1, 2, 3.
Мерой отклонения функции (р(х) от заданной функции f(x) на мно-
множестве точек (xi^yi) (i = 0,1,...,п) при среднеквадратичном прибли-
приближении является величина S, равная сумме квадратов разностей между
значениями аппроксимирующей и заданной функции в данных точках:
г=0
3 Л.И. Турчак, П.В. Плотников

34
Гл. 2. Аппроксимация функций
Аппроксимирующую функцию нужно подобрать так, чтобы величи-
величина S была наименьшей. В этом состоит метод наименьших квадратов,
который будет рассмотрен в § 4.
3. Непрерывная аппроксимация. Равномерное приближение. Во
многих случаях, особенно при обработке экспериментальных данных,
среднеквадратичное приближение вполне приемлемо, поскольку оно сгла-
сглаживает некоторые неточности функции f(x) и дает достаточно правильное
представление о ней. Иногда, однако, при построении приближения ста-
ставится более жесткое условие: требуется, чтобы во всех точках некоторого
отрезка [а, Ь] отклонение аппроксимирующей функции (р(х) от функции
f(x) было по абсолютной величине меньшим заданной величины е > 0:
\f(x) - <р(х)\ < ?, а ^х ^Ъ.
В этом случае говорят, что функция ip(x) равномерно приближает {ап-
{аппроксимирует) функцию f(x) с точностью е на отрезке [а, Ь\.
Понятие равномерного приближения предполагает сравнение задан-
заданной и аппроксимирующей функций на непрерывном множестве — отрез-
отрезке [а, Ъ]. Поэтому равномерное приближение относится к непрерывной
аппроксимации.
Введем понятие абсолютного отклонения А функции (р(х) от функции
f(x) на отрезке [а, Ъ]. Оно равно максимальному значению абсолютной
величины разности между ними на данном отрезке:
А = max \f(x) — ip(x)
a<x<b
B.5)
По аналогии можно ввести понятие среднеквадратичного отклонения
А = л/Sjn при среднеквадратичном приближении функций. На рис. 2.2
показано принципиальное различие двух рассматриваемых приближений.
У = fix)
= <р{х)
у=Л
Рис. 2.2. Приближения: среднеквадратичное (а); равномерное (б)

§ 1. Понятие о приближении функций	35
Возможность построения многочлена, равномерно приближающего
данную функцию, следует из теоремы Вейерштрасса об аппрокси-
аппроксимации.
Теорема. Если функция f(x) непрерывна на отрезке [а, Ь\
то для любого е > 0 существует многочлен Рш(х) степени га = т(е),
абсолютное отклонение которого от функции f(x) на отрезке [а, Ъ] мень-
меньше е.
В частности, если функция f(x) на отрезке [а, Ь\ разлагается в рав-
равномерно сходящийся степенной ряд, то в качестве аппроксимирующего
многочлена можно взять частичную сумму этого ряда. Такой подход ши-
широко используется, например, при вычислении на компьютере значений
элементарных функций.
Существует также понятие наилучшего приближения функции f(x)
многочленом Рш(х) фиксированной степени га. В этом случае коэффи-
коэффициенты многочлена B.1) следует выбрать так, чтобы на заданном отрез-
отрезке [а, Ь\ величина абсолютного отклонения B.5) была минимальной. Та-
Такой многочлен Рт(х) называется многочленом наилучшего равномерного
приближения. Существование и единственность многочлена наилучшего
равномерного приближения вытекает из следующей теоремы.
Теорем а. Для любой функции f(x), непрерывной на замкнутом огра-
ограниченном множестве G, и любого целого га ^ 0 существует много-
многочлен Рш(х) степени не выше т, абсолютное отклонение которого от
функции f(x) среди всех многочленов степени не выше га минимально,
т. е. А = Amiib причем такой многочлен единственный.
Множество G обычно представляет собой некоторый отрезок [а, Ь].
4. Вычисление многочленов. При аппроксимации функций, а также в
некоторых других задачах приходится вычислять значения много-
многочленов вида
Рп(х) = «о + a>ix + а2х2 + ... + апхп.	B.6)
Если проводить вычисления «в лоб», т. е. находить значения каждого
члена и суммировать их, то при больших п потребуется выполнить боль-
большое число операций (п2 + п/2 умножений и п сложений). Кроме того,
это может привести к потере точности за счет погрешностей округления. В
некоторых частных случаях удается выразить каждый последующий член
через предыдущий и таким образом значительно сократить объем вычис-
вычислений.
Анализ многочлена B.6) в общем случае приводит к тому, что для исклю-
исключения возведения х в степень в каждом члене многочлен целесообразно
переписать в виде
Рп(х) = а0 + х(а\ + х(а2 + ... + x(an-i + хап)...).
Прием, с помощью которого многочлен представляется в таком виде,
называется схемой Горнера. Соответствующий ему алгоритм вычисления
з*

36
Гл. 2. Аппроксимация функций
значения многочлена изображен на рис. 2.3. Этот метод требует п умноже-
умножений и п сложений. Использование схемы Горнера для вычисления значений
многочленов не только экономит машинное время, но и повышает точность
вычислений за счет уменьшения погрешностей округления.
§ 2. Использование рядов
1. Элементарные функции. Как правило, при решении различных за-
задач приходится вычислять значения элементарных функций (тригономет-
(тригонометрических, показательных, логарифмических и др.). При ручном счете для
этой цели могут быть использованы таблицы. Однако в вычислениях на ком-
компьютере ввод таблиц функций в машину потребовал бы больших затрат па-
памяти. Кроме того, поиск нужного значения функции в памяти
компьютера не простое для машины занятие.
Поэтому для вычисления значений функций
на компьютере используются разложения этих
функций в степенные ряды. Например, функ-
функция sin х вычисляется с помощью ряда
для г
до 0
Ввод п,
Р = с
от п — 1
Р = ui
с шагом -
Вывод
iCLzj , X
+ хР
-1
р
_ + _
х7
_
B.7)
Рис. 2.3. Схема Горнера
При известном значении аргумента х зна-
значение функции может быть получено с
точностью до погрешностей округления. Коли-
Количество используемых членов ряда B.7) зави-
зависит от значения аргумента. Напомним, что в
соответствии с правилами приближенных вы-
вычислений (см. гл. 1) для предотвращения влияния погрешностей округления
необходимо выполнение неравенства \х\ < 1.
С помощью степенных рядов вычисляются значения и других элемен-
элементарных функций. В частности, для вычисления значений функции cos ж
можно использовать ряд B.7) с учетом соотношения cos х = sin (тг/2 + х),
а можно и непосредственно воспользоваться разложением в ряд функ-
функции cos х:
х2
х4
х6
B.8)
Гиперболические синус и косинус можно вычислить с помощью
разложения в ряд экспоненты ех , поскольку
sh х = (ех - е~х)/2, chx = (ех + е~х)/2.
Можно также воспользоваться разложением в ряд самих функций sh x
и ch х. Так, при вычислении sh x для х и 0 в целях предотвращения по-
потери точности из-за вычитания близких величин полезно использовать ряд

§ 2. Использование рядов	Ъ1
ж3 ж5 х7
shx = x+- + - + - + ...
Для вычисления на компьютере логарифмических функций достаточно
иметь программу вычисления логарифма по одному основанию, например
натурального логарифма. Для вычисления логарифма по другому основа-
основанию можно воспользоваться соотношением loga х = In х/ In a.
В качестве примера построим алгоритм вычисления синуса с по-
помощью ряда B.7). Будем учитывать члены ряда, которые по абсолют-
абсолютной величине больше некоторого малого числа е > 0, характеризую-
характеризующего точность вычисления. На практике, когда используют стандартные
программы для вычисления функций, точность не задается. В этом слу-
случае учитываются все члены, большие машинного нуля, а точность резуль-
результата определяется погрешностями округлений.
Возможный вариант алгоритма вычисления синуса с помощью ря-
ряда B.7) изображен на рис 2.4 в виде структурограммы. Дадим некоторые
пояснения к ней.
В алгоритме предусматривается выход из программы при малом значе-
значении аргумента, поскольку в этом случае sin x и х. Выделяется абсолютная
величина аргумента с учетом соотношений
х = к\х\, smx = ksm\x\, k = signx.
Далее полагается, что х > 0.
При анализе точности вычислений отмечалось, что при суммировании
ряда погрешность значительно меньше, если \х\ < 1. Поэтому в структуро-
грамме аргумент должен удовлетворять неравенству х < тг/4. Это дости-
достигается последовательным уменьшением аргумента до значений х < 2тг,
х < тг, х < тг/2. Для этой цели использована функция Е(х), вычисляющая
целую часть аргумента, а также формулы приведения sin (тг ± х) = =р
=F sin x, sin (тг/2 — х) = cos х. Например, при х = 7.6 тг (к = 1) получим
следующий алгоритм:
х - 2тгп = 7.6 тг — бтг = 1.6 тг > тг, к = -1,
тг тг тг
х — тг = О.бтг > —, тг — х = 0.4тг > —,	ж = 0.1тг.
А	т-	Zi
Текущее значение члена ряда на рис. 2.4 обозначено через и, значе-
значение функции — через у.
Здесь первые два члена ряда вычисляются непосредственно, а каж-
каждый последующий выражается через предшествующий. Например,
—, U3 = — = -U2——.

38
Гл. 2. Аппроксимация функций
"
Да
Да
Да
Да
Да
и =
—¦—^
к
\
^^
у = >.
-- -х3/6
пока \и ^
У =
——
¦—
= sign х, х
х <
Т <^
X <
г,
, m =
г
= У +
и
2тг
^^
п
7Г
х = х
тг/2
х-
^^ ^<
= 2
1А, 771
и
т(п
У =
ВВОД Ж, ?
"	.
= X
= х — 2тгп
— тг, к = —к
^^
— 7Г — X
_^' >^^^^ Нет
Нет
Нет
Нет
^^
Нет
	 /Q
Ж — 7Г / ^ — X ;
1А = — X2 /2 , 171=1
= ш + 2
2
, + 1)
Вывод 2/
Рис. 2.4. Алгоритм вычисления синуса

§2. Использование рядов
39
Если тг/4 < х < тг/2, то проводится уменьшение аргумента до величи-
величины тг/2 — х и вычисление синуса сводится к вычислению косинуса, т. е.
используется ряд B.8). Для этого полагаем первоначально у = 1, и = 1,
т = 1. В остальном алгоритм не меняется.
В рассмотренном алгоритме аргумент предполагается заданным в ра-
радианах. Если он задан в градусах, то следует предусмотреть перевод его в
радианы, т. е. умножение на величину тг/180.
Погрешность функции у = sin x, полученной с помощью ряда B.7) с ис-
использованием приведенного на рис. 2.4 алгоритма, состоит из двух частей —
погрешности округления и погрешности ограничения, возникающей из-за
учета лишь ограниченного числа членов ряда.
Погрешности ограничений при фиксированном числе членов ряда зави-
зависят от значения аргумента. При \х\ < тг/4 они весьма малы и сравнимы
с погрешностями округлений даже при небольшом числе членов ряда, а
с увеличением х возрастают. В частности, если ограничиться первыми
четырьмя членами разложения B.7) и провести вычисления с одинарной
точностью, то погрешность при х = тг/4 составит около 3 • 10~7, а при
х = тг/2 —уже около 1.6 • 10~4 (порядка первого отброшенного члена —
в соответствии с признаком Лейбница).
2. Многочлены Чебышева. Из приведенного выше примера вычисле-
вычисления синуса с помощью ряда следует, что погрешности могут быть распреде-
распределены неравномерно по рассматриваемому
аргумента. Одним из способов совершенство-
совершенствования алгоритма вычислений, позволяющих
более равномерно распределить погрешность
по всему интервалу, является использование
многочленов Чебышева.
Многочлен Чебышева Тп(х) степени п оп-
определяется следующей формулой:
Тп(х) = \\
интервалу изменения
= Т0(х)
Рис. 2.5. Многочлены Чебы-
Чебышева
— 1 < т < 1 п — О 1	(Лл)
Легко показать, что B.9) на самом деле яв-
является многочленом. Действительно, при воз-
возведении двучленов в степень п выражения
л/1 — х2 будут возведены в четные и нечетные степени. Эти выраже-
выражения в четных степенях станут рациональными, а в каждой из нечетных
степеней они будут присутствовать дважды с коэффициентами противо-
противоположных знаков, вследствие чего взаимно уничтожатся.
Приведем многочлены Чебышева, полученные по формуле B.9)
при п = 0,1,2,3 (рис. 2.5):
Г0(ж) = 1, Тг(х)=х, Т2(х) = 2х2-1, Т3(х) =4х3 -Зх.

40
Гл. 2. Аппроксимация функций
Для вычисления многочленов Чебышева можно воспользоваться сле-
следующим рекуррентным соотношением:
Тп+1(х) = 2хТп(х) - Тп_1(ж), п = 1, 2,...
B.10)
В ряде случаев важно знать коэффициент ап при старшем члене мно-
многочлена Чебышева степени п
Тп(х) = а0 ¦
Разделив этот многочлен на ж" найдем
Тп{х) а0
апхп
fln-i
х
B.11)
Определим многочлен Чебышева при \х\ > 1 аналогично B.9) с за-
заменой 1 - х2 на х2 - 1 и воспользуемся этим определением, переходя к
пределу при х —> оо в B.11). Получим
ап = lim
ж—>-ос
Тп(х)
= 2
п-1
B.12)
Многочлены Чебышева можно также представить в тригонометри-
тригонометрической форме:
Тп(х) = cos (narccosx), n = 0,1,...
С помощью этих выражений могут быть получены формулы B.9), B.10).
Нули {корни) многочленов Чебышева на отрезке [—1,1] определяются
формулой
2А;-1
zk = cos
2n
-тг, k = 1, 2,... ,n.
Они расположены неравномерно на отрезке и сгущаются к его концам.
Вычисляя экстремумы многочлена Чебышева по обычным правилам (с
помощью производных), можно найти его максимумы и минимумы:
хк = cos (&тг/п), к = 1, 2,..., п — 1.
В этих точках многочлен принимает поочередно значения Тп(хк) = =Ы,
т, е, все максимумы равны 1, а минимумы равны — 1. На границах отрезка
значения многочленов Чебышева равны ±1.
Многочлены Чебышева широко используются при аппроксимации
функций. Рассмотрим их применение для улучшения приближения функ-
функций с помощью степенных рядов, а именно для более равномерного
распределения погрешностей аппроксимации B.7) по заданному отрезку
[-7г/2,тг/2].
Отрезок [—тг/2, тг/2] является не совсем удобным при использова-
использовании многочленов Чебышева, поскольку они обычно рассматриваются на

§2. Использование рядов
41
стандартном отрезке [—1,1]. Первый отрезок легко привести ко второму
заменой переменной х на тгж/2. В этом случае ряд B.7) для аппрокси-
аппроксимации синуса на отрезке [—1,1] примет вид
7ГЖ
BЛЗ)
При использовании этого ряда погрешность вычисления функции в
окрестности концов отрезка х = ±1 су-
существенно возрастает и становится зна-
значительно больше, чем в окрестности точ-
точки х = 0. Если вместо B.13) использо-
использовать ряд
. 7ГЖ
sin — = со
¦ с^х) + с2Т2(х)
B.14)
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0 0.2 0.4 0.6
1 х
Рис. 2.6. Погрешность аппрокси-
аппроксимации функции у = sin (тгж/2)
с помощью рядов Тейлора и Че-
бышева с учетом членов до ж3
включительно
членами которого являются многочлены
Чебышева, то погрешность будет рас-
распределена равномерно по всему отрезку
(рис. 2.6). В частности, при использо-
использовании многочленов Чебышева до седь-
седьмой степени включительно (т. е. при ис-
использовании четырех ненулевых членов
B.14)) погрешность находится в интер-
интервале (—5.9 -г- 5.9) • 10 ~7. Для сравнения напомним, что при использова-
использовании четырех членов ряда Тейлора погрешность для этой задачи на концах
отрезка составляет ±1.6 • 10~4.
Для нахождения коэффициентов ряда Чебышева нужно воспользо-
воспользоваться свойством ортогональности многочленов Чебышева на отрез-
отрезке [—1,1] свесом 1/л/1 — х2:
-1
0, т Ф п.
m = п = 0.
B.15)
B.16)
умножить обе части равенства на Тп(х)/л/1 — х2, проинтегрировать ряд
почленно по отрезку [—1,1] и воспользоваться B.15), то получим
Если записать ряд Чебышева для функции f(x) в виде
/(«) = J + ciT^x) + с2Т2(х) + ... ,
п = 0,1,2,...
B.17)

42	Гл. 2. Аппроксимация функций
Заметим, что коэффициенты ряда Чебышева находятся аналогично
тому, как находятся коэффициенты известного из курса высшей матема-
математики тригонометрического ряда Фурье.
Входящие в B.17) интегралы в общем случае можно найти только
численно. На практике часто находят коэффициенты ряда Чебышева при-
приближенно, заменяя функцию f(x) частичной суммой ее ряда Тейлора,
содержащей достаточное для требуемой точности количество членов ря-
ряда. Рассмотрим этот подход более подробно.
Пусть частичная сумма ряда Тейлора, представленная в виде много-
многочлена, используется для приближения функции f(x) на стандартном от-
отрезке [—1,1], т. е.
f(x) и (р(х) = а0 + а\х + ... + апхп.	B.18)
Если рассматриваемый отрезок [а, Ъ] отличается от стандартного, то его
всегда можно привести к стандартному заменой переменной
а + Ъ Ь-а
t = —	1	— ж, -1 ^ х ^ 1.
Каждая степень хш может быть выражена через многочлены Чебы-
Чебышева степеней, не превышающих т (в приложении Б приведены соот-
соответствующие формулы). В результате ряд Чебышева для функции (р(х)
будет являться конечной суммой:
Ф) = f + С1ТЛХ) + • • • + «U-iTn-iM + c*nTn(x),	B.19)
где Cq , с\, ... , с^ — приближенные значения коэффициентов со,
с\, ... , сп в B.16). Если в такой сумме привести подобные слагаемые,
содержащие х в одинаковых степенях, т. е. явно представить ее в виде
многочлена, то этот многочлен совпадет с многочленом B.18).
Поскольку, как нам известно, погрешность при использовании ря-
ряда Чебышева меньше, чем при использовании ряда Тейлора, требуемую
точность можно сохранить и при меньшем количестве слагаемых в B.19).
Несколько последних из них можно отбросить.
Рассмотрим, что происходит при отбрасывании слагаемого с Тп . Мно-
Многочлен Тп(х) входит в частичной сумме B.18) только в выражение для
хп. Поэтому чтобы найти многочлен, получающийся после отбрасыва-
отбрасывания слагаемого с Тп в B.19), не нужно вычислять коэффициенты Cq ,
с\, ... , с^1_1. Достаточно выразить хп через Тп(х) и меньшие степе-
степени х, подставить полученное выражение в B.18), а затем отбросить
в B.18) слагаемое с Тп .
Оценим допускаемую при этом погрешность. Из B.12) следует, что
многочлен Чебышева Тп (х) можно записать в виде
Тп(х) = Ъо + Ьгх + Ъ2х2 + ... + 2n~1xn.

§ 2. Использование рядов	43
Отсюда получаем
хп = _21-«(ь0 + ъ1Х + ... + бп-i^) + 21~пТп{х).	B.20)
Если отбросить последний член, то допущенную при этом погреш-
погрешность Ai легко оценить: |Ai| ^ 21~п, поскольку |Гп(ж)| ^ 1. Тогда
погрешность А, допущенная в B.18) и B.19), будет оцениваться как
|А| ^ К^1"".
Из B.20) получаем, что хп есть с точностью до Ai линейная комби-
комбинация более низких степеней х. Подставляя эту линейную комбинацию
в B.18), приходим к многочлену степени п — 1 вместо многолена степе-
степени п. Этот процесс может быть продолжен до тех пор, пока погрешность
не превышает допустимого значения.
Описанную процедуру можно трактовать как повышение точности
аппроксимации функции с помощью ряда Тейлора. Используем эту про-
процедуру для повышения точности аппроксимации B.13). Поставим задачу
примерно сохранить погрешность 5.9 • 10~7, имевшую место при исполь-
использовании многочленов Чебышева до седьмой степени. Будем учитывать
члены ряда B.13) до 11-й степени включительно. При этом допускается
погрешность примерно б • 10~8. Вычисляя коэффициенты при степе-
степенях х, получаем, округляя до восьми разрядов,
sin (тгж/2) и 1.5707963 х - 0.64596410 х3 + 0.079692626 ж5 -
- 0.0046817541 ж7 + 0.00016044118 ж9 + 0.0000035988432 ж11. B.21)
Многочлен Чебышева 11-й степени имеет вид
Гц = 1024 ж11 - 2816 ж9 + 2816 ж7 - 1232 ж5 + 220 ж3 - Их.
Выразим отсюда х11 через более низкие степени:
х11 = 2-10A1ж - 220ж3 + 1232ж5 - 2816ж7 + 2816х9 + Гц).
Подставляя в B.21) вместо х11 правую часть этого равенства и вы-
вычисляя новые значения коэффициентов, получаем
sin (тгж/2) и 1.5707963 х - 0.64596332 ж3 + 0.079688296 ж5 -
- 0.0046718573ж7 + 0.00015054437ж9 - 0.00000000351 Гц. B.22)
Отбрасывая последний член этого разложения, мы допускаем погреш-
погрешность |А| ^ 3.51 • 10~9. Из-за приближенного вычисления коэффици-
коэффициентов при степенях х и погрешности ограничения реальная погреш-
погрешность больше. Она составляет 4-10~8 , что все равно значительно меньше
чем 5.9 • 10~7 . Поэтому процедуру отбрасывания слагаемого можно по-
повторить.
Подставляя в B.22) (без последнего члена с Гц ) выражение х9 через
более низкие степени

44	Гл. 2. Аппроксимация функций
х9 = 2~8(-9х + 120 ж3 - 432 ж5 + 576 х7 + Г9),
получаем
sin (тгж/2) и 1.5707910 ж - 0.64589276 ж3 + 0.079434253 ж5 -
- 0.0043331325ж7 + 0.000000588 Т9. B.23)
Отбрасывая последний член этого разложения, мы допускаем погреш-
погрешность |А| ^ 5.88 • 10~7.
Суммарная погрешность B.23) без последнего члена с Тд соста-
составит 6.4 • 10~7, что лишь немногим больше погрешности 5.9 • 10~7, до-
допускаемой при вычислении коэффициентов ряда Чебышева по B.17).
В приложении Б приведены некоторые формулы, необходимые при
использовании многочленов Чебышева.
3. Рациональные приближения. Рассмотрим другой вид аппрокси-
аппроксимации функций — с помощью дробно-рационального выражения. Функ-
Функцию представим в виде отношения двух многочленов некоторой степени.
Пусть, например, это будут многочлены третьей степени, т. е. предста-
представим функцию f(x) в виде дробно-рационального выражения
_ 60 + hx + Ъ2х2 + Ъ3х3
ПХ)	B24)
Значение свободного члена в знаменателе со = 1 не нарушает общности
этого выражения, поскольку при с0 ф 1 числитель и знаменатель мож-
можно разделить на со. Если же со = 0, то f(x) будет иметь особенность
при х = 0. Такую аппроксимацию рассматривать не будем.
Перепишем выражение B.24) в виде
&о + Ьгх + Ъ2х2 + Ъ3х3 = A + сгх + с2х2 + c3x3)f(x).
Используя разложение функции f(x) в ряд Тейлора
f(x) = а0 + агх + а2х2 + ...	B.25)
и учитывая члены до шестой степени включительно, получаем
Ъо + Ъгх + Ь2х2 + Ъ3х3 = A + сгх + с2х2 + с3х3) х
х (ао + а\х + а2х2 + а3х3	А	ъ	6)
Преобразуем правую часть, этого равенства, записав ее разложение по
степеням х:
Ьо + bix + Ъ2х2 + Ъ3х3 = а0 + x(ai + aoci) + х2(а2 + aiCi + а0с2) +
3	а0с3) + х4(а4 + a3ci + а2с2 + aic3) +
а3с2 + а2с3) + ж6(а6 + a5ci + а4с2 + а3с3).

§ 2. Использование рядов	45
Приравнивая коэффициенты при одинаковых степенях х в левой и
правой частях, получаем следующую систему уравнений:
bi = ai -\-
62 = &2 + aiCi +
Ьз = аз + a2ci + aic2 + а0с3,	B.26)
О = CL4 + CI3C1 + C12C2 +
О = as + CI4C1 + азС2 +
О = а6 + a5ci + CL4C2 + а3с3.
Решив эту систему, найдем коэффициенты 60, Ь1, Ъ2, 63, ci, с2, с3
необходимые для аппроксимации B.24).
Пример. Рассмотрим рациональное приближение для функции
f(x) = sin (тгж/2). Воспользуемся представлением B.24), которое в дан-
данном случае упрощается, поскольку функция sin ж нечетная. В частнос-
частности, в числителе можем оставить только члены с нечетными степенями
х, а в знаменателе — с четными; коэффициенты при других степенях х
равны нулю: Ьо = 62 = ci = с3 = 0.
Коэффициенты &i, 63 , С2 найдем из системы уравнений B.26), при-
причем значения коэффициентов а$ , ai,... , аб разложения функции в ряд
Тейлора B.25) можем взять из выражения B.13), т. е.
7Г	7Г3	7Г5
ао = ^2 = CL4 = ^6 = 0, а
2' 3 8-3!' 5 32-5!"
Система уравнений B.26) в данном случае примет вид
7Г
8-
7Г5
3
3!
7Г
|_
1 2
7Г
С2,
3
Отсюда находим Ъ\ = тг/2, 63 = -7тг3/480, с2 = тг2/80.
Таким образом, дробно-рациональное приближение B.24) для функ-
функции sin (тгж/2) примет вид
тгх _ Gг/2)Я - Gтг3/480)х3
8 2 " 1 + (тг2/80х2) '	Ц'27)
Это приближение по точности равносильно аппроксимации B.13) с уче-
учетом членов до пятого порядка включительно.
На практике с целью экономии числа операций выражение B.24) пред-
представляется в виде цепной дроби. Представим в таком виде дробно-раци-
дробно-рациональное выражение B.27). Сначала перепишем это выражение, вынося

46	Гл. 2. Аппроксимация функций
за скобки коэффициенты при х3 и ж2. Получим
. тпк _ 7тг ж3-F0/7)B/тгJж
8Ш 2	6 Ж2 + 20B/тгJ '
Разделим числитель на знаменатель по правилу деления многочленов и
введем обозначения для коэффициентов. Получим
sin —- = кг [х + —
7тг	200 /2\2
*i = -T *» = -—(-;. *з=2о^-
Полученное выражение можно записать в виде
/	h
sin —— = кг ( ж +
2	V ж + Аз/
Для вычисления значения функции по этой формуле требуется меньше
операций (два деления, два сложения, одно умножение), чем для вычис-
вычисления с помощью выражения B.27) или усеченного ряда Тейлора B.13)
с использованием схемы Горнера (четыре умножения, два сложения).
Следует, однако, иметь в виду, что процессору на выполнение операции
деления обычно требуется гораздо больше времени, чем для выполне-
выполнения операции умножения. Поэтому вопрос об использовании той или
иной аппроксимирующей функции требует в каждом конкретном случае
дополнительного исследования.
Приведем формулы для приближения некоторых элементарных функ-
функций с помощью цепных дробей, указывая интервалы изменения аргумен-
аргумента и погрешности А:
-0.5Ж
+ к2х2
к0 = 1.0000000020967, кг = 0.0999743507186, к2 = 0.0166411490538,
-iln2^^i, |Д|О0-10;
т
In A + х) = к0 +
+
ко = 0.0000000894, кг = 1.0000091365, к2 = 2.0005859000,
к3 = 3.0311932666, к4 = 1.0787748225, къ = 0.1124191908,
0 ^ х ^ 1, А ^ 10;

§ 3. Интерполирование
47
tg —- = х [ ко +
+
/с0 = 0.7853980289, кг = 6.1922344479, к2 = -0.6545887679,
к3 = 0.0020366541, -1 ^ х ^ 1, |А| ^ 2 ¦ 10"8;
,2	\
arctg х = х
х*
+
\
)
к0 = 0.99999752, кг = -3.00064286, к2 = -0.55703890,
/с3 = -17.03715998, к4 = -4.86455143,
-1 ^х ^ 1, |А| ^ 2 • 10~7.
§ 3. Интерполирование
1. Линейная и квадратичная интерполяции. Простейшим и часто
используемым видом локальной интерполяции является линейная (или
кусочно линейная) интерполяция. Она состоит в том, что заданные точ-
точки (xi,yi) (г = 0,1,..., п) соединяются прямолинейными отрезками, и
функция f(x) приближается ломаной с вершинами в данных точках.
Уравнения каждого отрезка ломаной в общем случае разные. По-
Поскольку имеется п интервалов (жг-ъ^г)? то для каждого из них в ка-
качестве уравнения интерполяционного многочлена используется урав-
уравнение прямой, проходящей через две точки. В частности, для г-го
интервала можно написать уравнение прямой, проходящей через точки
Отсюда
У ~ Уг-1
Уг ~ Уг-1
У = CLiX + b{,
Уг - Уг-1
X — Х{ — \
U<i — \ <^ X <^ Х<1)
i = Уг-1 - uiXi-г-
B.28)
Следовательно, при использовании линейной интерполяции сна-
сначала нужно определить интервал, в который попадает значение аргу-
аргумента х, а затем подставить его в формулу B.28) и найти приближенное
значение функции в этой точке. Данный алгоритм представлен на рис. 2.7.
Попытайтесь разобраться, будет ли работать этот алгоритм, если окажет-
окажется, что х < xq или х > хп, и при необходимости модифицировать его.

48
Гл. 2. Аппроксимация функций
Рассмотрим теперь случай квадратичной интерполяции. В качестве
интерполяционной функции на отрезке [a^-b^i+i] принимается квад-
квадратный трехчлен. Такую интерполяцию называют также параболической.
Уравнение квадратного трехчлена
содержит три неизвестных коэффициента а«, 6«, с«, для определения
которых необходимы три уравне-
уравнения. Ими служат условия прохожде-
прохождения параболы B.29) через три точки
хг —15 Уг — \)г> \х1чУ1)г> y^i + li Уг + 1)- ^>1]Л
условия можно записать в виде
&г%г — 1* ^г^г — 1 ~г ^г ~Уг — 1)
CLiX,: + biXi + Ci = yij	B.30)
Если отвлечься от того, что
отрезок [xi-iiXi+i} является под-
подмножеством отрезка [жо,жп], и рас-
рассмотреть его отдельно, то квадра-
квадратичную интерполяцию B.29) можно
трактовать как глобальную с п = 2, а
систему B.30) — как частный случай
системы B.4).
ДО
b
х <
= Уг
Ввод
Xi
-1 —
{^г}, {Уг
г = 0
г = г +
Уг ~ Уг-1
axi-i , у -
Вывод у
} 5 Х
1
= ах + b
Рис. 2.7. Алгоритм линейной интер-
интерполяции
Алгоритм вычисления приближенного значения функции с помощью
квадратичной интерполяции можно записать в виде структурограммы,
как и для случая линейной интерполяции (см. рис. 2.7). Вместо форму-
формулы B.28) нужно использовать B.29) с учетом решения системы линейных
уравнений B.30). Интерполяция для любой точки х е [жо,жп] проводит-
проводится по трем ближайшим к ней узлам.
Пример. Найти приближенное значение функции у = f(x) при х =
= 0.32, если известна следующая таблица ее значений:
0.15
2.17
0.30
3.63
0.40
5.07
0.55
7.78
Воспользуемся сначала формулой линейной интерполяции B.28).
Значение х = 0.32 находится между узлами xi-1 = 0.30 и Xi = 0.40.
В этом случае
= K-tt-i = 5.07-3.63
1 х{ - Хг-г 0.40 - 0.30	" '
bi = Уг-i ~ a+Xi-! = 3.63 - 14.4 • 0.30 = -0.69,
у и 14.4 х - 0.69 = 14.4 ¦ 0.32 - 0.69 = 3.92.
Найдем теперь приближенное значение функции с помощью формулы
квадратичной интерполяции B.29). Составим систему уравнений B.30)

§ 3. Интерполирование	49
с учетом ближайших к точке х = 0.32 узлов: Х{-\ = 0.15, xi = 0.30,
xi+i = 0.40. Соответственно yi-\ = 2.17, yi = 3.63, yi+\ = 5.07. Систе-
Система B.30) запишется в виде
0.152^ + 0.15 Ъ{ + Ci = 2.17,
О.ЗО2^ + 0.30 Ь{ + Ci = 3.63,
0.402^ + 0.40 Ь{ + Ci = 5.07.
Решая эту систему, находим а« = 18.67, hi = 1.33, с« = 1.55. Искомое
значение функции у и 18.67 • 0.322 + 1.33 • 0.32 + 1.55 и 3.89.
2. Многочлен Лагранжа. Перейдем к случаю глобальной интерпо-
интерполяции, т. е. к построению интерполяционного многочлена, единого для
всего отрезка [жо, хп\. Как уже отмечалось в § 1, для построения такого
многочлена достаточно решить систему B.4), однако этот путь неэффек-
неэффективен.
Будем искать интерполяционный многочлен в виде линейной комби-
комбинации многочленов степени п:
Цх) = yolo(x) + yih(x) + ... + ynln(x).	B.31)
При этом потребуем, чтобы каждый многочлен h(x) обращался в нуль
во всех узлах интерполяции, за исключением одного (г -го), где он дол-
должен равняться единице. Легко проверить, что этим условиям при г = 0
отвечает многочлен вида
, / ч (х-х1)(х-х2)...(х-хп)
/о (ж) = т	^7	\	1	\ •	B.32)
[Хо - Xi){Xq — Х2) . . . (Хо - Хп)
Действительно, Iq(xq) = 1. При х = х\, х2, ... , хп числитель
выражения B.32) обращается в нуль. По аналогии с B.32) получим
(х - xq)(x - х2)... (х - хп)
h(x) =
1 - ХО)(Х1 -Х2)... (Ж1 - Хп) '
(х - хр)(х - xi)(x - х3)... (х - хп)
1"П\Х)
(X2 ~ XO)(X2 ~ X\)(X2 - Ж3) . . . (Жз - Ж„) '
	 B.33)
\Xi ^0J . . . \Xi JLi — \)\Xi eL^_|_^J • • • yXi Ju^}
	 7»л ) f 7» 	 7>1 I	I rp 	 rp 1 I
и / V	X / *** V	ть — X у
/	ч /	ч	/	ч •
{xn - xo){xn - x\)... {xn - xn-i)
Подставляя в B.31) выражения B.32), B.33), находим
'
Эта формула определяет интерполяционный многочлен Лагранжа.
4 Л.И. Турчак, П.В. Плотников

50	Гл. 2. Аппроксимация функций
Единственность представления многочлена B.31) в виде B.34)
следует из единственности решения системы B.4).
Из формулы B.34) можно получить выражения для линейной (п = 1)
и квадратичной (п = 2) интерполяций:
X — Xi	X - Хо
Цх) = 	уо Н	2/i, n = 1;
Х Х	Х Х
х-х1)(х-х2) , (жжо)(жж2)
2/ +	\ 2/1 +	B.35)
(х0 -
+ ^42/2
Существует несколько обобщений интерполяционного многочлена Ла-
гранжа. Например, довольно широко используются интерполяционные
многочлены Эрмита. Здесь наряду со значениями функции yi в узлах Х{
задаются значения ее производной у[. Задача состоит в том, чтобы най-
найти многочлен ip(x) степени 2п + 1, значения которого и значения его
производной в узлах Xi удовлетворяют соответственно соотношениям
<Р(Хг)=Уг, <р'(Хг)=у'ц % = 0, 1, . . . , П.
В этом случае также существует единственное решение, если все Xi раз-
различны.
3. Многочлен Ньютона. До сих пор не делалось никаких предполо-
предположений о законе распределения узлов интерполяции. Теперь рассмотрим
случай равноотстоящих значений аргумента, т. е. х\ — х\-\ = h = const
(г = 1,2,...,п). Величина h называется шагом.
Введем также понятие конечных разностей. Пусть известны значения
функции в узлах Xi: yi = f(xi). Составим разности значений функции:
А^/о =У1~Уо= /(ж0 + h) - /(ж0),
Aj/i =У2~У1 = /(ж0 + 2ft) - f(xo + ft),
Aj/n-i =Уп- Уп-i = /(ж0 + raft) - /(ж0 + (п - l)ft).
Эти значения называются первыми разностями (или разностями первого
порядка) функции.
Можно составить вторые разности функции:
Аналогично составляются разности порядка к :
Afc№ = Д*-1!^! - A^yi, г = 0,1,..., п - 1.
Конечные разности можно выразить непосредственно через значения
функции. Например,
А2у0 = Aj/i - А^/о = B/2 - 2/1) - (г/1 - 2/о) = 2/2 - 2j/i + 2/0,
А32/о = AV - А22/О = • • • = Уз ~ 3^/2 + 3yi - 2/0 ¦

§ 3. Интерполирование	51
Аналогично для любого к можно написать
Afe2/o =ук- kyk-i + к^~ ^ Ук-2 + • • • + (-1)^2/0-	B.36)
Эту формулу можно записать и для значения разности в узле xi:
к(к — 1)
A*2/i = Ук+i ~ кук+i-l +	2,	Ук+г-2 + • • • + (-l)*2/t-
Используя конечные разности, можно определить у к :
ук=Уо + кАУо + 2~ А22/о + • • • + Аку0.
Перейдем к построению интерполяционного многочлена Ньютона.
Этот многочлен будем искать в следующем виде:
N(x) = ао + di(x — xq) + а2(х — xq)(x — xi) + ...
... + an(x - xq)(x - xi)... (x - xn-i). B.37)
График многочлена должен проходить через заданные узлы, т. е.
N(x{) = уi (г = 0,1,... ,п). Эти условия используем для нахождения
коэффициентов многочлена:
N(xQ) = а0 = 2/о,
N(xi) = ао + а\(х\ — xq) = ао + a\h = г/i,
N(x2) = ао + d\{x2 — xq) + a2(x2 — xq)(x2 — x\) =
= ао H~ 2oL\lfi -\- 2iCL2h = т/2.
Найдем отсюда коэффициенты ао , а\, а2 :
2/1 - «о 2/1 - 2/о А^/о
«о = 2/о
/г	/г	h '
у2 - clq - 2a\h y2 - уо -
2/г2	2/г2	2/г2 '
Аналогично можно найти и другие коэффициенты. Общая формула имеет
вид
= А*у0 ^ = () х
Подставляя эти выражения в формулу B.37), получаем следующий
вид интерполяционного многочлена Ньютона:
лг/_\	 I ^2/0 /	\
{х - хо)(х -хг)...(х- хп-г). B.38)
Конечные разности Аку0 могут быть вычислены по формуле B.36).

52
Гл. 2. Аппроксимация функций
Формулу B.38) часто записывают в другом виде. Для этого вводится
переменная t = (х — xq) / h; тогда
х = xq + th,
X — ;
Ж — X\
h
х — xq — h
= t-2,
ЮЕ
7V(x) = N(x0 +th)=yo+t Ay0
h
JL JLyi — 1
= t - 1,
= t - n + l.
h --»•¦•> h
С учетом этих соотношений формулу B.38) можно переписать в виде
Д22/о
-+^-1)-nf-n + 1)A",o. B.39)
Полученное выражение называется первым интерполяционным много-
многочленом Ньютона для интерполирования вперед. Оно может аппроксими-
аппроксимировать данную функцию у = f\x) на всем отрезке изменения аргумен-
аргумента [хо, хп]. Однако с точки зрения повышения точности расчетов (путем
уменьшения погрешностей округления) более целесообразно использо-
использовать B.39) для вычисления значении функции в точках левой половины
рассматриваемого отрезка.
Для правой половины отрезка [жо,жп] разности лучше вычислять
справа налево. В этом случае t = (х — xn)/h, т. е. t < 0, и интерпо-
интерполяционный многочлен Ньютона можно получить в виде
N(x) = N(xn +th)=yn+
2!
/n_2 + ...
... + '(«+1)-п(«+)Дп№. B.40)
Полученная формула называется вторым интерполяционным много-
многочленом Ньютона для интерполирования назад.
Рассмотрим пример применения интерполяционной формулы Ньюто-
Ньютона при ручном счете.
Таблица 2. 1
X
0
0.2
0.4
0.6
0.8
1.0
У
1.272
4.465
5.644
5.809
3.961
2.101
Ау
3.193
1.179
0.165
-1.848
-1.860
А2у
-2.014
-1.014
-2.013
-0.012
А3у
1.000
-0.999
2.001
А4у
-1.999
3.000
А5у
4.999
П р и м е р. Вычислить в точке х = 0.1 значение функции у = f(x),
заданной табл. 2. 1.

§ 3. Интерполирование	53
Процесс вычислений удобно свести в ту же табл. 2.1. Каждая после-
последующая конечная разность получается путем вычитания в предыдущей
колонке верхней строки из нижней. При х = 0.1 имеем t = (х — xo)/h =
= @.1 — 0)/0.2 = 0.5. Проводя вычисления с пятью разрядами по форму-
формуле B.39), получим
/@.1) и TV(O.l) = 1.272 + 0.5 • 3.193 +
+ 0-5@-5-1) . BШ4) + 0-5@-5-1К0-5-2) .
+	BШ4) +
+ 0.5@.5 -
+ 0.5@.5-1)@.5-2)@.5-3)@.5-4) 4 ggg =
= 1.272 + 1.597 + 0.2518 + 0.06249 + 0.07806 + 0.1367 = 3.398.
Для сравнения проведем аналогичные вычисления по формуле B.40).
В этом случае t = (х — xn)/h = @.1 — 1)/0.2 = —4.5. Тогда
/@.1) и TV(O.l) = 2.101 - 4.5 • (-1.860) +
+ -4-5(-f + 1) • (-0.012) -45(-45 + 1)(-45 + 2)
+ -4.Б(-4.5 + 1)(-4.Б + 2)(-4.Б + 3) .
4!
-4.5(-4.5 + 1)(-4.5 + 2)(-4.5 + 3)(-4.5 + 4)
4 ggg =
+
= 2.101 + 8.370 - 0.09450 - 13.13 + 7.383 - 1.231 = 3.402.
Видно, что здесь происходит заметная потеря точности. Если проводить
вычисления более точно, то формулы B.39) и B.40) приведут к одному
результату /@.1) и 3.3975.
Мы рассмотрели построение интерполяционного многочлена Ньюто-
Ньютона для равноотстоящих узлов. Можно построить многочлен Ньютона
и для произвольно расположенных узлов, как и в случае многочлена
Лагранжа. Однако этот случай мы рассматривать не будем.
В заключение отметим, что разные способы построения многочленов
Лагранжа и Ньютона дают тождественные интерполяционные формулы
при заданной таблице значений функции. Это следует из единственности
интерполяционного многочлена заданной степени (при отсутствии сов-
совпадающих узлов интерполяции). Многочлен Ньютона удобно применять,
например, если количество узлов интерполяции постепенно увеличива-
увеличивается. При этом учет нового узла требует лишь вычисления одного допол-
дополнительного слагаемого в B.39) и B.40), в то время как при использовании
многочлена Лагранжа требуется пересчитывать все слагаемые в B.34).

54	Гл. 2. Аппроксимация функций
4. Точность интерполяции. График интерполяционного многочлена
у = (р(х) проходит через заданные точки, т. е. значения многочлена и дан-
данной функции у = f(x) совпадают в узлах х = Х{ (г = 0,1,..., п). Если
функция f(x) сама является многочленом степени п, то имеет место
тождество f(x) = (р(х). В общем случае в точках, отличных от узлов
интерполяции, R(x) = f(x) — (р(х) ф 0. Эта разность есть погрешность
интерполяции и называется остаточным членом интерполяционной фор-
формулы. Оценим его значение.
Предположим, что заданные числа у{ являются значениями некоторой
функции у = f(x) в точках х = Х{. Пусть эта функция непрерывна и
имеет непрерывные производные до (п + 1) - го порядка включительно.
Можно показать, что в этом случае остаточный член интерполяционного
многочлена Лагранжа имеет вид
=(x-x0)(x-Xl) (x-xn)f
=	f
(п + 1).
Здесь /(п+1^(ж*) — производная (п + 1) - го порядка функции f(x)
в некоторой точке х = ж* , ж* е [жо,жп]. Если максимальное значение
этой производной равно
max |/("+1)(Ж)|=М„+1,
то можно записать формулу для оценки остаточного члена:
где функция ujn (x) определена как
ujn(x) = (х - хо)(х - xi)... (х - хп).	B.42)
Проанализировав поведение функции ujn(x), можно сделать вывод о
том, что погрешность интерполяции Rl (x) в среднем будет тем выше,
чем ближе точка х лежит к концам отрезка [xq ,xn]. Если же использовать
интерполяционный многочлен для аппроксимации функции вне отрез-
отрезка [жо, хп] (экстраполяция), то погрешность возрастет существенно.
Вид остаточного члена интерполяционного многочлена Ньютона в
случае равноотстоящих узлов можно легко получить из B.41):
Если предположить, что разность Ап+1уп постоянна, то можно запи-
записать следующую формулу остаточного члена первого интерполяцион-
интерполяционного многочлена Ньютона:
Следует еще раз подчеркнуть, что существует один и только один
интерполяционный многочлен при заданном наборе узлов интерполя-
интерполяции. Формулы Лагранжа, Ньютона и другие порождают один и тот же

§ 3. Интерполирование	55
многочлен (при условии, что вычисления проводятся точно). Разница
лишь в алгоритме их построения. Правда, интерполяционный многочлен
Лагранжа не содержит явных выражений для коэффициентов.
Выбор способа интерполяции определяется различными соображени-
соображениями: точностью, временем вычислений, погрешностями округлений и др.
В некоторых случаях более предпочтительной может оказаться локальная
интерполяция, в то время как построение единого многочлена высокой
степени (глобальная интерполяция) не приводит к успеху.
Такого рода ситуацию в 1901 г. обнаружил К. Рунге. Он строил на отрезке
— 1 ^ х ^ 1 интерполяционные многочлены с равномерным распределени-
распределением узлов для функции у = 1/A + 25ж2). Оказалось, что при увеличении
степени интерполяционного многочлена последовательность его значений
расходится для любой фиксированной точки х при 0.7 < ж < 1.
Положение в некоторых случаях может быть исправлено специаль-
специальным распределением узлов интерполяции (если они не зафиксированы).
Доказано, что если функция f(x) имеет непрерывную производную на
отрезке [—1,1], то при выборе значений Х{, совпадающих с корнями
многочленов Чебышева степени п + 1, интерполяционные многочлены
степени п сходятся к значениям функции в любой точке этого отрезка.
Наглядно пояснить сделанное утверждение можно следующим образом.
Как было отмечено в § 2, корни многочленов Чебышева расположены
неравномерно на отрезке и сгущаются к его концам. Такое сгущение
компенсирует увеличение погрешности интерполяции при приближении
к концам отрезка, которое имеет место при равномерном распреде-
распределении узлов.
Таким образом, повышение точности интерполяции целесообразно
производить за счет уменьшения шага и специального расположения то-
точек xi. Повышение степени интерполяционного многочлена при локаль-
локальной интерполяции также уменьшает погрешность, однако здесь не все-
всегда ясно поведение производной f(n+1\x) при увеличении п. Поэтому
на практике стараются использовать многочлены малой степени (линей-
(линейную и квадратичную интерполяции, сплайны).
5. Сплайны. Сейчас широкое распространение для интерполяции
получило использование кубических
сплайн-функций — специальным образом
построенных многочленов третьей сте-
степени. Они представляют собой некото-
некоторую математическую модель гибкого
тонкого стержня из упругого материала.
Если закрепить его в двух соседних уз-
узлах интерполяции с заданными углами
наклонов а и /3 (рис. 2.8), то между
точками закрепления этот стержень (ме- Рис. 2.8. Механический сплайн
ханический сплайн) примет некоторую
форму, минимизирующую его потенциальную энергию.

56	Гл. 2. Аппроксимация функций
Пусть форма этого стержня определяется функцией у = S(x). Из
курса сопротивления материалов известно, что уравнение свободного
равновесия имеет вид 5IV (х) = 0. Отсюда следует, что между каждой па-
парой соседних узлов интерполяции функция S(х) является многочленом
степени не выше третьей. Запишем ее в виде
V ) - х\ ) ~ г г\	г-Ц гУ.	г-1)	^ ^
+ уЧ . / /уг 	 /у* .	\	Т* • -t ^^ Т* ^^ Т* •
Для определения коэффициентов ai, bi, Ci, ^ на всех п элементар-
элементарных отрезках необходимо получить An уравнений. Часть из них вытекает
из условий прохождения графика функции S(х) через заданные точки,
т. е. Si(xi-i) = yi-\, Si(xi) = yi. Эти условия можно записать в виде
SiiXi-!) = СЦ = Уг-и	B.44)
Si(xi) = cti + bih + c«/i2 + rfi/г3 = 2/i,	B.45)
hi = Xi — Xi-i, i = 1, 2,..., n.
Эта система содержит 2п уравнений. Для получения недостающих урав-
уравнений зададим условия непрерывности первых и вторых производных во
внутренних узлах интерполяции, т. е. условия гладкости второго порядка
кривой во всех точках:
Qf /	\ 	 Qf	/	\	Qff 1	\ 	 Qll	/	\	• 	 1О	п-,	1	/О Л?\\
Di\xi) — Di+l\xi)i Di \xi) — Di+l\xi)i Z — ±, Z, . . . , П — ±. {/,.40)
Вычислим производные многочлена B.43):
Si(x) = bi + 2a(x - xi-г) + 3di(x - х{-г) ,	^^
Si (x) = 2ci + бс/дж — Xi-i).
Подставляя эти выражения в B.46), получаем 2п — 2 уравнений
c«+i = a -\- Shidi,	B.49)
i = 1, 2,... ,n — 1.
Недостающие два уравнения получаются из условий закрепления кон-
концов сплайна. Обычно эти условия представляют собой соотношения, в
которые входят значения первой и второй производных функции S(х)
в точках xq и хп . Поэтому указанные значения должны входить в рас-
рассматриваемую систему уравнений. Из B.47) следует, что S[(xi-i) = bi,
Si(xi-i) = 2ci. Отсюда S[(xq) = 6i, 5"(#o) = 2ci, т. е. значения про-
производных в точке xq присутствуют в системе. Значения же производных
в точке хп в системе отсутствуют. Введем их в систему с помощью
дополнительных неизвестных 6n+i и cn+i:
Ь \Хп) = on+i, о \хп) = zcn_^i.
Из условий непрерывности производных в точке хп следует, что
6n+i = bn + 2/incn + 3h^dn, cn+i = cn + 3hndn.

§ 3. Интерполирование	57
Таким образом, соотношения B.48), B.49) можно рассматривать для диа-
диапазона индексов
г = 1,2,..., п.	B.50)
Система B.44), B.45), B.48), B.49) с учетом B.50) содержит 4п+2
неизвестных и An уравнений и может быть дополнена, например, сле-
следующими условиями закрепления концов сплайна:
S'(x0) = Ьг = ки S'(xn) = Ъп+1 = к2	B.51)
или
S"(x0) = 2d = пц, S"(xn) = 2cn+1 = m2,	B.52)
где ki, /c2 , mi, rri2 — заданные числа.
В частности, при заданных углах наклона а и /3 (см. рис. 2.8)
кг =tga, k2 = tg/3.
При свободном закреплении концов можно приравнять нулю кривизну
линии в точках закрепления. Получаемая таким образом функция назы-
называется свободным кубическим сплайном. Из условия нулевой кривизны
на концах следуют равенства нулю вторых производных в этих точках.
Отсюда
mi = rri2 = 0.
Соотношения B.44), B.45), B.48), B.49), а также B.51) или B.52) со-
составляют систему линейных алгебраических уравнений для определения
коэффициентов п{, d{ (i = 1, 2,..., п) и hi, с« (i = 1, 2,..., п + 1). Ее
можно решить одним из методов, изложенных в гл. 4.
Однако с целью экономии памяти компьютера и машинного времени
эту систему можно привести к более удобному виду. Из условия B.44)
сразу можно найти все коэффициенты а«. Далее, из B.49) получим
C\ г = 1,2,...,п.	B.53)
Подставим эти соотношения, а также значения а« = у{-\ в B.45) и най-
найдем отсюда коэффициенты
bi = yi-yi-1 _fc?( +2c<) i = 12,...,n.	B.54)
hi	6
Учитывая выражения B.53) и B.54), исключаем из уравнения B.48) ко-
коэффициенты di и hi. Окончательно получим следующую систему урав-
уравнений только для коэффициентов с«:
2(fti_1 + ftOci + /ВД+1 = 3 (
\ H	ti	B.55)
г = 2,3,..., n;
=fc2 B.56)

58	Гл. 2. Аппроксимация функций
или
ci = ^, сп+1 = ^.	B.57)
Здесь уравнения B.56) используются при применении условий B.51),
а уравнения B.57) — при применении условий B.52).
Матрица системы B.55) трехдиагональная, т. е. ненулевые элементы
находятся лишь на главной и двух соседних с ней диагоналях, распо-
расположенных сверху и снизу. Для ее решения целесообразно использовать
метод прогонки (см. гл. 4). По найденным из системы B.55), B.56) или
B.57) коэффициентам с« легко вычислить коэффициенты d{, hi.
Заметим, что кубическая сплайн-функция, удовлетворяющая услови-
условиям B.51) или B.52), обладает наименьшей (в некотором смысле) кри-
кривизной среди всех дважды непрерывно дифференцируемых функций на
отрезке [жо,жп] с заданными значениями в узлах интерполяции, удов-
удовлетворяющих условиям B.51) или B.52). А именно,
[f"(x)fdx
для всех f(x) из указанного класса функций.
6. О других формулах интерполяции. Ранее уже упоминалась одна
из модификаций многочлена Лагранжа — интерноляционный многочлен
Эрмита. При построении этого многочлена требуется, чтобы в узлах Х{
совпадали с табличными данными не только его значения, но и значения
его производных до некоторого порядка. В общем случае выражение
для многочлена Эрмита очень громоздко, и пользоваться им на прак-
практике трудно. Поэтому ограничиваются лишь некоторыми простейшими
случаями. Например, многочлен Эрмита, который сохраняет в двух точ-
точках (х = хо, х\) значения заданной функции у = f(x) и ее первой
производной у' = f'(x), имеет вид
тт{ \	, /	\ / / , х~х0 \(	2/0 -2/1
Н(х) =уо + (х- х0) < 2/о +	 2/о
[ х x |Д xx
Хо ~ X!
Иногда при выводе интерполяционных формул удобнее использовать
не односторонние разности, как для многочлена Ньютона, а центральные.
На этом основаны интерполяционные формулы Стирлинга и Бесселя.
Они могут быть получены путем преобразования формулы Ньютона.
Рассмотрим интерполирование функций специального вида, а именно
периодических функций. Для функции с периодом 2тг можно построить

§ 3. Интерполирование	59
интерполяционную формулу по аналогии с формулой Лагранжа:
— х\) sin(x — Х2) • • • sin(x — хп)
Т/0 Н
, . s
Г [X) = —
;	Т/0
"" Х2) • • • Sin(#o "" хп)
sin(x — Х2) • • • sin(x — хп
2/1
+ г7тг7чг77
sin(#i — жо) sm(xi — Х2)... sm(xi — хп)
sin(x — жо) sin(x — Ж1)... sin(x — xn_i)
sin(xn — жо) sin(xn — Ж1)... sin(xn — xn_i)
7. Функции двух переменных. До сих пор мы рассматривали ин-
интерполирование функций одной независимой переменной у = f(x). На
практике возникает также необходимость построения интерполяционных
формул для функций нескольких переменных. Для простоты ограничим-
ограничимся функцией двух переменных z = /(ж, у). Пусть ее значения заданы на
множестве равноотстоящих узлов (xi,yi) (г, j = 0,1, 2). Введем обозна-
обозначения Zij = f(xi,yj), fti = xi+1 -xi9 h2= 2/7+1 - Уз .
Построим многочлен, аналогичный многочлену Ньютона для случая
одной переменной. Здесь нужно вычислять разности двух видов — по на-
направлениям х и у. Эти частные разности первого порядка определяются
формулами
Запишем также выражения для частных разностей второго порядка:
Интерполяционный многочлен второй степени можно записать в виде
х — xq у — 2/о
F[x, у) = ?оо "I	г	Аж^оо Н	г— А^^оо +
(х-хо)(х-х1) 2	(ж-жо)B/-2/о) д2 „ ,
2/^
(У-Уо)(у-У1) А2
Можно также построить линейную интерполяционную формулу. Гео-
Геометрически это означает, что нужно найти уравнение плоскости, про-
проходящей через три заданные точки (x«,2/i,^), (i = 1, 2, 3), где zi =
= f(xi, y^ . Из курса аналитической геометрии известно, что уравнение
плоскости, проходящей через три точки, можно записать в виде
х-xi у - 2/1 2 - z\
Х2 — Х1 2/2 "" 2/1 ^2 "" 2i
?з - xi Уз- У\ ^з - zx
= 0.

60
Гл. 2. Аппроксимация функций
Отсюда можно найти
z= —(D0-D1x-D2y),
B.58)
Х\
Х2
хз
XI
Х2
Хз
У\
У2
Уз
1
1
1
Z\
%2
Z3
Z\
Z2
Z3
, D\ =
, D3 =
1
1
1
xi
Х2
хз
У\
2/2
Уз
У\
2/2
Уз
Z\
%2
Z3
1
1
1
Do =
D2 =
Пример. Вычислить приближенное значение функции z = f(x,y)
в точке A,0), если известны ее значения z\ = /@,0) = 0, ?2 = /B,4) =
= -3, *з = /D,-2) = 1.
Воспользуемся формулой линейной интерполяции B.58). Вычислим
значения определителей
= -2,
= -20.
Таким образом, z и — Bж + 14г/)/20, или 2; и —0.1ж — 0.72/. Это и есть
формула линейной интерполяции для данного примера. При ж = 1, 2/ = 0
получим z и —0.1.
0
2
4
0
2
4
0
4
-2
1
1
1
0
-3
1
0
-3
1
= 0, D1 =
= -14, D3 =
1
1
1
0
2
4
0
4
-2
0
4
-2
0
-3
1
1
1
1
§ 4. Подбор эмпирических формул
1. Характер опытных данных. При интерполировании функций мы
использовали условие равенства значений интерполяционного много-
многочлена и данной функции в известных точках — узлах интерполяции.
Это предъявляет высокие требования к точности данных значений функ-
функции. В случае обработки опытных данных, полученных в результате
наблюдений или измерений, нужно иметь в виду ошибки этих дан-
данных. Они могут быть вызваны несовершенством измерительного при-
прибора, субъективными причинами, различными случайными факторами
и т. д. Ошибки экспериментальных данных можно условно разбить на
три категории по их происхождению и величине: систематические, слу-
случайные и грубые.
Систематические ошибки обычно дают отклонение в одну сторону от
истинного значения измеряемой величины. Они могут быть постоянны-
постоянными или закономерно изменяться при повторении опыта, и их причина и
характер известны. Систематические ошибки могут быть вызваны усло-
условиями эксперимента (влажностью, температурой среды и др.), дефектом

§ 4. Подбор эмпирических формул	61
измерительного прибора, его плохой регулировкой (например, смеще-
смещением указательной стрелки от нулевого положения) и т. д. Эти ошибки
можно устранить наладкой аппаратуры или введением соответствующих
поправок.
Случайные ошибки определяются большим числом факторов, кото-
которые не могут быть устранены либо достаточно точно учтены при измере-
измерениях или при обработке результатов. Они имеют случайный (несистема-
(несистематический) характер, дают отклонения от средней величины в ту и другую
стороны при повторении измерений и не могут быть устранены в экс-
эксперименте, как бы тщательно он ни проводился. С вероятностной точки
зрения математическое ожидание случайной ошибки равно нулю. Стати-
Статистическая обработка экспериментальных данных позволяет определить
величину случайной ошибки и довести ее до некоторого приемлемого
значения повторением измерений достаточное число раз.
Грубые ошибки явно искажают результат измерения; они чрезмерно
большие и обычно пропадают при повторении опыта. Грубые ошибки
существенно выходят за пределы случайной ошибки, полученные при
статистической обработке. Измерения с такими ошибками отбрасыва-
отбрасываются и в расчет при окончательной обработке результатов измерений
не принимаются.
Таким образом, в экспериментальных данных всегда имеются слу-
случайные ошибки. Они, вообще говоря, могут быть уменьшены до сколь
угодно малой величины путем многократного повторения опыта.
Однако это не всегда целесообразно, поскольку могут потребоваться
большие материальные или временные ресурсы. Значительно дешевле
и быстрее можно в ряде случаев получить уточненные данные хорошей
математической обработкой имеющихся результатов измерений.
В частности, с помощью статистической обработки результатов изме-
измерений можно найти закон распределения ошибок измерений, наиболее
вероятный диапазон изменения искомой величины (доверительный ин-
интервал) и другие параметры. Рассмотрение этих вопросов выходит за
рамки данного пособия; их изложение можно найти в некоторых книгах,
приведенных в списке литературы. Здесь ограничимся лишь определе-
определением связи между исходным параметром х и искомой величиной у на
основании результатов измерений.
2. Эмпирические формулы. Пусть, изучая неизвестную функцио-
функциональную зависимость между у и х, мы в результате серии экспериментов
произвели ряд измерений этих величин и получили таблицу значений
Хо
Уо
XI
2/1
Х-п
Уп
Задача состоит в том, чтобы найти приближенную зависимость
У = Я*),	B.59)

62	Гл. 2. Аппроксимация функций
значения которой при х = xi (г = 0,1,...,п) мало отличаются от
опытных данных yi. Приближенная функциональная зависимость B.59),
полученная на основании экспериментальных данных, называется эмпи-
эмпирической формулой.
В § 1 отмечалось, что задача построения эмпирической формулы отли-
отличается от задачи интерполирования. График эмпирической зависимости,
вообще говоря, не проходит через заданные точки (ж^ 2/г), как в случае
интерполяции. Это приводит к тому, что экспериментальные данные
в некоторой степени сглаживаются, в то время как интерполяционная
формула повторила бы все ошибки, имеющиеся в экспериментальных
данных.
Построение эмпирической формулы состоит из двух этапов: подбора
общего вида этой формулы и определения наилучших значений содер-
содержащихся в ней параметров. Общий вид формулы иногда известен из
физических соображений. Например, для упругой среды связь между
напряжением а и относительной деформацией е определяется законом
Гука: а = Ее, где Е — модуль упругости; задача сводится к опре-
определению одного неизвестного параметра Е.
Если характер зависимости неизвестен, то вид эмпирической форму-
формулы может быть произвольным. Предпочтение обычно отдается наиболее
простым формулам, обладающим достаточной точностью. Они первона-
первоначально выбираются из геометрических соображений: эксперименталь-
экспериментальные точки наносятся на график, и примерно угадывается общий вид
зависимости путем сравнения полученной кривой с графиками извест-
известных функций (многочлена, показательной или логарифмической функ-
функций и т. п.).
Успех здесь в значительной мере определяется опытом и интуицией
исследователя.
Простейшей эмпирической формулой является линейная зависимость
у = ах + Ь.	B.60)
Близость экспериментального распределения точек к линейной зависи-
зависимости легко просматривается после построения графика данной экспери-
экспериментальной зависимости. Кроме того, эту зависимость можно проверить
путем вычисления значений ki:
ki = Ayi/Axi, Ayi = yi+1 - уи Ахг = xi+1 - хи
i = 0,1,... ,п - 1.
Если при этом к{ и const, то точки (х^ у г) расположены приблизитель-
приблизительно на одной прямой, и может быть поставлен вопрос о применимости
эмпирической формулы B.60). Точность такой аппроксимации определя-
определяется отклонением величин ki от постоянного значения. В частном случае
равноотстоящих точек Xi (т. е. Axi = const) достаточно проверить по-
постоянство разностей А^.

§ 4. Подбор эмпирических формул	63
Пример. Проверим возможность использования линейной зависи-
зависимости для описания следующих данных:
X
У
0
1.17
0.5
1.81
1.0
2.50
1.5
3.15
2.0
3.79
2.5
АЛА
Поскольку здесь х\ —равноотстоящие точки (Ах{ = a^+i — хг — 0.5), то
достаточно вычислить разности Ayi: 0.64, 0.69, 0.65, 0.64, 0.65. Так как
эти значения близки друг к другу, то в качестве эмпирической формулы
можно принять линейную зависимость.
В ряде случаев к линейной зависимости могут быть сведены и другие
экспериментальные данные, когда их график в декартовой системе ко-
координат не является прямой линией. Это может быть достигнуто путем
введения новых переменных ?, г] вместо х, у:
? = (р(х,у), Г1 = ф(х,у).	B.61)
Функции (р(х,у) и г] = ф(х,у) выбираются такими, чтобы точ-
точки (&, rji) лежали на некоторой прямой линии в плоскости (?, rj). Такое
преобразование называется выравниванием данных.
Для получения линейной зависимости
г] = а? + Ъ
с помощью преобразования B.61) исходная формула должна быть запи-
записана в виде
ф(х,у) = CL(f(x,y) +6.
К такому виду легко сводится, например, степенная зависимость у = ахъ,
(х > 0, у > 0). Логарифмируя эту формулу, получаем lgу = Ьlg x + lg a.
Полагая \ = lg х, г] = lg у, находим линейную связь: г] = h\ + с (с = lg a).
Другой простейшей эмпирической формулой является квадратный
трехчлен
у = ах2 + Ьх + с.	B.62)
Возможность использования этой формулы для случая равноотстоящих
точек Xi можно проверить путем вычисления вторых разностей А2^ =
= yi+2 — tyi+i +yi. При A2yi и const в качестве эмпирической формулы
может быть выбрана B.62).
3. Определение параметров эмпирической зависимости. Будем
считать, что тип эмпирической формулы выбран, и ее можно представить
в виде
у = (р(х, а0, аь ..., ат),	B.63)
где (р — известная функция, а0, а1:..., аш — неизвестные постоянные
параметры. Задача состоит в том, чтобы определить такие значения этих
параметров, при которых эмпирическая формула дает хорошее прибли-
приближение данной функции, значения которой в точках Х{ равны у{ (г =
= 0,1,...,п).

64	Гл. 2. Аппроксимация функций
Как уже отмечалось выше, здесь не ставится условие (как в случае
интерполяции) совпадения опытных данных щ со значениями эмпири-
эмпирической функции B.63) в точках Х{. Разность между этими значениями
(отклонения) обозначим через ?«. Тогда
Si = (p(xi, а0, аъ ..., ат) -уи г = 0,1,..., п.	B.64)
Задача нахождения наилучших значений параметров ао, ai,..., аш сво-
сводится к некоторой минимизации отклонений б{. Существует несколько
способов решения этой задачи.
Простейшим из них является метод выбранных точек. Он состоит
в следующем. По заданным экспериментальным значениям на коорди-
координатной плоскости OXY наносится система точек. Затем проводится
простейшая плавная линия (например, прямая), которая наиболее близ-
близко примыкает к данным точкам. На этой линии выбираются точки, ко-
которые, вообще говоря, не принадлежат исходной системе точек. Число
выбранных точек должно быть равным количеству искомых параметров
эмпирической зависимости. Координаты этих точек (х®,у®) тщательно
измеряются и используются для записи условия прохождения графика
эмпирической функции B.63) через выбранные точки:
ip(xp a0, ai,..., ат) = ?/•, j = 0,1,..., га.	B.65)
Из этой системы уравнений находятся значения параметров ao, ai,...
• • • 5«m-
В частности, если в качестве эмпирической формулы принята линей-
линейная зависимость у = ах+Ъ, то на этой прямой выбираются точки (х$, у^)
и (#5, у\), и уравнения B.65) примут вид
ах^ + 6 = 2/2,
п	п	B-66)
ах°1+Ъ = у°1.
Можно также сразу записать уравнение прямой, проходящей через эти
выбранные точки. В этом случае не нужно решать систему B.66)
Рассмотрим еще один способ определения параметров эмпиричес-
эмпирической формулы — метод средних. Он состоит в том, что параметры ао,
ai, ... , аш зависимости B.63) определяются с использованием условия
равенства нулю суммы отклонений B.64) во всех точках xf.
ivfai a<b ai, ¦ ¦ ¦, Q>m) - Vi] = 0-	B.67)
г=0	г=0
Полученное уравнение служит для определения параметров a0, ai,...
..., аш . Ясно, что из одного уравнения нельзя однозначно определить
все га + 1 параметров. Однако, поскольку других условий нет, равен-
равенство B.67) путем группировки отклонений е\, разбивается на систему,
состоящую из га + 1 уравнений. Например,

§ 4. Подбор эмпирических формул
65
е3 + е4
?2 = О,
?6 = О,
?n-i +еп=О.
Решая эту систему уравнений, можно найти неизвестные параметры.
Пример. Тело, движущееся прямолинейно с неизвестной скоростью
г>о, в момент времени t = 0 начинает тормозить. Измерялось расстояние
х от начала торможения в следующие моменты времени t:
t, с
ж, м
0
0
5
106
10
182
15
234
20
261
25
275
Считая движение тела равнозамедленным с постоянным замедле-
замедлением —а (т. е. с ускорением а), найти приближенные значения пара-
параметров г>о и а.
Решение. Искомые параметры могут быть найдены из уравнения
движения тела, которое представим с помощью эмпирической формулы,
используя результаты измерений. Вид эмпирической формулы в данном
случае известен из физических соображений — при равнозамедленном
движении тела пройденное расстояние является квадратичной функцией
времени:
х = At2 + Bt + С.
Легко установить, что G = 0, поскольку х = 0 при t = 0. Эмпирическая
формула принимает вид	о
х = At2 + Bt.	B.68)
Для определения параметров А, В нужно получить два уравнения. Вос-
Воспользуемся методом средних и запишем уравнение B.67) для всех точек
(кроме начальной):
?]_ + ?2 + ?з + ?4 + ?5 = 0.
Запишем вместо этого уравнения систему двух уравнений путем его рас-
расщепления:
?\ + ?2 + S3 = 0,
г4 + еъ = 0.
Используя выражение B.68) и табличные данные, получаем
(А • 52 + В • 5 - 106) + (А • 102 + В • 10 - 182) +
+ (А • 152 + В • 15 - 234) = 0,
(А • 202 + В • 20 - 261) + (А • 25 + В • 25 - 275) = 0.
Или окончательно
1025А + 45Б = 536.
Решая эту систему уравнений, находим А = -0.30, В = 39.07.
5 Л.И. Турчак, П.В. Плотников

66	Гл. 2. Аппроксимация функций
Следовательно, эмпирическую формулу B.68), которая дает прибли-
приближенную связь между пройденным расстоянием и временем, можно запи-
записать в виде
х = -0.30*2 + 39.07*.
Сравнивая это уравнение с уравнением х = at2/2 -\-vot, получаем оценки
для среднего ускорения тела и его начальной скорости:
а = 2А = -0.60 м/с2, vo= В = 39.07 м/с.
Рассмотренные методы определения параметров эмпирической фор-
формулы являются сравнительно простыми, однако в ряде случаев получае-
получаемые с их помощью аппроксимации не обладают достаточной точностью.
4. Метод наименьших квадратов. Запишем сумму квадратов от-
отклонений B.64) для всех точек ж0, #ь ..., хп:
S = ^?2 = ^2 [<p(xi,ao,au...,am) - Уг}2.	B.69)
г=0	г=0
Параметры ао, ai,..., аш эмпирической формулы B.63) будем находить
из условия минимума функции S = 5(ao, ai,..., аш). В этом состоит
метод наименьших квадратов.
В теории вероятностей доказывается, что полученные таким методом
значения параметров наиболее вероятны, если отклонения Е{ подчиня-
подчиняются нормальному закону распределения.
Поскольку здесь параметры a0, ai,..., аш выступают в роли неза-
независимых переменных функции S, то ее минимум найдем, приравнивая
нулю частные производные по этим переменным:
|^0, 1^=0, ..., |^-=0.	B.70)
Полученные соотношения — система уравнений для определе-
определения ao,ai,... ,am.
Рассмотрим применение метода наименьших квадратов для широко
используемого на практике частного случая, когда функция <р является
линейной по неизвестным параметрам ао , а\, ... , аш:
<p(#,ao,ai,... , аш) =
з=о
где ipo,(pi,---, Фт — известные функции х. Формула B.69) для опреде-
определения суммы квадратов отклонений S примет вид
ТЬ	ТП	г»
г=0 j=0

§ 4. Подбор эмпирических формул
67
Для составления системы B.70) продифференцируем S по перемен-
переменным аи (к = 0,1,... ,га):
~ Уг
[
п	m
^[
г=0
3=0
Приравнивая найденные производные нулю, получим следующую
систему уравнений:
B.71)
3=0
Система B.71) является системой линейных алгебраических уравне-
уравнений, ее можно записать в наглядном векторно-матричном виде (см. гл. 4).
Для этого введем векторы опытных данных у и неизвестных парамет-
параметров а, а также матрицу Ф следующим образом
У =
а =
dl
(Р1(хг)
Здесь векторы у и а имеют размерности п + 1 и га + 1 соответственно,
а матрица Ф имеет размерность (п + 1)х(га + 1). Для ее элементов
справедливо выражение
()
Нетрудно убедиться, что выражение, стоящее в квадратных скобках
в B.71), является г-й компонентой вектора Фа — у, а каждое уравне-
уравнение B.71) представляет собой равенство нулю к-й компоненты векто-
вектора Фт(Фа - у), где Фт — транспонированная матрица. Таким образом,
систему B.71) можно записать в виде
или
Фт(Фа-у) = 0
(ФтФ)а = Фту.
B.72)
Матрица этой системы ФТФ имеет размерность (га + 1) х (га + 1), век-
вектор а является искомым.
Пример. Используя метод наименьших квадратов, вывести эмпи-
эмпирическую формулу для функции у = /(ж), заданной в табличном виде:
X
У
0.75
2.50
1.50
1.20
2.25
1.12
3.00
2.25
3.75
4.28

68
Гл. 2. Аппроксимация функций
Решение. Если изобразить данные табличные значения на графике
(рис. 2.9), то легко убедиться, что в качестве эмпирической формулы для
аппроксимации функции y = f(x) можно принять квадратный трехчлен,
графиком которого является парабола:
у и (р(х) = а0 + а\х + а2х2.
В данном случае имеем
т = 2, п = 4,
У =
/2.50\
1.20
1.12
2.25
а =
Ф =
(\
1
1
1
U
0.75
1.50
2.25
3.00
3.75
0.75 \
1.502
2.252
З.ОО2
3.752У
После вычисления матрицы ФТФ и вектора Фту (приведены округлен-
округленные значения)
/ 5
ФТФ = 11.25
V30.94
11.25
30.94
94.92
30.94 >
94.92
309.76;
Фту = 29.00
\90.21>
система уравнений B.72) принимает вид
5а0 + 11.25а! + 30.94а2 = 11.35,
11.25а0 + 30.94ai + 94.92а2 = 29.00,
30.94а0 + 94.92ai + 309.7ба2 = 90.21.
Отсюда находим значения параметров эмпирической формулы: а0 =
= 4.82, ai = —3.88, а2 = 1.00. Таким образом, получаем следующую
аппроксимацию функции, заданной в таб-
ж	личном виде:
4
Оценим относительные погрешности
полученной аппроксимации в заданных точ-
точках, т. е. найдем значения
1	2
Рис. 2.9
Si
—
Уг
(f(Xj)
	¦
Уг

§ 4. Подбор эмпирических формул
69
Результаты вычислений представим в виде таблицы
X
0.75
1.50
2.25
3.00
3.75
(р(х)
2.47
1.25
1.15
2.17
4.32
У
2.50
1.20
1.12
2.25
4.28
?
-0.03
0.05
0.03
-0.08
0.04
ду
-0.012
0.042
0.027
-0.036
0.009
На рис. 2.9 построен график найденной эмпирической функции. Точ-
Точками, как уже отмечалось, нанесены заданные табличные значения функ-
функции у = f(x).
В заключение сделаем несколько замечаний, касающихся метода наи-
наименьших квадратов.
Замечание 1. Можно доказать, что если столбцы матрицы Ф ли-
линейно независимы, то система B.72) имеет единственное решение.
Замечание 2. Как видно из рассмотренного примера, матрица ФТФ
является симметрической, т. е. (ФТФ)^ = (ФТФ)^ (i,j = 0,1,...,га).
Замечание 3. В случае использования в качестве эмпирической
функции многочлена
(р(х) = а0 + а\х + ... + ашхт
элементы матрицы ФТФ и компоненты вектора Фту можно вычислить
по формулам
п	п
E4+J	$>	i,j = 0,l,...,m. B.73)
k=0
k=0
В частности, равны все элементы (ФТФ)^ при г + j = const.
5. Локальное сглаживание данных. Как отмечалось в п. 1, опыт-
опытные данные содержат случайные ошибки, что является причиной раз-
разброса этих данных. Во многих случаях бывает целесообразно провести
их сглаживание для получения более плавного характера исследуемой
зависимости. Существуют различные способы сглаживания. Рассмотрим
один из них, основанный на методе наименьших квадратов.
Пусть в результате экспериментального исследования зависимос-
зависимости у = f(x) получена таблица значений искомой функции уо, г/i, ...
... , уп в точках xq, х\, ... , хп. Значения аргумента xi предполагают-
предполагаются равноотстоящими, а опытные данные щ — имеющими одинаковую
точность. Предполагается также, что функция у = f(x) на произвольной
части отрезка [жо,жп] может быть достаточно хорошо аппроксимирована
многочленом некоторой степени т.
Рассматриваемый способ сглаживания состоит в следующем. Для
нахождения сглаженного значения у{ в точке Х{ выбираем по обе сторо-
стороны от нее к + 1 значение аргумента из имеющихся в таблице (к четное):

70	Гл. 2. Аппроксимация функций
Xi-k/2 >••• , xi-i, Хг, ж»+1 ? • • • ? ^i+fe/2 • По опытным значениям рассмат-
рИВаеМОЙ фуНКЦИИ В ЭТИХ ТОЧКаХ уг-к/2 , • • • , Уг-1 , Уг , Уг+1 , • • • , Уг+к/2
строим многочлен степени га с помощью метода наименьших квадратов
(при этом га ^ к). Значение полученного многочлена щ в точке xi и
будет искомым (сглаженным) значением. Процесс повторяется для всех
внутренних точек. Сглаживание значений, расположенных вблизи концов
отрезка [жо, хп], производится с помощью крайних точек.
Опыт показывает, что сглаженные значения щ, как правило, с до-
достаточной степенью точности близки к истинным значениям. Иногда
сглаживание повторяют. Однако это может привести к существенному
искажению истинного характера рассматриваемой функциональной за-
зависимости.
Приведем в заключение несколько формул для вычисления сгла-
сглаженных значений опытных данных при различных га, к :
га = 1:
Уг = ^(Уг-1 +Уг+ 2/t+l)> к = 2,
Уг = тЫ-2 + Уг-1 + Уг + Уг+1 + 2/t+2), & = 4,
5
Уг = -=(У{-3 + Уг-2 + 2/t-l + Уг + 2/i+l + 2/г+2 + 2/г+з), fc = 6;
га = 3:
Уг = ^(-^Уг-2 + 122/i_i + 17уг + 12^+1 - 3^+2), А^ = 4,
2/t = ^т(~22/г-3 + 3^_2 + 62/i_i + lyi + 6^+1 + 3^+2 - 2^+3), /С = б,
Bl + 14 + 3% + 54_! + 5% + 54yi+1 +
+2 + 14^+з - 21у{+4), к = 8;
га = 5:
Уг = ^
к = 6,
;+2 - 557/i+3 + 1%+4), /С = 8.
Упражнения
1. Записать алгоритмы вычислений с помощью разложений в ряды значений
функций: а) у = cos ж; б) у = е~ж; в) у = shz; г) у = л/Г+ж.

Упражнения
71
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
Преобразовать приближенно данные многочлены в многочлены третьей
степени: а) Р{х) = х5 - Зх4 + 4; б) Р{х) = х4 + 5ж3 - 1. Оценить
допущенные погрешности.
Записать по схеме Горнера алгоритм вычисления первых пяти членов сте-
степенных рядов при разложении функций: а) у = sin х\ б) у = chx.
Используя цепные дроби, вычислить значения: a) In 2; б) tgGr/8);
в) е01; г) arctg0.5.
Дана таблица значений функции
0
1.763
0.2
1.917
0.4
2.143
0.6
2.362
а)	С помощью линейной и квадратичной интерполяций найти приближен-
приближенное значение функции при х = 0.25.
б)	Вычислить, при каком значении аргумента справедливо равенство у =
= 2.000.
Убедиться, что формула линейной интерполяции B.35) полученная как
частный случай многочлена Лагранжа, эквивалентна выведенной ранее
формуле B.28).
Записать алгоритм вычисления функции с помощью квадратичной интер-
интерполяции. Алгоритм должен работать при любом значении аргумента х .
Построить интерполяционный многочлен Лагранжа для функции, заданной
таблицей в упр. 5.
Вычислить значение функции, заданной в упр. 5, при х = 0.1, используя
интерполяционный многочлен Ньютона. Оценить погрешность результата.
Построив график функции шп(х), заданной формулой B.42), проверить,
как изменяется погрешность интерполяции при изменении х.
Найти величину ускорения при равноускоренном движении тела, если из-
известны значения пройденного им пути S в некоторые моменты вре-
времени t:
t, с
S, м
0
5
5
150
10
560
15
1200
20
2100
25
3300
Закон Гука можно записать в виде сг = Ее, где сг — напряжение, Е — мо-
модуль упругости, е — относительная деформация. При испытаниях образца
произвели п измерений значений а и е. Написать алгоритм вычисления
параметра Е.
Изучается зависимость между, электродвижущей силой Е и температу-
температурой нагревания Т термопары. Данные измерений приведены в следующей
таблице:
т, °с
Е, мВ
500
3.23
750
4.52
1000
5.71
1250
10.17
1500
18.49
Найти приближенную зависимость Е(Т) в виде квадратного трехчлена.
14*. Получить соотношения B.73).

ГЛАВА 3
ДИФФЕРЕНЦИРОВАНИЕ
И ИНТЕГРИРОВАНИЕ
§1. Численное дифференцирование
1. Аппроксимация производных. Напомним, что производной функ-
функции у = f(x) называется предел отношения приращения функции Ау к
приращению аргумента Ах при стремлении Ах к нулю:
——, Ay = f(x + Ах) — f(x).	C.1)
Обычно для вычисления производных используют готовые формулы
(таблицу производных) и к выражению C.1) не прибегают. Однако в чис-
численных расчетах на компьютере использование этих формул не всегда удоб-
удобно и возможно. В частности, функция у = f(x) может быть задана в
виде таблицы значений (полученных, например, в результате численного
расчета). В таких случаях производную можно найти, опираясь на форму-
формулу C.1). Значение шага Ах полагают равным некоторому конечному числу,
и для вычисления значения производной получают приближенное ра-
равенство
у' « Ау/Ах.	C.2)
Это соотношение называется аппроксимацией (приближением) произ-
производной с помощью отношения конечных разностей (значения Ау, Ах
в формуле C.2) конечные в отличие от их бесконечно малых значе-
значений в C.1)).
Рассмотрим аппроксимацию производной для функции у = f(x), за-
заданной в табличном виде: у = у о , у\, ... , в узлах х = xq , х\, ... Пусть
шаг — разность между соседними значениями аргумента — постоянный и
равен h. Запишем выражения для производной у[ в узле х = х\, который
слева отмечен крестиком. Используемые при этом узлы (шаблон) отмечены
кружочками. В зависимости от способа вычисления конечных разностей
получаем разные формулы для вычисления производной в одной и той же
точке:
О ®	Ау1=у1-у0, Ax = h, у[& г	C.3)
а
с помощью левых разностей;
® О Ayi =2/2 — 2/i 5 ^х = h-> Hi~ —г—	C.4)
h
с помощью правых разностей;

§ 1. Численное дифференцирование	1Ъ
О х О Аух = 2/2 - 2/о, Ах = 2/г, у[ и 2^	C.5)
с помощью центральных разностей.
Можно найти также выражения для старших производных. Например,
ПбЪП 7/" - (v'Y ~ ^ ~ ^ ~ B/2 - У1)М ~ B/1 ~ 2/о)Л _
U 09 U 2/1 — U/iJ	7	7	 —
2/2 - 22/i + 2/о /о ^ч
Л2	* W'^
Таким образом, используя формулу C.2), можно найти приближенные
значения производных любого порядка. Однако при этом остается откры-
открытым вопрос о точности полученных значений. Кроме того, как будет пока-
показано ниже, для хорошей аппроксимации производной нужно использовать
значения функции во многих узлах, а в формуле C.2) это не предусмотрено.
2. Погрешность численного дифференцирования. Аппроксимируем
функцию f(x) некоторой функцией (р(х), т. е. представим ее в виде
f(x) =(p(x) + R(x).	C.7)
В качестве аппроксимирующей функции (р(х) можно принять частичную
сумму ряда или интерполяционную функцию. Тогда погрешность аппрок-
аппроксимации R(x) определяется остаточным членом ряда или интерполяцион-
интерполяционной формулы.
Аппроксимирующая функция (р(х) может быть использована также
для приближенного вычисления производной функции f(x). Дифференци-
Дифференцируя равенство C.7) необходимое число раз, можно найти значения произ-
производных f'(x), f"(x),...:
В качестве приближенного значения производной порядка к функ-
функции f(x) можно принять значение соответствующей производной функ-
функции ip(x), т. е. f(k\x) и ср(к\х). Величина
характеризующая отклонение приближенного значения производной от
ее истинного значения, называется погрешностью аппроксимации произ-
производной.
При численном дифференцировании функции, заданной в виде таб-
таблицы с шагом h, эта погрешность зависит от h, и ее записывают в виде
()	(о большое от hr) х\ Показатель степени г называется
Напомним, это означает, что \R^ | < Chr, где С > 0 и не зависит от h.

74	Гл. 3. Дифференцирование и интегрирование
порядком погрешности аппроксимации производной (или порядком точ-
точности данной аппроксимации). При этом предполагается, что значение
шага по модулю меньше единицы.
Оценку погрешности легко проиллюстрировать с помощью ряда Тейлора
f(x + Ах) = f(x) + f'(x)Ax + J-^- Ax2 + ^p Ax3 + ...
Пусть функция f(x) задана в виде таблицы f(xi) = y^ (г = 0,1,..., п).
Запишем ряд Тейлора при х = х\, Ах = — h с точностью до членов по-
порядка h2:
Уо =У1 -y[h^O(h2).
Отсюда найдем значение производной в точке х = х\\
Это выражение совпадает с формулой C.3), которая, как видно, является
аппроксимацией первого порядка (г = 1). Аналогично, записывая ряд
Тейлора при Ах = h, можно получить аппроксимацию C.4). Она также
имеет первый порядок.
Используем теперь ряд Тейлора для оценки погрешностей аппроксима-
аппроксимаций C.5) и C.6). Полагая Ах = h и Ах = — h соответственно, получаем
C-8)
Вычитая эти равенства одно из другого, после очевидных преобразова-
преобразований получаем
Это аппроксимация производной C.5) с помощью центральных разностей.
Она имеет второй порядок.
Складывая равенства C.8), находим оценку погрешности аппроксима-
аппроксимации производной второго порядка вида C.6):
Таким образом, эта аппроксимация имеет второй порядок. Аналогично
можно получить аппроксимации производных более высоких порядков и
оценку их погрешностей.
Мы рассмотрели лишь один из источников погрешности численного
дифференцирования — погрешность аппроксимации (ее также называют
погрешностью усечения). Она определяется величиной остаточного члена.

§ 1. Численное дифференцирование	75
Анализ остаточного члена нетривиален, некоторые сведения по этому
вопросу приведены в п. 3. Отметим, лишь, что погрешность аппроксимации
при уменьшении шага /г, как правило, уменьшается.
Погрешности, возникающие при численном дифференцировании, опре-
определяются также неточными значениями функции щ в узлах и погрешностя-
погрешностями округлений при проведении расчетов на компьютере. Обусловленные
этими причинами погрешности, в отличие от погрешности аппроксимации,
возрастают с уменьшением шага h. Действительно, если при вычислении
значений функции у = f(x) абсолютная погрешность составляет d, то
при вычислении дробей в C.3) и C.4) она составит 2d/h. Поэтому сум-
суммарная погрешность численного дифференцирования может убывать при
уменьшении шага лишь до некоторого предельного значения, после чего
дальнейшее уменьшение шага не повысит точности результатов.
Потеря точности аппроксимации производных может быть предотвра-
предотвращена за счет регуляризации процедуры численного дифференцирования.
Простейшим способом регуляризации является такой выбор шага h, при
котором справедливо неравенство \f(x + К) — f(x) \ > е, где е > 0 — неко-
некоторое малое число. При вычислении производной это исключает вычитание
очень близких по величине чисел, которое обычно приводит к увеличению
погрешности. Это тем более опасно при последующем делении прираще-
приращения функции на малое число h. Другой способ регуляризации заключается в
оценке суммарной погрешности численного дифференцирования и выборе
такого шага h, который минимизировал бы эту суммарную погрешность.
Возможен и еще один подход — сглаживание табличных значений функ-
функции подбором некоторой гладкой аппроксимирующей функции, например,
многочлена.
3. Использование интерполяционных формул. Предположим, что
функция f(x), заданная в виде таблицы с постоянным шагом h = Х{ —
— Xi-i (i = 1, 2,..., n), может быть аппроксимирована интерполяцион-
интерполяционным многочленом Ньютона B.39):
у и N(x0 + th) =yo+ tAy0 + ^ ~ ' А2у0 + ...
ф-1)...(*-п + 1)	Х-Хо
+А2/ t
п\
Дифференцируя этот многочлен по переменной х с учетом правила диф-
дифференцирования сложной функции:
dx dt dx h dt
можно получить формулы для вычисления производных любого порядка:

76
Гл. 3. Дифференцирование и интегрирование
4i3 - 18i2 + 22i - б . 4
+	^	A42/o +
— 6
•24
5!
20i3 - 120i2 + 210* - 100
5!
. + ••• ,
+ ••• ,
Число слагаемых в этих формулах зависит от количества узлов, исполь-
используемых для вычисления производных. Как и при построении многочлена
Ньютона, добавление к шаблону нового узла означает добавление к сумме
одного слагаемого.
Пример. Вычислить в точке х = 0.1 первую и вторую производные
функции, заданной табл. 3.1 .
Таблица 3.1
X
0
0.1
0.2
0.3
0.4
0.5
У
1.2833
1.8107
2.3606
2.9577
3.5969
4.2833
Ау
0.5274
0.5599
0.5971
0.6392
0.6864
А2у
0.0325
0.0372
0.0421
0.0472
А3у
0.0047
0.0049
0.0051
А4у
0.0002
0.0002
А5у
0.0000
Здесь h = 0.1, t = @.1 - 0)/0.1 = 1. Заполняя табл. 3.1 аналогич-
аналогично табл. 2.1 и используя полученные выше формулы, находим
/
у' и 10 ( 0.5274 +
у" и 100 (о.О325
2-1 — 1 Я • 1 — 6-1
	0.0325 +
• 0.0047 +
4-1-18-1 +22-1-6
• 0.0047
24
• 0.0002 | = 3.25.

§ 1. Численное дифференцирование	11
Интерполяционные многочлены Ньютона (а также Стирлинга и Бессе-
Бесселя) дают выражения для производных через разности Аку (к = 1,2,...).
Однако на практике часто выгоднее выражать значения производных не
через разности, а непосредственно через значения функции в узлах. Для по-
получения таких формул удобно воспользоваться формулой Лагранжа с рав-
равномерным расположением узлов (xi - xi-i = h = const, i = 1, 2,..., n).
Запишем интерполяционный многочлен Лагранжа L(x) и его остаточ-
остаточный член Rl(x) (cm. B.34), B.41)) для случая трех узлов интерполяции
(п = 2) и найдем их производные:
L{x) =
1
- хо)(х - x2)yi + (х- хо)(х - х1)у2],
у'"
= -^-(x - хо)(х -
L'(x) = — [Bх ~Х1~ х2)Уо - 2Bж - х0 - х2)у\ + Bх - х0 - х1)у2\,
У1"
^ [(Х X)(X Х) + (Х Х)(Х - Х2) + (Х ~ ХО)(Х - Хг)} .
Здесь yl — значение производной третьего порядка в некоторой внутрен-
внутренней точке ж* G [жо, хп].
Запишем выражение для производной у'о при х = xq :
'0 = Lf(x0)+RfL(x0) =
[B
- х0 - х2)ух + Bх0 - х0 -
-хо)(хо -х2) + (х0 -хо)(хо -х1)} =
C + 4) + '"
Аналогичные соотношения можно получить и для значений у[ и у2 при
х = xi, x2:
У[ = ^(У2 ~ 2/о) - у?/*', 2/2 = 2ftB/о - *У1 + З2/2) + -^у'"-
Для каждой из этих формул значения у'", вообще говоря, различны.
Записывая интерполяционный многочлен Лагранжа и его остаточный
член для случая четырех узлов (п = 3), получаем следующие аппроксима-
аппроксимации производных:

78	Гл. 3. Дифференцирование и интегрирование
2/о = ^("П2/о + 182/1 - 9^/2 + 2у3) - ^ у™,
1	h3
2/i = ^(2/о " 32/1 + б2/2 - 2/з) + у^ 2/*\
1	/г3
2/2 = ^B/о - 62/i + 32/2 + 2у)	V
1	/г3
2/з = ^(2/о + 92/1 " 182/2 + П2/з) + у 2/*\
В случае пяти узлов (п = 4) получим
1	/г4
2/о = y^(-252/o + 482/i - З62/2 + I62/3 - Зу4) + у 2/У,
1	/г4
2/i = Y^(2/o - Ю2/1 + I82/2 - 62/3 + 2/4) - go 2/У,
2/2 = Y2h^yo ~ Syi + 8?/3 ~ У^ + 30
1	/г4
2/з = 12Л (/0 + 62/1 - 182/2 + Щз + З2/4) + ^ 2/У,
1	/г4
2/4 = 77^(ЗУо ~ 162/1 + 362/2 - 482/з + 252/4) + — 2/У •
12/г	5
Таким образом, используя значения функции в п+1 узле, получаем аппрок-
аппроксимацию производных п-го порядка точности. Эти формулы можно исполь-
использовать Не ТОЛЬКО ДЛЯ уЗЛОВ X = Xq, Х\, ... , НО И ДЛЯ ЛЮбыХ уЗЛОВ X = Х{,
Жг+ь • • • •> соответствующим образом изменяя значения индексов.
Обратим внимание на то, что при четных п наиболее простые выра-
выражения и наименьшие коэффициенты в остаточных членах получаются для
производных в средних (центральных) узлах (у[ при п = 2, у'2 при п = 4
и т. д.). Выпишем аппроксимации производных для узла с произвольным
номером г, считая его центральным:
У'г = ^(Уг+1 ~ Уг-l) ~ у 2/1",	П = 2,	C.9)
1	h4
У'г = ГГ-г(Уг-2 ~ Syi-г + 8^+1 - yi+2) + ^ 2/У,	П = 4.
Они называются аппроксимациями производных с помощью центральных
разностей и широко используются на практике. Аппроксимация C.9) —
это не что иное, как уже встречавшаяся нам аппроксимация C.5).
С помощью интерполяционных многочленов Лагранжа можно получить
аппроксимации для старших производных. Приведем аппроксимации для
вторых производных.

§ 1. Численное дифференцирование	79
В случае трех узлов интерполяции (п = 2) имеем
B+)
р(»о -2»1 + й) + 0{h).
В случае четырех узлов (п = 3) имеем
B 5 + 4
УЗ = Д2 ('У» + 42/1 " %2 + 2УЗ) + O(h2).
В случае пяти узлов (п — 4) имеем
У" = ^2"(П2/о - 2°2/1 + 62/2 + 42/з - 2/4)
2/2 = ^(-2/0 + 162/1 - 30J/2 + 1б2/з - 2/4) + 0{h4),
2/з = Ju?(~yo + 42/i + 62/2 - 202/3
2/4 = y^2 (Иг/о " 562/1 + П42/2 - Ю42/3 + 35щ) +
Аппроксимации вторых производных с помощью центральных разнос-
разностей при четных п также наиболее выгодны.
4. Метод неопределенных коэффициентов. Аналогичные формулы
можно получить и для случая произвольного расположения узлов. Исполь-
Использование многочлена Лагранжа в этом случае приводит к вычислению гро-
громоздких выражений, поэтому удобнее применять метод неопределенных
коэффициентов. Он заключается в следующем. Искомое выражение для
производной k-то порядка в некоторой точке х = Х{ представляется в виде
линейной комбинации заданных значений функции в узлах xq, х\, ... , хп:
у^к) и со2/о + ci2/i + ... + спуп.	C.10)
Предполагается, что это соотношение выполняется точно, если функция у
является многочленом степени не выше п, т. е. может быть представлена

80	Гл. 3. Дифференцирование и интегрирование
в виде
у = bo + bi(x - xq) + ... + bn(x - xo)n.
Отсюда следует, что соотношение C.10), в частности, должно выполняться
точно для многочленов у = 1, у = х - х0, ... , у = (х - хо)п. Подставляя
последовательно эти выражения в C.10) и требуя выполнения точного ра-
равенства, получаем систему п + 1 линейных алгебраических уравнений для
определения неизвестных коэффициентов cq, ci, ... , сп.
Пример. Найти выражение для производной у[ в случае четырех
равноотстоящих узлов (п = 3).
Приближение C.10) запишется в виде
Уг ~ соУо + сил + С22/2 + с3уз-	C.11)
Используем следующие многочлены:
2/ = 1, y = x-xQ, у = (х-х0J, у = (х-х0K.	C.12)
Вычислим их производные:
у' = 0, у' = 1, у' = 2(х-х0), у'=3(х-х0J.	C.13)
Подставляем последовательно соотношения C.12) и C.13) соответственно
в правую и левую части C.11) при х = х\, требуя выполнения точного
равенства:
2(Я1
3(a?i-
0 = со • 1 + С]
1 = со(жо - х
- хо) = со(хо - х
- Хо) = Co(xq — X
L -1 +
'о) + с
^оJ +
^оK +
с2 • 1 + с2
ci(xi — х
ci(xi — х
1-1,
) + с2(ж2 -хо
оJ +с2(ж2 -с
оK +с2(ж2 -:
) + с3(х3 -хо
соJ + с3(ж3 -
соK + с3(ж3 -
Получаем окончательно систему уравнений в виде
со + с\ + с2 + с3 = 0,
hc\ + 2hc2 + З/1С3 = 1,
hc\ + 4hc2 + 9/гсз = 2,
ftci + 8/гс2 + 27/гсз = 3.
Решая эту систему, получаем
1	111
СО = "ЗХ' Cl = ~2h> C2 = h> C3 = ~6h'
Подставляя эти значения в C.11), находим выражение для производной:
у[ « ^(~22/о - 3?/1 + б2/2 - 2/з)¦

§ 1. Численное дифференцирование	81
5. Улучшение аппроксимации. Как видно из конечно-разностных со-
соотношений для аппроксимаций производных (см. п. 3), порядок их точности
возрастает с увеличением числа узлов, используемых при аппроксимации.
Однако при большом числе узлов эти соотношения становятся весьма гро-
громоздкими, что приводит к существенному возрастанию объема вычисле-
вычислений. Усложняется также оценка точности получаемых результатов. Вместе
с тем существует простой и эффективный способ уточнения решения при
фиксированном числе узлов, используемых в аппроксимирующих конечно-
разностных соотношениях. Это метод Рунге—Ромберга. Изложим вкрат-
вкратце его сущность.
Пусть F(x) — производная, которая подлежит аппроксимации; /(ж, h) —
конечно-разностная аппроксимация этой производной на равномерной сетке
с шагом h; R — погрешность (остаточный член) аппроксимации, главный
член которой можно записать в виде hp^p{x), ъ е.
Тогда выражение для аппроксимации производной в общем случае можно
представить в виде
F(x) = /(я, h) + hp<p(x) + O(hp+1).	C.14)
Запишем это соотношение в той же точке х при другом шаге h\ = kh.
Получим
F(x) = f(x, kh) + (kh)p(p(x) + O((kh)p+1).	C.15)
Приравнивая правые части равенств C.14) и C.15), находим выражение для
главного члена погрешности аппроксимации производной:
Подставляя найденное выражение в равенство C.14), получаем формулу
Рунге
F(x) = f(x, h) + /(*.*W(*.M) + 0(wh-i}.	(ЗЛ6)
Эта формула позволяет по результатам двух расчетов значений производ-
производной /(ж, К) и /(ж, kh) (с шагами h и kh) с порядком точности р найти ее
уточненное значение с порядком точности р + 1.
Пример. Вычислить производную функции у = х3 в точке х = 1.
Очевидно, что у' = Зх2; поэтому у'A) = 3. Найдем теперь эту произ-
производную численно. Составим таблицу значений функции:
6 Л.И. Турчак, П.В. Плотников
X
У
0.8
0.512
0.9
0.729
1.0
1.0

82
Гл. 3. Дифференцирование и интегрирование
Воспользуемся аппроксимацией производной с помощью левых разностей,
имеющей первый порядок (р = 1). Примем шаг равным 0.1 и 0.2, т. е.
к = 2. Получим
f(xh) -У A0 1) M-W-9) ^O-729 2 71
0.1
0.1
По формуле Рунге найдем уточненное значение производной:
о 71 — 9 44
F(x) = z/(l) « 2.71 + 31 _ ^ = 2.98.
Таким образом, формула Рунге дает более точное значение производной. В
общем случае порядок точности аппроксимации увеличивается на единицу.
Мы рассмотрели уточнение решения, полученного при двух значениях
шага. Предположим теперь, что расчеты могут быть проведены с шага-
шагами h 1, /12, • • • , hq. Тогда можно получить уточненное решение для произ-
производной F(x) по формуле Ромберга, которая имеет вид
F(x) =
... ^+9-2
f(x,hq
... Л5+«
1 ft? hp+1
1 Ы
,p hp+1
2 
C.17)
Таким образом, порядок точности возрастает на q - 1. Заметим, что для
успешного применения уточнения исходная функция должна иметь непре-
непрерывные производные достаточно высокого порядка.
6. Частные производные. Рассмотрим функцию двух переменных и =
= /(ж, у), заданную в табличном виде: щ^ = f(xi,yj), где xi = xq + ih\
(i = 0,1,..., /), yj = уо + jh2 (j = 0,1,..., J). В табл. 3.2 представле-
представлена часть данных, которые нам в дальнейшем понадобятся.
Используя понятие частной производной, можем приближенно записать
для малых значений шагов hi и h2
ди
дх
- f(x,y)
hi
ди
ду
- f(x,y)
Л2
Воспользовавшись введенными выше обозначениями, получим следую-
следующие приближенные выражения (аппроксимации) для частных производных

§ 1. Численное дифференцирование
83
Таблица 3.2
^\ X
У 3-2
Уз-г
Уз
Уз+i
У 3+2
Xi-2
Ui-2,j-2
Щ-2,з-1
Ui-2,j
Ui-2,j + l
Ui-2,j+2
Xi-1
Ui-l,j-l
ui-h3
Щ-1,з+1
Ui-l,j + 2
Xi
Ui,j-2
Uij
Ui,j + 1
Ui,j + 2
Xi+1
Ui+l,j-2
Ui+l,j-l
Ui+l,j
Ui+l,j + l
Ui+l,j+2
Xi+2
Ui+2,j-2
Ui+2,j-l
ui+2,j
Ui+2,j + l
Ui+2,j+2
в узле (xi, г/j) с помощью отношений конечных разностей:
- Ui
fcl
ди
Для численного дифференцирования функций многих переменных
можно, как и ранее, использовать интерполяционные многочлены. Од-
Однако рассмотрим здесь другой способ — разложение в ряд Тейлора функ-
функции двух переменных:
f(x + Дж, у + Ау) =
дх	ду
Ах2 + 2^ Ах Ay+^,
Используем эту формулу дважды:
1)найдем ttj+ij = /(ж, + h\,yj) при Дж =
2)найдем Щ-ij = f(xi - h\,yj) при Дж =
Получим
ди
д2и
0
/н, Ду = 0;
-hi, Ay = 0.
1 {д3и\ ,
Вычитая почленно из первого равенства второе, получаем
C.18)

84	Гл. 3. Дифференцирование и интегрирование
Отсюда найдем аппроксимацию производной с помощью центральных раз-
разностей:
Она имеет второй порядок.
Аналогично могут быть получены аппроксимации производной ди/ду,
а также старших производных. В частности, для второй производной можно
получить
д2и\ _ ui+lij - 2ui:j + щ-^j
)
п1
Записывая разложения в ряд C.18) при разных значениях Ах и Ау,
можно вывести формулы численного дифференцирования с необходимым
порядком аппроксимации.
Приведем окончательные формулы для некоторых аппроксимаций част-
частных производных. Слева указывается используемый шаблон. Значения про-
производных вычисляются в узле (ж^г/j), отмеченном крестиком (напомним,
что на шаблонах и в табл. 3.2 по горизонтали изменяются переменная х и
индекс г, по вертикали — переменная у и индекс j):
00
X
о	,dyjtj	2h2
2
d2u \ ^ Uj+ij+i - Uj+ij-i - Uj-i
oxo
о о

§ 2. Численное интегрирование	85
о
°	((Йи
лл \ду2
U
О
ООО /д2и
О <8> О ( -5^ ) » ^i2 («<+i j+i - 2щ
ООО
OOO
o®o
OOO
Приведенные аппроксимации производных могут быть использованы
при построении разностных схем для решения уравнений с частными про-
производными (см. гл. 8).
§ 2. Численное интегрирование
1. Вводные замечания. Напомним некоторые понятия, необходи-
необходимые для дальнейшего изложения.
Пусть на отрезке [а, Ъ] задана функция у = f(x). С помощью точек
xo,xi,...,xn разобьем отрезок [а, Ь] на п элементарных отрезков [^_i, X{]
(г = 1, 2,... ,п), причем хо = а, хп = Ъ. На каждом из этих отрезков
выберем произвольную точку & (Жг-i ^ ?г ^ ж«) и найдем произведе-
произведение Si значения функции в этой точке /(&) на длину элементарного от-
отрезка Axi = Xi — Xi-i:
8{=/(Ь)Ах{.	C.19)
Составим сумму всех таких произведений:
п
Y^ Axi.	C.20)
г=1
Сумма Sn называется интегральной суммой. Определенным интег-
интегралом от функции f(x) на отрезке [а, Ъ] называется предел интегральной
суммы при таком неограниченном увеличении числа точек разбиения, при
котором длина наибольшего из элементарных отрезков стремится к нулю:
\f(x)dx= lim V/te)Axi.	C.21)

86
Гл. 3. Дифференцирование и интегрирование
Теорема (существования определенного интеграла). Если функ-
функция f(x) непрерывна на [а, Ь], то предел интегральной суммы существует
и не зависит ни от способа разбиения отрезка [а, Ъ] на элементарные
отрезки, ни от выбора точек &.
Геометрический смысл введенных понятий для случая f(x) > 0 про-
проиллюстрирован на рис. 3.1. Абсциссами точек М{ являются значения &,
ординатами — значения /(&). Выражения C.19) при г = 1, 2,..., п опи-
описывают площади элементарных прямоугольников (штриховые линии),
интегральная сумма C.20) — площадь ступенчатой фигуры, образуемой
этими прямоугольниками. При неограниченном увеличении числа точек
деления и стремлении к нулю всех элементов Axi верхняя граница фигуры
(ломаная) переходит в линию у = f(x). Площадь полученной фигуры,
которую называют криволинейной трапецией, равна определенному интег-
интегралу C.21).
Во многих случаях, когда подынтегральная функция задана в аналити-
аналитическом виде, определенный интеграл удается вычислить непосредственно
с помощью неопределенного интеграла (вернее, первообразной) по форму-
формуле Ньютона—Лейбница. Она состоит в том, что определенный интеграл
равен приращению первообразной F(x) на отрезке интегрирования:
ъ
\f(x)dx=F(x)\ba=F(b)-F(a).	C.22)
Однако на практике этой формулой часто нельзя воспользоваться по
двум основным причинам:
У
Мп_х
О хо
X\	X2
Xn-2 Xn-i
Рис. 3.1

§ 2. Численное интегрирование	87
1)	вид функции f(x) не допускает непосредственного интегрирования,
т. е. первообразную нельзя выразить в элементарных функциях;
2)	значения функции f(x) заданы только на фиксированном конечном
множестве точек Xi, т. е. функция задана в виде таблицы.
В этих случаях используются приближенные методы интегрирования.
Они основаны на аппроксимации подынтегральной функции некоторыми
более простыми выражениями, например многочленами.
Одним из таких способов, который может быть использован для вычис-
вычисления интегралов в первом случае, является представление подынтеграль-
подынтегральной функции в виде степенного ряда (ряда Тейлора). Это позволяет свести
вычисление интеграла от сложной функции к интегрированию многочлена,
представляющего первые несколько членов ряда.
1
я
Пример. Вычислить интеграл / = е~ж dx с погрешностью 10~4.
j
о
Воспользуемся разложением экспоненты в ряд:
еж = 1 + ж+ — + — + ...
Используя последнее выражение и заменяя х на -х2, записываем интег-
интеграл в виде
о
= х- — +
х3 хъ	х7
5-2! 7-3!
Более универсальными методами, которые пригодны для обоих слу-
случаев, являются методы численного интегрирования, основанные на ап-
аппроксимации подынтегральной функции с помощью интерполяционных
многочленов. Такая аппроксимация позволяет приближенно заменить
определенный интеграл конечной суммой
оцуи	C.23)
где yi — значения функции в узлах интерполяции, с^ — числовые коэф-
коэффициенты. Соотношение C.23) называется квадратурной формулой, а его
правая часть — квадратурной суммой. В зависимости от способа ее
вычисления получаются разные методы численного интегрирования
(квадратурные формулы) — методы прямоугольников, трапеций, па-
парабол, сплайнов и др.

Гл. 3. Дифференцирование и интегрирование
Квадратурную сумму можно вычислить по аналогии с интегральной
суммой C.20)
г=0	г=1
где di — приближенное значение площади элементарной криволинейной
трапеции, соответствующей элементарному отрезку [а^_х,Жг]. Например,
можно положить &i = Si при некотором выборе точки & в C.19). В даль-
дальнейшем при вычислении квадратурной суммы будем аппроксимировать
подынтегральную функцию, используя кусочную (локальную) интерпо-
интерполяцию.
Следует отметить, что к вычислению определенного интеграла сводят-
сводятся многие практические задачи: вычисление площади фигур, определение
работы переменной силы и т. д. Решение задач с использованием кратных
интегралов также обычно может быть в конечном итоге сведено к вычис-
вычислению определенных интегралов.
2. Методы прямоугольников и трапеций. Простейшим методом чис-
численного интегрирования является метод прямоугольников. Он непосред-
непосредственно использует замену определенного интеграла интегральной сум-
суммой C.20). В качестве точек & могут выбираться левые (& = ж«_х) или
правые (& = Xi) границы элементарных отрезков. Обозначая f(xi) = yi,
Axi = hi, получаем следующие формулы метода прямоугольников соот-
соответственно для этих двух случаев:
ь
f(x) dx и Лх2/о + h2yi + ... + Куп-!,	C.24)
ъ
(х) dx и Л-12/i + h2y2 + ... + hnyn.	C.25)
о
Широко распространенным и более точным является вид формулы пря-
прямоугольников, использующий значения функции в средних точках элемен-
элементарных отрезков (в полуцелых узлах):
C.26)
= хг-1 + ^г/25 2 = 1, 2, . . . , П.
В дальнейшем под методом прямоугольников будем понимать последний
алгоритм (он еще называется методом средних).
В рассмотренных методах прямоугольников используется кусочно пос-
постоянная интерполяция: на каждом элементарном отрезке функция f(x)

§ 2. Численное интегрирование
89
приближается функцией, принимающей постоянные значения (констан-
(константой). При этом площадь всей фигуры (криволинейной трапеции) при-
приближенно складывается из площадей элементарных прямоугольников. На
рис. 3.2 верхняя, средняя и нижняя горизонтальные штриховые линии от-
относятся к элементарным прямоугольникам, которые соответствуют форму-
формулам C.25), C.26) и C.24).
Метод трапеций использует линейную интерполяцию, т. е. график функ-
функции у = f(x) представляется в виде ломаной, соединяющей точ-
точки (xi,yi).B этом случае площадь всей фигуры приближенно складывается
из площадей элементарных прямолинейных
трапеций (рис. 3.2). Площадь каждой такой
трапеции равна произведению полусуммы
оснований на высоту:
У
(Уг =
Уг-1 +Уг
% = 1, 2, . . . , П.
Складывая все эти равенства, получаем фор-
формулу трапеций для численного интегрирова-
интегрирования:
C-27)
О
I I
I I
Уг-1
i
I IK
I
I
Xi-i Xi_1/2Xi X
Рис. 3.2. Вычисление cri в ме-
Важным частным случаем рассмотрен- тодах прямоугольников и тра-
ных формул является их применение при	пеций
численном интегрировании с постоянным
шагом hi = h = const (г = 1,2, ...,п). Формулы прямоугольников и
трапеций в этом случае принимают соответственно вид
r	"
*—1
Г	/	_L	"-1 \
j f(x) dxxihl	+ 2^ Vi I •
C.28)
C.29)
Рассмотрим пример использования этих формул при ручном счете для
простейшего интеграла, допускающего также непосредственное вычисле-
вычисление. Такой пример позволит сравнить результаты расчетов, полученные
различными способами.
1
Г dx
Пример. Вычислить интеграл / = 	^.
J 1 + хЛ
о

90
Гл. 3. Дифференцирование и интегрирование
Этот интеграл легко вычисляется по формуле C.22):
/ = arctgzlj = j и 0.785398.
Используем теперь для вычисления данного интеграла формулы прямо-
прямоугольников и трапеций. Разобьем отрезок интегрирования [0,1] на десять
равных частей: п = 10, h = 0.1. Вычислим значения подынтегральной
функции yi = 1/A + ж?) в точках разбиения xi = a^-i + К а также
в полуцелых точках xi_1/2 = #«-i + h/2 (i = 1, 2,..., 10) (табл. 3.3).
Таблица 3.3
г
0
1
2
3
4
5
6
7
8
9
10
Xi
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Уг
1.000000
0.990099
0.961538
0.917431
0.862069
0.800000
0.735294
0.671141
0.609756
0.552486
0.500000
Жг-1/2
0.05
0.15
0.25
0.35
0.45
0.55
0.65
0.75
0.85
0.95
Уг-1/2
0.997506
0.977995
0.941176
0.890868
0.831601
0.767754
0.702988
0.640000
0.580552
0.525624
По формуле прямоугольников C.28) получим
10
= 0.1 • @.997506 + ... + 0.525624) = 0.785606.
Погрешность в вычислении интеграла составляет А/ = / —1\ = —0.00021
(около 0.027 %). Используя формулу трапеций C.29), находим
h = 0.1 • @.750000 + 0.990099 + ... + 0.552486) = 0.784981.
Погрешность здесь равна А/2 = 0.00042 (около 0.054%).
Таким образом, в рассмотренном примере лучшую точность вычис-
вычисления интеграла дает формула прямоугольников. Это, на первый взгляд,
неожиданный результат, поскольку формула прямоугольников использует
интерполяцию нулевого порядка (кусочно - постоянную), в то время как
формула трапеций использует кусочно - линейную интерполяцию. Повы-
Повышение точности здесь объясняется способом вычисления элементарных

§ 2. Численное интегрирование
91
площадей сг^, использующим значения функции в центральной точке #i_i/2
отрезка [ж^-ъ ?«]• Заметим, что использование формул прямоугольников в
виде C.24) или C.25) приведет к погрешности более 3 %.
В п. 5 показано, что погрешность численного интегрирования определя-
определяется шагом разбиения. Уменьшая этот шаг, можно добиться большей точ-
точности. Правда, увеличивать число точек не всегда возможно. Если функция
задана в табличном виде, приходится, как правило, ограничиваться дан-
данным множеством точек. Повышение точности может быть в этом случае
достигнуто за счет повышения степени используемых интерполяционных
многочленов. Рассмотрим два таких способа численного интегрирования:
использование квадратичной интерполяции (метод Симпсона) и интерпо-
интерполирование с помощью сплайнов.
3. Метод Симпсона. Разобьем отрезок интегрирования [а, Ъ] на чет-
четное число п равных частей с шагом h. Па каж-
каждом отрезке [ж0, х2], [ж2, Х4, •-, [xi-uxi+1],...
... , [хп-2,хп] подынтегральную функцию f(x)
заменим интерполяционным многочленом вто-
второй степени:
f(x)
(fi(x) =
+ biX +
Коэффициенты этих квадратных трехчленов мо-
могут быть найдены из условий равенства много-
члена в точках xi соответствующим табличным
данным yi. В качестве (fi(x) можно принять ин-
интерполяционный многочлен Лагранжа второй степени, проходящий через
точки Mi-^Xi-uyi-!), Mi(xi,yi), Mi+1(xi+1, 2/t+i):
<Рг(х) =
(х -
¦2/t-i
~r
(x — Xi-\
У г
{x - Xi-{){x -
Уг+1-
Сумма элементарных площадей сг« и ai+\ (рис. 3.3) может быть вычис-
вычислена с помощью определенного интеграла. Учитывая равенства Xi+\ —Xi =
= Xi — Xi-i = h, получаем
Жг + 1	Жг + 1
Г	1 Г
=	(fi(x)dx=—	[(ж-Ж<)
Xi-1
- 2{х - Xi-!
yi-\-(x- Xi-!
] dx =
= з (Уг-l

92	Гл. 3. Дифференцирование и интегрирование
Проведя такие вычисления для каждого элементарного отреза [жг-ъ
просуммируем полученные выражения:
S = ~(уо+ 4^/1 + 2у2 + 4у3 + 2у4 + ... + 22/п_2 + 42/гг_1 + уп).
о
Данное выражение для S принимается в качестве значения определенного
интеграла:
ь
Г	h
f(x) dx и - [(уо + 4B/! + 2/з + • • • + 2/n-i) +
+ 2B/2 + 2/4 + • • • + 2/n-2) + уп] • C.30)
Полученное соотношение называется формулой Симпсона или формулой
парабол.
Эту формулу можно получить и другими способами, например двукрат-
двукратным применением метода трапеций при разбиениях отрезка [а, Ъ] на части с
шагами h и 2/г или комбинированием формул прямоугольников и трапеций
(см. п. 5).
Иногда формулу Симпсона записывают с применением полуцелых ин-
индексов. В этом случае число отрезков разбиения п произвольно (не обяза-
обязательно четно), и формула Симпсона имеет вид
ъ
С	h
f(x) dx и - [B/о + 4B/1/2 + 2/3/2 + ¦ ¦ ¦ + Уп-1/2) +
+ 2B/i+2/2+ ... + 2/n-i) + 2/n]. C.31)
Легко видеть, что формула C.31) совпадет с C.30), если формулу C.30)
применить для числа отрезков разбиения 2п и шага h/2.
1
Г dx
Пример. Вычислить по методу Симпсона интеграл / = 	^.
J 1 + х
о
Значения функции при п = 10, h = 0.1 приведены в табл. 3.3.
Применяя формулу C.30), находим
1= -^-[Уо+ %i + 2/з + 2/5 + 2/7 + 2/9) +
+ 2B/2 + 2/4 + 2/6 + 2/8) + 2/ю] = • • • = 0.785398.
Результат численного интегрирования с использованием метода Симп-
Симпсона оказался совпадающим с точным значением (шесть значащих цифр).
Один из возможных алгоритмов вычисления определенного интеграла
по методу Симпсона представлен на рис. 3.4. В качестве исходных данных
задаются границы отрезка интегрирования [а, Ь], погрешность е, а также
формула для вычисления значений подынтегральной функции у = f(x).

§ 2. Численное интегрирование
93
Первоначально отрезок [а, Ь] разби-
разбивается на две части с шагом h =
= (Ъ — а)/2. Вычисляется значение
интеграла 1\. Потом число шагов
удваивается, вычисляется значение
/2 с шагом h/2. Условие окончание
счета принимается в виде |/i — /21 <
< е. Если это условие не выполне-
выполнено, происходит новое деление шага
пополам и т. д.
Отметим, что представленный
на рис. 3.4 алгоритм не являет-
является оптимальным: при вычислении
каждого приближения /2 не исполь-
используются значения функции f(x), уже
найденные на предыдущем этапе.
Более экономичные алгоритмы бу-
будут рассмотрены в п. 6.
4. Использование сплайнов.
Одним из методов численного ин-
интегрирования, особенно эффектив-
эффективным при строго ограниченном чис-
числе узлов, является метод сплай-
сплайнов, использующий интерполяцию
сплайнами (см. гл. 2, § 3, п. 5).
Разобьем отрезок интегрирова-
интегрирования [а, Ь\ на п частей точками
',-КП
ДО
для г от
X -
до п -
Ввод а, Ъ, е
-- 2, h = (Ъ -
h = /2, h
п = 2п, h -
1
= а -\- ih, 1ч-
1 с шагом 2
для г от 2
X -
до п -
h =
- а-\- ih, h -
2 с шагом 2
Uf(a) + 2/2
Вывод /2
а)/2,
= 0,
= ft/2
= /2 +
= /2 +
+ /F)
;ь)]/з
/(*)
]/з
Рис. 3.4. Метод Симпсона
Пусть Axi = hi (i = 1, 2,..., п). На каждом элементарном отрезке интерпо-
интерполируем подынтегральную функцию f(x) с помощью кубического сплайна:
ж«_1 ^ ж ^ Xi, i = 1, 2,..., п.
Выражение для интеграла представим в виде
C.32)
C.33)
Используя выражение C.32), в результате вычисления интегралов находим
1
hi + 2 Ь
4 *ЛП •
C.34)
Способ вычисления коэффициентов ai, Ь{, С{, &{ описан в гл. 2. Здесь
лишь отметим, что а« = Уг-\-

94	Гл. 3. Дифференцирование и интегрирование
Подставив в C.34) выражения для а«, bi, di, эту формулу можно пред-
представить в удобном для практических расчетов виде
i\
Отметим, что во всех предыдущих методах (см. п. 2, 3) формулы чис-
численного интегрирования можно записать в виде линейной комбинации
табличных значений функции, т. е. в виде квадратурной формулы C.23)
с постоянными коэффициентами с^- При использовании сплайнов такое
представление невозможно, поскольку коэффициенты с^ зависят в этом
случае от всех значений yi%
5. Погрешность численного интегрирования. При вычислении
приближенного значения интеграла по квадратурной формуле C.23) до-
допускается погрешность
ь
Г
=
R f()
I	i=0
Она зависит от шага разбиения, и ее можно представить в виде R = O(hr).
В случае переменного шага можно принять h = max hi, hi = Ax{. Как и в
случае численного дифференцирования, показатель степени г называют
порядком точности данной квадратурной формулы (или данного метода).
Квадратурная формула должна быть составлена таким образом, чтобы
для любой интегрируемой на отрезке [а, Ъ] функции f(x) при h —> О
(п —> ею) значения интеграла, получаемые путем численного интегри-
интегрирования, сходились к его точному значению. Это означает выполнение
неравенства г > 0.
Погрешность R, допускаемую при интегрировании функции по отрез-
отрезку [а, Ь], можно представить в виде суммы погрешностей г«, допускаемых
на каждом элементарном отрезке:
п	Хр
, п=\ f(x)dx-(Ji.	C.36)
Получим выражения для погрешностей формул прямоугольников и тра-
трапеций. Запишем разложение функции у = f(x) в ряд Тейлора на отрез-
отрезке [xi-UXi\:
f(x) = Уг-г/2 + у[-1/2(Х - xi-l/2) + у У"-1/2(Х
C.37)
Заметим, что разложение C.37) можно оборвать на k-м слагаемом, ес-
если в этом слагаемом производную 2/г-_1/2 заменить на f(k\x*), где ж* —

§ 2. Численное интегрирование
95
некоторая точка отрезка [xi-i,Xi]. Например,
/0*0 = Уг-1/2 + У'г-1/2(х ~ Xi-\li) + у /"(Ж*)(Ж ~ ^
Проинтегрировав обе части равенства C.37) по отрезку
лучим
г-ъ xi\> по-
поf(x)dx = 2/i_i
Xi-!
Xi
Xi-!
1
24
Для метода прямоугольников yi-\j^hi — &i\ отсюда
f). C.38)
C.39)
Для метода трапеций О{ = (уг-i + yi)hi/2. Вычислим сг^, для чего
найдем ^/г—1 и У г из C.37) (подставив в C.37) х = xi-1 и ж = ж«
соответственно), сложим полученные выражения и домножим их сумму
на hi/2:
Из этого соотношения выразим Уг-1/2Ы и подставим в C.38):
Отсюда
C.40)
Выражение для погрешности метода трапеций C.40) можно получить и
другим способом, проинтегрировав по отрезку [а^-ъ^г] интерполяцион-
интерполяционный многочлен с учетом остаточного члена B.41).
Используя формулы прямоугольников и трапеций, можно получить
уточненные значения интегралов, если учесть характер погрешностей этих
формул. Как следует из C.39) и C.40), главный член погрешности формулы
трапеций вдвое больше по модулю и имеет другой знак. На основании
этого можно записать уточненную формулу для вычисления определенного

96	Гл. 3. Дифференцирование и интегрирование
интеграла с использованием значений /пр и /тр, вычисленных по методам
прямоугольников и трапеций:
I « B/пр + 7тр)/3.	C.41)
Погрешность полученной формулы составит
Нетрудно убедиться, что формула C.41) совпадает с C.31) — формулой
Симпсона, записанной с помощью полуцелых индексов.
Предположим теперь шаг разбиения постоянным, hi = h = const (г =
= 1, 2,..., п), и оценим погрешность метода прямоугольников на отрез-
отрезке [а, Ъ] согласно C.36) и C.39):
fi E i
|| й |! fi
г=1	г=1	г=1
г=1
где
М2 = max \f"(x)
При выводе C.42) мы воспользовались неравенством
а\ +а2 + ... + ап\ ^ |ai| + \а2\ + ... + \ап\.
Аналогично C.42), для метода трапеций имеем
F-а)/*2
' р 	12
Для метода Симпсона C.30) можно получить следующую оценку по-
погрешности:
где М4 — максимум модуля четвертой производной функции f(x).
Сравнив методы прямоугольников и трапеций с методом Симпсона,
отметим, что последний обладает более высоким порядком точности —
четвертым, в то время как методы прямоугольников и трапеций — вторым.
Для анализа погрешности интегрирования с использованием сплайнов
рассмотрим формулу C.35). Первый член в ее правой части совпадает с
правой частью формулы C.27) для метода трапеций. Следовательно, второй
член характеризует поправку к методу трапеций, которую дает использова-
использование сплайнов.
Как следует из формулы C.32), коэффициенты С{ выражаются через
вторые производные (р"(х):

§ 2. Численное интегрирование	97
Это позволяет оценить второй член правой части формулы C.35):
Полученная оценка показывает (см. C.40)), что добавка к формуле тра-
трапеций, которую дает использование сплайнов, компенсирует погреш-
погрешность самой формулы трапеций. Можно показать, что метод сплайнов, как
и метод Симпсона, имеет четвертый порядок точности.
Рассмотрев разные методы численного интегрирования, трудно срав-
сравнивать их достоинства и недостатки. Любая попытка такого сравнения
непременно поставит перед нами альтернативный вопрос: что больше, h2y"
или h4yw? Все зависит от самой функции у = f(x) и поведения ее произ-
производных.
Уточнение результатов численного интегрирования можно проводить
по-разному. В частности, в представленном на рис. 3.4 алгоритме с ис-
использованием метода Симпсона проводится сравнение двух значений ин-
интеграла 1\ и /2, полученных при разбиениях отрезка [а, Ь\ соответственно
с шагами huh/2. Аналогичный алгоритм можно построить и для других
методов.
Здесь мы упомянем другую схему уточнения значения интеграла —
процесс Эйткена. Он дает возможность оценить погрешность O(hr) ме-
метода и указывает алгоритм уточнения результатов. Расчет проводится
последовательно три раза при различных шагах разбиения hi, h2, ^з?
причем их отношения постоянны: h2/hi = h3/h2 = q (например, при
делении шага пополам q = 0.5). Пусть в результате численного интегриро-
интегрирования получены значения интеграла /i, 12, /3- Тогда уточненное значение
интеграла вычисляется по формуле
il — Zl2 -\- Is
а порядок точности используемого метода численного интегрирования оце-
оценивается соотношением
1-/3-/2
г и 	ш	.
In q I2 — Ii
Уточнение значения интеграла можно также проводить методом Рунге -
Ромберга (см. § 5, п. 5).
6. Адаптивные алгоритмы. Из анализа погрешностей методов чис-
численного интегрирования следует, что точность получаемых результатов за-
зависит как от характера изменения подынтегральной функции, так и от шага
интегрирования. Будем считать, что величину шага мы задаем. При этом
ясно, что для достижения сравнимой точности при интегрировании слабо
меняющейся функции шаг можно выбирать большим, чем при интегриро-
интегрировании резко меняющихся функций.
7 Л.И. Турчак, П.В. Плотников

98	Гл. 3. Дифференцирование и интегрирование
На практике нередко встречаются случаи, когда подынтегральная функ-
цияменяется по-разному на отдельных участках отрезка интегрирования.
Это обстоятельство требует такой организации экономичных численных
алгоритмов, при которой они автоматически приспосабливались бы к
характеру изменения функции. Такие алгоритмы называются адаптив-
адаптивными {приспосабливающимися). Они позволяют вводить разные значения
шага интегрирования на отдельных участках отрезка интегрирования.
Это дает возможность уменьшить машинное время без потери точности
результатов расчета. Подчеркнем, что этот подход используется обычно
при задании подынтегральной функции у = f(x) в виде формулы, а не в
табличном виде.
Программа, реализующая адаптивный алгоритм численного интегри-
интегрирования, входит обычно в виде стандартной подпрограммы в математиче-
математическое обеспечение компьютера. Пользователь готовой программы задает
границы отрезка интегрирования а, 6, допустимую абсолютную погреш-
погрешность е и составляет блок программы для вычисления значения подын-
подынтегральной функции f(x). Программа вычисляет значение интеграла / с
заданной погрешностью е, т. е.
f(x) dx —
s.	C.43)
Разумеется, не для всякой функции можно получить результат с за-
заданной погрешностью. Поэтому в программе может быть предусмотре-
предусмотрено сообщение пользователю о недостижимости заданной погрешности.
Интеграл при этом вычисляется с максимально возможной точностью, и
программа выдает эту реальную точность.
Рассмотрим принцип работы адаптивного алгоритма. Первоначаль-
Первоначально отрезок [а, Ъ] разбиваем на п частей. В дальнейшем каждый такой
элементарный отрезок делим последовательно пополам. Окончательное
разбиение отрезка зависит от подынтегральной функции и допустимой
погрешности е.
К каждому элементарному отрезку [ж^-ь^] применяем формулы
численного интегрирования при двух различных его разбиениях. Полу-
Получаем приближения ц и /| для интеграла по этому отрезку:
Xi
1{= Г f(x)
dx.	C.44)
Полученные значения сравниваем и проводим оценку их погрешности.
Если погрешность находится в допустимых границах, то одно из этих
приближений принимается за значение интеграла по данному элементар-
элементарному отрезку. В противном случае происходит дальнейшее деление от-
отрезка и вычисление новых приближений. С целью экономии машинного

§ 2. Численное интегрирование	99
времени точки деления располагаются таким образом, чтобы использова-
использовались вычисленные значения функции в точках предыдущего разбиения.
Например, при вычислении интеграла C.44) по формуле Симпсона от-
отрезок [жг-ъ хг\ сначала разбиваем на две части с шагом hi/2 и вычисляем
значение ц\ Потом вычисляем If* с шагом hi/4, If** с шагом hi/8 и т. д.
Получим выражения
[/() + 4/ (xi-г + |) + /(х,)] ,	C.45)
2/ (xt-! + у) + 4/ (xi-г + ^j+ /(**)] , C.46)
4/ (*«_! + |) + 2/
4/ (*,_! + ^) + 2/ (*«_! + у) + 4/
i of I T. . I l 1 \ Af I т i I	1 I fir \	И 47^
\ z/ I хг-1 \ д /	V г~1	Я / J\Xl) ' W'^'J
Процесс деления отрезка пополам и вычисления уточненных значений
/{ * и /} +1^ (к = 1,2,...) продолжается до тех пор, пока их разность
станет не больше некоторой заданной величины Si, зависящей от е и h:
: Si.	C.48)
Аналогичная процедура проводится для всех п элементарных отрезков.
п
Величина / = J] 1г принимается в качестве искомого значения интеграла.
Условия C.48) и соответствующий выбор величин Si обеспечивают выпол-
выполнение условия C.43).
Замечание. Нетрудно увидеть, что в формулах C.45)-C.47) зна-
значения функции f(x) во внутренних точках отрезка [я^_1,ж*] первый раз
появляются в сумме, заключенной в квадратные скобки, с коэффициентом 4.
В формулах для последующих приближений те же слагаемые входят в сум-
сумму с коэффициентом 2. Этот факт позволяет организовать процесс последо-
последовательного деления отрезка [жг-ъ хг\ пополам так, чтобы при нахождении
приближения l\ +1* значения функции /(ж), вычисленные ранее при на-
нахождении предыдущих приближений, заново не вычислялись. Например,
можно отдельно запоминать сумму значений функции, вычисленных на те-
текущем шаге, сумму значений функции, вычисленных на предыдущем шаге,
а также сумму прочих слагаемых, стоящих в квадратных скобках.

100	Гл. 3. Дифференцирование и интегрирование
7. О других методах. Особые случаи. Кроме рассмотренных выше ме-
методов численного интегрирования существует ряд других. Дадим краткий
обзор некоторых из них.
Формулы Ньютона - Котеса получаются путем замены подынтеграль-
подынтегральной функции интерполяционным многочленом Лагранжа с разбиением
отрезка интегрирования на п равных частей. Получающиеся формулы ис-
используют значения подынтегральной функции в узлах интерполяции и яв-
являются точными для всех многочленов некоторой степени х\ зависящей от
числа узлов. Точность формул растет с увеличением степени интерполяци-
интерполяционного многочлена 2).
Метод Гаусса не предполагает разбиения отрезка интегрирования на
равные промежутки. Формулы численного интегрирования интерполяци-
интерполяционного типа ищутся такими, чтобы они обладали наивысшим порядком
точности при заданном числе узлов. Узлы и коэффициенты формул числен-
численного интегрирования находятся из условий обращения в нуль их остаточных
членов для всех многочленов максимально высокой степени.
Формула Эрмита, являющаяся частным случаем формул Гаусса, исполь-
использует многочлены Чебышева для вычисления интегралов вида
Получающаяся формула характерна тем, что все коэффициенты щ в C.23)
равны.
Метод Маркова состоит в том, что при выводе формул Гаусса вводятся
дополнительные предположения о совпадении точек разбиения отрезка по
крайней мере с одним из его концов.
Формула Чебышева представляет интеграл в виде
C.49)
-1
г=0
При этом решается следующая задача:
найти точки хо, xi, ... , хп и коэффициент к такие, при которых фор-
формула C.49) точна для всех многочленов как можно большей степени.
Формула Эйлера использует не только значения подынтегральной функ-
функции в точках разбиения, но и значения ее производных до некоторого по-
порядка на границах отрезка.
*' Это означает, что если подынтегральная функция — многочлен, то квадратур-
квадратурная формула дает точное значение интеграла (естественно, без учета погрешностей
округления).
2) Заметим, что формулы прямоугольников, трапеций и Симпсона являются част-
частными случаями формул Ньютона - Котеса.

§ 2. Численное интегрирование	101
Рассмотрим особые случаи численного интегрирования'.
а)	подынтегральная функция разрывна на отрезке интегрирования;
б)	несобственные интегралы.
а) В ряде случаев подынтегральная функция f(x) или ее производные
в некоторых внутренних точках с& (к = 1,2,...) отрезка интегрирова-
интегрирования [а, Ъ] терпят разрыв. В этом случае интеграл вычисляют численно для
каждого участка непрерывности и результаты складывают. Например, в
случае одной точки разрыва х = с (а < с < Ь) имеем
с	о
f(x) dx = f(x) dx = f(x) dx.
Для вычисления каждого из стоящих в правой части интегралов можно
использовать рассмотренные выше методы.
б) Не так просто обстоит дело с вычислением несобственных интегра-
интегралов. Напомним, что к такому типу относятся интегралы, которые имеют
хотя бы одну бесконечную границу интегрирования или подынтегральную
функцию, обращающуюся в бесконечность хотя бы в одной точке отрезка
интегрирования.
Рассмотрим сначала интеграл с бесконечной границей интегрирования,
например интеграл вида
+ ОО
f(x)dx, 0 < а < +оо.
а
Существует несколько приемов вычисления таких интегралов.
Можно попытаться ввести замену переменных х = а/A — i), которая
превращает промежуток интегрирования [0,+оо) в отрезок [0,1]. При
этом подынтегральная функция и первые ее производные до некоторого
порядка должны оставаться ограниченными.
Еще один прием состоит в том, что бесконечная граница заменяется
некоторым достаточно большим числом А так, чтобы принятое значение
интеграла отличалось от исходного на некоторый малый остаток, т. е.
+ос	А	+ос
f f(x) dx = f f(x) dx + R, R= f f(x)dx.	C.50)
a	a	A
Если функция обращается в бесконечность в некоторой точке х = с
конечного отрезка интегрирования, то можно попытаться выделить особен-
особенность, представив подынтегральную функцию в виде суммы двух функций:
f(x) = (р(х)-\-ф(х). При этом (р(х) ограничена, а ф(х) имеет особен-
особенность в данной точке, но интеграл (несобственный) от нее может быть
вычислен аналитически. Тогда численный метод используется только для
интегрирования ограниченной функции (р(х). Можно также использовать

102	Гл. 3. Дифференцирование и интегрирование
представление интеграла, аналогичное C.50):
Ъ	С\	Ъ	С	С2
f(x) dx = f(x) dx + f(x) dx + R, R = f(x) dx + f(x) dx.
a	a	c\	c\	с
Здесь с\ и С2 — некоторые близкие к с числа.
Еще один вид несобственных интегралов (сингулярные интегралы),
имеющий важное прикладное значение, будет рассмотрен в дальнейшем в
разделе, посвященном сингулярным интегральным уравнениям (см, гл, 9, § 3),
8. Кратные интегралы. Численные методы используются также для
вычисления кратных интегралов. Ограничимся здесь рассмотрением двойных
интегралов вида
jj f(x, y)dxdy.	C.51)
G
Одним из простейших способов вычисления этого интеграла является ме-
метод ячеек.. Рассмотрим сначала случай, когда областью интегрирования G
является прямоугольник: а ^ х ^ 6, с ^ у ^ d. По теореме о среднем найдем
среднее значение функции f(x,y):
(x,y)dxdy, S = (b-a)(d-c).	C.52)
G
Будем считать, что среднее значение приближенно равно значению функции
в центре прямоугольника, т. е, f(x,y) = f(x,y). Тогда из C.52) получим
выражение для приближенного вычисления двойного интеграла:
я
& Sf(x,y),
g	C.53)
а-\-Ь _ с + d
Х= 2 ' У= 2 "
Точность этой формулы можно повысить, если разбить область G на
прямоугольные ячейки AG«j (рис, 3,5): xi-i ^ х ^ Xi (г = 1,2, ...,М),
Уз-1 ^ 2/ ^ 2/j (j = 15 2,..., TV), Применяя к каждой ячейке формулу C.53),
получаем
f(x,y)dxdyttf(xi,yj)AxiAyj.
AGi:j
Суммируя эти выражения по всем ячейкам, находим значение двойного инте-
интеграла:
М N
(ж, у) dx dy и S2 У2 /(^ь Уз)АхгАу3-.	C.54)
JJ-

§ 2. Численное интегрирование
103
У >
d
Уз
Уз-i
с
AVj{
G
\xi У)
чд^Г
1—*¦
Xi-l
Рис. 3.5
b x
В правой части стоит интегральная сумма; поэтому при неограниченном
уменьшении периметров ячеек (или стягивании их в точки) эта сумма стре-
стремится к значению интеграла для любой непрерывной функции f(x,y).
Можно показать, что погрешность такого приближения интеграла для од-
одной ячейки оценивается соотношением
24
b — а
М
f"
J XX
- С
f"
Суммируя эти выражения по всем ячейкам и считая все их площади одинако-
одинаковыми, получаем оценку погрешности метода ячеек в виде
R = ОA/М2 + 1/7V2) = О(Ах2 + Ау2).
Таким образом, формула C.54) имеет второй порядок точности. Для умень-
уменьшения погрешности можно использовать обычные методы сгущения уз-
узлов сетки.. При этом по каждой переменной шаги уменьшают в одинаковое
число раз, т, е, отношение M/N остается постоянным.
Если область G непрямоугольная, то в ряде случаев ее целесообразно
привести к прямоугольному виду путем соответствующей замены перемен-
переменных.. Например, пусть область задана в виде криволинейного четырехуголь-
четырехугольника: а ^ х ^ 6, (fi(x) ^ у ^ ч>2{х). Данную область можно привести к
прямоугольному виду с помощью замены
t =
(f2(x) -
Кроме того, формула C.54) может быть обобщена и на случай более слож-
сложных областей.

104	Гл. 3. Дифференцирование и интегрирование
Другим довольно распространенным методом вычисления кратных ин-
интегралов является их сведение к последовательному вычислению опреде-
определенных интегралов.. Интеграл C.51) для прямоугольной области можно
записать в виде
ь	d
/(ж, у) dxdy = F{x) dx, F{x) = f(x, у) dy.
G	а	с
Для вычисления обоих определенных интегралов могут быть использованы
рассмотренные ранее численные методы.
Если область G имеет более сложную структуру, то она либо приводится
к прямоугольному виду с помощью замены переменных, либо разбивается на
простые элементы.
Для вычисления кратных интегралов используется также метод замены
подынтегральной функции многомерным интерполяционным многочленом..
Вычисление коэффициентов этих многочленов для простых областей обыч-
обычно не вызывает затруднений.
Существует ряд других численных методов вычисления кратных инте-
интегралов. Среди них особое место занимает метод статистических испытаний,
который мы вкратце изложим,
9. Метод Монте-Карло. Во многих задачах исходные данные носят
случайный характер, поэтому для их решения должен применяться статисти-
ко-вероятностный подход. На основе таких подходов построен ряд численных
методов, которые учитывают случайный характер вычисляемых или изме-
измеряемых величин, К ним принадлежит и метод статистических испыта-
испытаний, называемый также методом Монте-Карло 1\ который применяется к
решению некоторых задач вычислительной математики, в том числе и для
вычисления интегралов.
Метод Монте-Карло состоит в том, что рассматривается некоторая слу-
случайная величина ?, математическое ожидание которой равно искомой вели-
величине х:
Проводится серия п независимых испытаний, в результате которых полу-
получается (генерируется) последовательность п случайных чисел ?ь ?2, • • •
... ,?п (выборка), имеющих то же распределение, чтои ?, и по совокуп-
совокупности этих значений находится выборочное среднее ?, которое является
статистической оценкой М?. Искомая величина х полагается приближен-
приближенно равной этой оценке
?(?+6 + +?)
1 ' Название метода произошло от названия города Монте-Карло, знаменитого
своими казино, в которых играют в рулетку, являющуюся одним из простейших
генераторов случайных чисел.

§ 2. Численное интегрирование	105
Пусть г] — равномерно распределенная па отрезке [0,1] случайная
величина. Это означает, что ее плотность распределения задается соотно-
соотношением
Тогда любая функция ? = f(rj) также будет случайной величиной, и ее
математическое ожидание равно
+ ОО	1
г	г
= f(x)pr](x)dx= \f(x)dx.
j	j
Следовательно, читая это равенство в обратном порядке, приходим к выво-
1
ду, что интеграл f(x) dx может быть вычислен как оценка математичес-
о
кого ожидания некоторой случайной величины ?, которая является функ-
функцией случайной величины г] с равномерным законом распределения,
причем оценка М? определяется независимыми реализациями щ случайной
величины г]:
|	1-
Аналогично могут быть вычислены и кратные интегралы. Для двойного
интеграла получим
1
я-
G
где G — квадрат 0 ^ ж ^ 1, 0 ^ г/ ^ 1; rji, Ci — независимые реализа-
реализации случайных величин г], ?, равномерно распределенных на отрезке [0,1],
Для использования метода Монте-Карло при вычислении определенных
интегралов, как и в других его приложениях, необходимо вырабатывать
последовательности случайных чисел с заданным законом распределения.
Существуют различные способы генерирования таких чисел.
Можно построить некоторый физический процесс (генератор) для вы-
выработки случайных величин, однако при использовании компьютера этот
способ не применяется, поскольку, во-первых, трудно дважды получить
одинаковые совокупности случайных чисел, которые необходимы при от-
отладке программ, а во-вторых, такой физический генератор существенно
усложнил бы конструкцию компьютера.
Известны многие таблицы случайных чисел, которые вычислялись незави-
независимо. Их можно ввести в компьютер и при необходимости обращаться к ним.

106
Гл. 3. Дифференцирование и интегрирование
В настоящее время наиболее распространенный способ выработки
случайных чисел на компьютере состоит в том, что в памяти хранится
некоторый алгоритм получения таких чисел по мере потребности в них
(подобно тому как вычисляются значения элементарных функций, а не
хранятся их таблицы). Поскольку эти числа генерируются по наперед
заданному алгоритму, то они не совсем случайны (псевдослучайны), хотя
и обладают свойственными случайным числам статистическими характе-
характеристиками. В современных языках программирования такие алгоритмы
реализованы в виде подпрограмм — датчиков случайных чисел.
Упражнения
1. Функция у = f(x) задана в табличной форме:
X
У
0
1.24
0.2
1.03
0.4
1.36
0.6
1.85
0.8
2.43
1.0
3.14
Вычислить:
а)	значения производной в точках х = 0, 0.4, 1.0 с первым и вторым порядками
точности;
б)	вторую производную в этих же точках со вторым в третьим порядками точ-
точности.
2. Записать алгоритм вычисления производной функции, заданной таблицей с по-
постоянным шагом на некотором отрезке.
Т. Получить формулу Ромберга C.17).
Указание. Воспользоваться представлением погрешности в виде
R = v()	\()
4.	Для функции f(x,y) = sin (x + у2) вычислить все частные производные до вто-
второго порядка включительно в точке @,0), используя различные аппроксимации,
приведенные в п. 6 § 1. Принять hi = hi = 0.1. Сравнить полученные результаты
с точными значениями производных.
1
5.	Вычислить ех dx, используя методы прямоугольников, трапеций и Симпсона.
о
Отрезок интегрирования разделить на десять равных частей.
6.	Используя процесс Эйткена и метод трапеций, вычислить
2
— sin — dx,
J x 2
Чис-
i
ло шагов интегрирования принять равным 4, 8, 46.
7.	Записать алгоритм решения упр. 9.
8.	Записать адаптивный алгоритм вычисления интеграла по методу Симпсона, вос-
воспользовавшись замечанием, сделанным в п. 6 § 2.
9.	С помощью метода Монте-Карло вычислить площадь фигуры, заданной
/~	/^	05 106
уравнением у/х? + у/у^ = 1. Принять п равным 104, 105, 106. Сравнить
ответы с точным значением площади.