Как найти погрешность дифференцированием

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

(5.16)

и

, (5.17)

а
для центральной

. (5.18)

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

. (5.19)

Здесь
и ниже

и

– некоторые точки, расположенные на
интервалах

и

соответственно. Откуда погрешности
этих методов будут иметь вид

, (5.20)

а
оценка абсолютных погрешностей будет
удовлетворять неравенству , (5.21)

где . (5.22)

Таким
образом, формулы (5.2) и (5.3) вычисления
правой и левой разностных производных
имеют первый порядок точности по
.

Для
центральной разностной производной
соответствующие разложения функций

в ряд Тейлора должны учитывать и
производную третьего порядка (вторая
производная при вычитании исчезает)

. (5.23)

Отсюда
получим

(5.24)

и
для оценки абсолютной погрешности будет
справедливо неравенство

(5.25)

где

(5.26)

Таким
образом, производная

вычисляется при помощи формул центральной
разностной производной со вторым
порядком точности по
,
т.е. точнее, чем по формулам (5.2) и (5.3).

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

Оценим
совместное влияния погрешностей
вычисления и метода для вычисления
первой производной. Пусть значение
производной вычисляется по формуле
(5.3), тогда погрешность метода можно
оценить по соотношению (5.21). Если значения
функции

известны с некоторой погрешностью

(),
то погрешность вычисления

будет содержать дополнительное слагаемое

,(5.27)

Пренебрегая
для простоты погрешностью округления,
имеем оценку погрешности в следующем
виде

.(5.28)

Из
формулы (5.28) очевидно, что уменьшение

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

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

будет равно

.(5.29)

Отметим,
что повышение точности метода лишь
отчасти повышает точность вычисления
производной.

Рисунок
5.1. Зависимость погрешности вычисления
первой производной от величины шага.

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

,(5.29)

где

и

– нижний и верхний пределы интегрирования,
а функция

– непрерывная на отрезке
.

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

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

,(5.30)

где

– погрешность вычисления.

Используемые
на практике методы численного
интегрирования можно условно сгруппировать
в зависимости от способа аппроксимации
подынтегральной функции.

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

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

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

4. Методы,
   основанные на случайном выборе узлов,
   где определяются значения интегрируемой
   функции. Эти методы (методы Монте-Карло)
   используются для вычисления кратных
   интегралов, для которых они наиболее
   эффективны.

В
дальнейшем ограничимся методами первой
и второй группы, а также рассмотрим
алгоритмы метода Монте-Карло для кратных
интегралов.

Наиболее
широко на практике для вычисления
определенных интегралов используются
квадратурные формулы – приближенные
равенства вида

, (5.31)

что
соответствует разбиению полного
интеграла на сумму интегралов от функции
по элементарным отрезкам. Здесь

– некоторые точки (узлы) из отрезка
,
разбитого на

элементарных отрезков
,
при этом
,
а
,
а

– некоторые коэффициенты. Сумма, стоящая
в правой части выражения (5.31) называется
квадратурной суммой. Погрешность
вычисления по этой формуле определяется
выражением

. (5.32)

Заметим,
что интеграл вида (5.29) определяет площадь
фигуры ограниченной прямыми
,

,

и графиком функции
.
Для вычисления этой и площади и
соответственно интеграла ее можно
разбить на простые фигуры. На этом и
основаны простейшие методы вычисления
интеграла (5.29).

Формулы
прямоугольников.
Рассмотрим разбиение отрезка

на

элементарных отрезков
,
пусть для простоты
,
тогда каждому значению

можно поставить в соответствие значение
функции
.
Рассматривая

как высоту прямоугольника с основанием

можно получить две формулы: формулу
левых прямоугольников (рисунок 5.2)

(5.33)

и
формулу правых прямоугольников

.(5.34)

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

отрезков и рассматривают только нечетные
узлы, что эквивалентно рассмотрению
прямоугольников высотой, равной значению
функции в середине элементарного
отрезка. Квадратурная формула называется
формулой средних прямоугольников
(рисунок 5.3) и имеет вид

.(5.35)

Рисунок
5.2. Интегрирование по формуле левых
прямоугольников

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

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

и

(рисунок 5.4). В этом случае площадь
-ой
элементарной трапеции с основанием

будет иметь вид

. (5.36)

Рисунок
5.3. Интегрирование по формуле средних
прямоугольников

Приближенное
значение интеграла в свою очередь будет
определяться как

. (5.37)

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

Рисунок
5.4. Интегрирование по формуле трапеций

В
качестве примера приведена программа
вычисления интеграла от функции

на отрезке [-5,5] по формуле правых
прямоугольников.

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

Пример 1. Длина карандаша измерена линейкой с миллиметровыми делениями. Измерение показало 17,9 см. Ошибка неизвестна, но она заведомо меньше чем 0,1 см. Поэтому можно принять 0,1 см за предельную погрешность. Относительная предельная погрешность равна . Округляя в сторону увеличения, находим 0,6%.

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

где есть предельная погрешность аргумента. Величина округляется в сторону увеличения (ввиду неточности самой формулы).

Пример 2. Сторона квадрата, найденная измерением, оказалась равной 46 м. Предельная погрешность

равна 0,1м. Найти предельную погрешность для площади квадрата.

Решение. Имеем: (х – сторона квадрата, у — площадь). Отсюда . В нашем примере Значит, Предельная абсолютная погрешность равна (округленно) Предельная относительная погрешность равна .

Предельную относительную погрешность можно найти также с помощью логарифмического дифференцирования (§ 244) по формуле

В частности, при тогда имеем:

т. е. предельная относительная погрешность степени хправна n-кратной предельной относительной погрешности аргумента.

Пример 3. При условиях примера 2 предельная относительная погрешность площади равна

Пример 4. Измерение ребра куба дало см. Предельная погрешность 0,05 см. Какова предельная относительная погрешность для объема куба?

Решение. Предельная относительная погрешность для равна ; для она равна .

Правило 1. Предельная относительная погрешность произведения двух или нескольких сомножителей равна сумме предельных относительных погрешностей сомножителей.

Лекция 1
Введение. Типы и источники погрешностей
при приближенных вычислениях
Целью курса «Численные методы» является рассмотрение вопросов, связанных с приближенными методами решения различных прикладных задач.
Численные методы разрабатывают и исследуют, как правило, высококвалифицированные специалисты по вычислительной математике. Что же касается прикладников, использующих численные методы как готовый инструмент в своей практической деятельности, то для них главной задачей является понимание основных идей методов, их особенностей и областей их применения.
Обычно решение прикладной задачи состоит из нескольких этапов.
На первом этапе формулируется содержательная постановка задачи, которую требуется решить. При этом решается вопрос, какие факторы следует учесть, а какими можно пренебречь.
На втором этапе содержательной модели ставится в соответствие математическая модель, т.е. математическое описание естественного процесса с помощью различного рода уравнений и неравенств.
Третий этап состоит в построении приближенного метода решения математической задачи, т.е. в выборе вычислительного алгоритма. Специалисту-прикладнику для решения задачи, как правило, необходимо выбрать из имеющегося арсенала методов тот, который наиболее пригоден в каждом конкретном случае.
На четвертом этапе осуществляется программирование вычислительного алгоритма для ЭВМ.
На пятом этапе готовятся исходные данные задачи и проводятся расчеты по составленной и отлаженной программе.
В курсе численного анализа мы в основном будем иметь дело с вопросами, касающимися третьего этапа, т.е. с численными методами решения математических задач.
При решении задачи на ЭВМ мы всегда получаем не точное, а приближенное решение. Мерой точности решения является допущенная при его получении погрешность. Можно выделить три основные причины возникновения погрешностей при численном решении задачи.
Прежде всего, исходные данные задачи (начальные и граничные условия, коэффициенты и правые части уравнений и т.п.) всегда задаются с некоторой погрешностью. Эту погрешность принято называть неустранимой.
Численный метод, используемый при решении математической задачи, также является источником погрешностей. Это связано, например, с заменой интегралов конечными суммами, усечением бесконечных рядов при вычислении значений функций и т.п. Как правило, погрешность метода регулируема, т.е. может быть уменьшена до некоторого разумного предела путем изменения параметров численного метода, например, уменьшения шага при численном интегрировании и т.д.
Наконец, в силу ограниченности разрядной сетки ЭВМ, т.е. числа разрядов, предназначенных для хранения в ЭВМ чисел, неизбежными являются погрешности округления. В процессе выполнения алгоритма погрешности округления накапливаются. Результирующая погрешность, вносимая в решение за счет погрешностей округления, называется вычислительной. Величина вычислительной погрешности определяется двумя факторами: точностью представления чисел в ЭВМ и чувствительностью алгоритма к погрешностям округления.
В связи со сказанным введем понятие устойчивости алгоритма. Алгоритм называется устойчивым, если в процессе его выполнения вычислительная погрешность возрастает незначительно, и неустойчивым – в противном случае.
В качестве примера неустойчивого алгоритма рассмотрим следующий метод вычисления определенного интеграла
.
Интегрируя по частям, находим:
,
.
Итак, для вычисления интеграла при всех натуральных получаем рекуррентную формулу
Пользуясь этой формулой, находим:
Очевидно, что значение интеграла не может быть отрицательным, поскольку подынтегральная функция неотрицательна на всем промежутке интегрирования . В чем же причина возникновения столь большой вычислительной погрешности? Погрешность при вычислении не превосходит 5 единиц последнего отброшенного разряда, т.е. . На каждом следующем шаге эта погрешность умножается на число, модуль которого больше . В итоге начальная погрешность, не превосходящая величины , умножается на , что и приводит к результату, не имеющему смысла. Таким образом, причиной возникновения столь большой вычислительной погрешности в данном примере является неустойчивость предложенного алгоритма к погрешностям округления.
Отметим, что в данном примере на каждом шаге алгоритма вычисления выполнялись абсолютно точно, и только один операнд содержал наследственную погрешность, которая появилась в результате неточного вычисления числа . В общем случае все операнды содержат погрешности, т.е. являются приближенными числами, и все операции выполняются приближенно. Поэтому естественным представляется вопрос о том, как погрешности отдельных операндов сказываются на точности результата той или иной операции. Чтобы ответить на этот вопрос, введем понятия абсолютной и относительной погрешности приближенного числа.
Пусть – точное, вообще говоря, неизвестное, значение некоторой величины, а – известное его приближение. Абсолютной погрешностью приближения называют такую величину , относительно которой известно, что
.
Таким образом, под абсолютной погрешностью приближенного значения величины понимается некоторая оценка модуля разности между приближенным и точным значениями этой величины.
Относительной погрешностью приближения называется такая величина , относительно которой известно, что
.
Очевидно, что в качестве можно взять, например,
.
Сформулируем основную задачу теории погрешностей. Пусть известны погрешности некоторой совокупности величин; требуется определить погрешность заданной функции этих величин. Конкретнее, пусть – дифференцируемая функция своих аргументов; – известные приближения точных значений аргументов соответственно, причем абсолютные погрешности этих приближений заданы:
.
Требуется оценить абсолютную погрешность значения как приближения к точному значению . Эта задача имеет следующее решение:
, (1)
где
,
а
,
т.е. – некоторое промежуточное значение частной производной функции по аргументу в точке, расположенной на отрезке, соединяющем точную точку и приближенную точку .
Так как точная точка и значение параметра нам неизвестны, то на практике пользуются приближенными оценками погрешности, следующими из формулы (1).
Приведем конкретные оценки погрешностей для простейших функций.
1. . В этом случае . Следовательно, для абсолютной погрешности суммы и разности имеет место оценка
.
Таким образом, абсолютная погрешность алгебраической суммы приближенных слагаемых не превосходит суммы абсолютных погрешностей слагаемых.
2. . Легко показать, что в этом случае из (1) и определения относительной погрешности следует оценка
.
Следовательно, относительная погрешность произведения и частного приближенно равна сумме относительных погрешностей аргументов.
3. . Для этой функции имеем оценку
,
т.е. абсолютная погрешность синуса от приближенного значения аргумента не превосходит абсолютной погрешности аргумента.
4. . В этом случае имеет место приближенная оценка
,
т.е. относительная погрешность квадратного корня из приближенного значения не превосходит относительной погрешности аргумента.
Лекции 2-4
Глава 1. Приближение функций
§ 1. Постановка задачи приближения функций
Пусть некоторая величина является функцией аргумента . На практике явная связь между и зачастую неизвестна, т.е. невозможно записать эту связь в виде некоторой формульной зависимости . В некоторых случаях даже при известной зависимости она настолько сложна, например, содержит трудно вычисляемые выражения, специальные функции, сложные интегралы и т.п., что ее использование в практических расчетах затруднительно.
Наиболее распространенным и практически важным случаем, когда зависимость между и неизвестна, является задание этой зависимости в виде таблицы , где – значения неизвестной функции в точках . Эти значения являются либо результатами предшествующих расчетов, либо экспериментальными данными. Пусть требуется найти значение величины в точках, отличных от точек , которые мы будем называть узлами. Этой цели и служит задача о приближении функций: неизвестную функцию требуется заменить некоторой функцией такой, что ее отклонение от функции является в некотором смысле наименьшим. Функция называется аппроксимирующей функцией.
Одним из основных способов приближения функций является так называемое интерполирование. Оно состоит в следующем: аппроксимирующая функция подбирается из условия совпадения ее значений в узлах с заданными значениями :
.
На практике наибольшее распространение получил случай интерполирования алгебраическими многочленами, т.е. случай, когда функция выбирается в виде
.
Условия интерполяции в этом случае принимают вид:
. (1.1)
Условия (1.1) представляют собой систему из линейных алгебраических уравнений относительно неизвестных . Мы знаем, что в случае, когда число уравнений больше числа неизвестных, система, как правило, несовместна (она является переопределенной). В случае, когда число уравнений меньше числа неизвестных, система является неопределенной. Обычно такие системы имеют бесчисленное множество решений. Поэтому естественно рассмотреть случай, когда число уравнений в системе (1.1) совпадает с числом неизвестных, т.е. случай, когда , или . В этом случае аппроксимирующая функция является многочленом степени , т. е.
а система (1) принимает вид
. (1.2)
Для выяснения вопроса о разрешимости системы (1.2) выпишем ее основной определитель:
.
Это известный из курса алгебры определитель Вандермонда, он равен
.
Ясно, что тогда и только тогда, когда все узлы различны.
Таким образом, если все узлы интерполяции попарно различны, то система (1.2) имеет и притом единственное решение. Тем самым доказаны существование и единственность интерполяционного многочлена степени , совпадающего в заданных точках с заданными значениями .
Поскольку многочлен используется для приближения функции на всем рассматриваемом интервале изменения переменной , то в этом случае говорят о глобальной интерполяции.
Если количество узлов интерполяции велико, то глобальная интерполяция приводит к многочленам высокой степени, что не совсем удобно с вычислительной точки зрения, так как вычисление значений многочленов высокой степени сопровождается большой вычислительной погрешностью. Поэтому часто интерполяционные многочлены строят отдельно для различных небольших частей рассматриваемого интервала изменения аргумента . В этом случае говорят о кусочной, или локальной, интерполяции.
В любом случае при интерполировании основным условием является прохождение графика аппроксимирующей функции через данные табличные точки. Однако в ряде случаев выполнение этого условия затруднительно или даже нецелесообразно. Допустим, например, что табличные значения получены путем измерений и, следовательно, содержат погрешности. Построение аппроксимирующей функции с условием обязательного прохождения ее графика через экспериментальные точки означало бы тщательное повторение допущенных при измерениях погрешностей.
В подобных случаях часто используется другой подход к построению аппроксимирующей функции, называемый методом наименьших квадратов. Мы рассмотрим идею этого метода для случая, когда в качестве аппроксимирующей функции выбирается многочлен
степени . Мерой отклонения многочлена от функции на множестве точек назовем величину
,
равную сумме квадратов разностей между значениями аппроксимирующего многочлена и приближаемой функции в узлах. Согласно методу наименьших квадратов коэффициенты многочлена подбираются так, чтобы функция
принимала наименьшее возможное значение. На практике обычно ограничиваются многочленами небольшой степени, чаще всего 1, 2 или 3.
Упражнение. Методом наименьших квадратов построить аппроксимирующие многочлены первой и второй степени для функции , заданной таблицей , где – значения функции в узлах .
§ 2. Интерполяционный многочлен Лагранжа
Мы доказали существование и единственность интерполяционного многочлена степени
, (1.3)
совпадающего с заданными значениями в узлах интерполяции . Его коэффициенты являются решением системы
.
На практике нет необходимости решать эту систему, тем более что ее решение является очень чувствительным к погрешностям округления. Как мы сейчас увидим, для интерполяционного многочлена (1.3) можно выписать ряд явных представлений.
Рассмотрим многочленов , каждый из которых имеет степень и удовлетворяет следующим условиям:
(1.4)
Условия (1.4) однозначно определяют многочлены . Действительно, так как при , то точки являются корнями многочлена . По теореме Безу это означает, что многочлен делится на двучлены , т.е. его можно представить в виде
.
Заметим, что степень многочлена по определению должна быть равна . Так как степень произведения , очевидно, уже равна , то . Значение этой постоянной определим из условия . Имеем:
.
Тогда
.
Покажем теперь, что многочлен
(1.5)
является искомым интерполяционным многочленом. Действительно, он имеет степень не выше чем , так как является линейной комбинацией многочленов степени . Далее для всех имеем
,
т.е. многочлен совпадает с заданными значениями функции в узлах интерполяции. Полученное выражение (1.5) интерполяционного многочлена называется интерполяционным многочленом Лагранжа.
В дальнейшем нам потребуется оценка остаточного члена
интерполяционного многочлена. Пусть функция раз непрерывно дифференцируема. Тогда имеет место следующее представление:
.
Здесь
,
а – некоторая промежуточная точка, расположенная между точками и .
§ 3. Разделенные разности. Интерполяционный многочлен Ньютона
с разделенными разностями
Пусть дана совокупность узлов . Значения функции в этих узлах будем называть разделенными разностями нулевого порядка.
Разделенными разностями первого порядка назовем величины
.
Разделенными разностями второго порядка называются величины
.
Разделенными разностями -го порядка назовем величины
.
Таким образом, разделенные разности более высокого порядка определяются рекуррентно через разделенные разности порядка, меньшего на . Вместе с тем можно выписать выражение разделенной разности любого порядка непосредственно через значения функции в узлах, т.е. через разности нулевого порядка. А именно, имеет место формула
. (1.6)
Поскольку значение суммы не зависит от порядка слагаемых, то из формулы (1) следует, что разделенная разность не изменяется при перестановке ее аргументов, например,
.
Все разделенные разности, построенные по совокупности узлов , принято располагать в виде следующей треугольной таблицы:
Перейдем теперь к другой форме записи интерполяционного многочлена, использующей разделенные разности. Преобразуем остаточный член многочлена Лагранжа:
=
.
Согласно формуле (1.6) выражение в скобках есть разделенная разность , построенная по узлам . Таким образом,
. (1.7)
Обозначим через интерполяционный многочлен степени , построенный по узлам . Представим многочлен в виде
(1.8)
и рассмотрим разность . Она представляет собой многочлен степени , обращающийся в в точках . Действительно, по определению интерполяционного многочлена
,
.
Следовательно, при имеем:
.
Отсюда по теореме Безу следует, что многочлен делится на двучлены . Значит, многочлен можно записать в виде
.
Заметим, что степень многочлена равна , т.е. совпадает со степенью многочлена . Поэтому . Для определения этой постоянной вычислим значение многочлена в узле :
.
Но , поэтому полученное равенство можно переписать в виде
. (1.9)
С другой стороны, воспользуемся равенством (1.7) при и :
,
или
. (1.10)
Сравнивая равенства (1.9) и (1.10), получаем:
.
Следовательно,
.
Замечая, что , и подставляя полученные выражения в (1.8), приходим к интерполяционному многочлену Ньютона с разделенными разностями:
.
В заключение этого параграфа сравним формулу (1.7) и представление остаточного члена многочлена Лагранжа, приведенное в конце предыдущего параграфа. В результате приходим к формуле
,
где . Таким образом, оказывается, что разделенная разность порядка лишь числовым множителем отличается от значения производной порядка в некоторой промежуточной точке. Отсюда, в частности, следует, что если функция является многочленом степени , то разделенные разности порядка постоянны, а все разделенные разности более высоких порядков равны нулю.
§ 4. Интерполирование с кратными узлами
Иногда приходится решать более общую задачу интерполирования, когда известны не только значения некоторой функции в узлах , но и значения ее производных до некоторого порядка. Строго эта задача формулируется следующим образом. Пусть требуется построить многочлен степени , удовлетворяющий условиям
(1.11)
где все узлы попарно различны и
. (1.12)
Числа называются кратностями узлов соответственно.
Можно показать, что поставленная задача имеет и притом единственное решение. Это решение можно записать в следующем виде:
.
Выписанное выражение называется интерполяционным многочленом Эрмита.
При его записи мы использовали так называемые разделенные разности с кратными узлами, например, . Ясно, что обычное определение разделенных разностей для вычисления не подходит. Разделенные разности с кратными узлами определяются следующим образом. Если все аргументы разделенной разности одинаковы, то следует воспользоваться формулой
.
Если же среди аргументов есть различные, то разделенную разность можно свести к разделенным разностям более низкого порядка с помощью таких же формул, что и в случае, когда все узлы различны:
=
.
В частности,
,
,
.
Можно показать, что для остаточного члена многочлена Эрмита имеет место представление, аналогичное представлению остаточного члена многочлена Ньютона:
,
где
.
Имеет место и представление, аналогичное представлению остаточного члена многочлена Лагранжа:
,
где .
Пример. Построить многочлен Эрмита , удовлетворяющий условиям
Решение. Здесь . Имеем:
.
Так как
,
то окончательно получаем:
.
§ 5. Конечные разности. Интерполяционные формулы для случая
равноотстоящих узлов
Пусть узлы расположены друг от друга на равном расстоянии :
.
Число называется шагом таблицы узлов.
Соответствующие этим узлам значения функции будем обозначать через и называть конечными разностями нулевого порядка.
Конечными разностями первого порядка называются величины
,
которые в зависимости от обстановки обозначаются и именуются по-разному:
– центральная разность;
– разность вперед;
– разность назад.
В общем случае конечные разности -го порядка определяются через разности порядка следующим образом:
;
;
.
Конечные разности принято располагать в виде треугольной таблицы:
Имеют место следующие свойства конечных разностей.
1. Конечные разности -го порядка выражаются через значения функции в узлах по формуле
,
где при четном число – целое, а при нечетном число – полуцелое.
В частности,
,
,
.
2. Если , то имеет место следующая связь между разделенными и конечными разностями:
. (1.13)
Воспользуемся указанной в конце § 3 связью между разделенной разностью и производной
, (1.14)
где . Сравнивая (1) и (2), получаем:
.
Следовательно, конечная разность порядка лишь числовым множителем отличается от значения производной порядка в некоторой промежуточной точке.
Отсюда, в свою очередь, следует, что если – многочлен степени , то конечные разности порядка и выше будут равны нулю.
Перейдем к построению интерполяционных формул, использующих конечные разности. Если данная таблица состоит из большого количества узлов, нам нет необходимости использовать их все. Это привело бы к увеличению степени интерполяционного многочлена, а значит, и к увеличению объема вычислений. На практике чаще всего достаточно использовать интерполяционный многочлен невысокой степени. Однако с целью достижения приемлемой точности интерполирования желательно выбирать те узлы таблицы, которые ближе всего расположены к точке, в которой выполняется интерполирование. Действительно, в выражении остаточного члена
фигурирует величина , равная произведению расстояний от точки , в которой выполняется интерполирование, до узлов таблицы. Очевидно, что величина будет минимальной, если изо всех имеющихся узлов использовать ближайших к точке . Если точка находится внутри таблицы узлов, то выгодно взять половину узлов слева от и половину – справа от . Если же точка находится в начале или в конце таблицы, то такой совокупности узлов выбрать нельзя. Сейчас мы рассмотрим как раз такие случаи.
Пусть точка расположена в начале таблицы узлов . Построим интерполяционный многочлен Ньютона с разделенными разностями по этой совокупности узлов:
.
Сделаем замену переменной
и учтем, что
.
Тогда
.
Переходя в выражении для от разделенных разностей к конечным, получаем:
.
Полученная формула называется первой интерполяционной формулой Ньютона, или формулой Ньютона для интерполирования вперед.
Рассмотрим теперь случай, когда точка находится в конце таблицы узлов . В этом случае удобнее привлекать узлы в следующей последовательности: . Построим по этой совокупности интерполяционный многочлен Ньютона с разделенными разностями:
.
Сделаем замену переменной
и учтем, что
.
Тогда
.
Перейдем от разделенных разностей к конечным:
,
,
,
.
Подставляя эти формулы в выражение для , получаем вторую интерполяционную формулу Ньютона, или формулу Ньютона для интерполирования назад:
.
Пусть теперь точка расположена внутри таблицы узлов. Обозначим через узел, ближайший к точке . Если , то оптимальным набором узлов будет совокупность
.
Построим по этой совокупности интерполяционный многочлен Ньютона с разделенными разностями:
.
Сделаем замену переменной
и учтем, что
.
Тогда
.
Воспользуемся связью между разделенными и конечными разностями:
,
,
,
,
,
.
Следовательно,
.
Полученная формула называется первой интерполяционной формулой Гаусса, или формулой Гаусса для интерполирования вперед.
Если же , то оптимальной будет совокупность узлов
.
В этом случае получаем формулу
,
которая называется второй интерполяционной формулой Гаусса, или формулой Гаусса для интерполирования назад.
§ 6. Интерполирование с помощью кубических сплайнов
Пусть функция определена на отрезке . На этом отрезке рассмотрим сетку , т.е. множество узлов
.
Функция называется интерполяционным кубическим сплайном для функции относительно сетки , если она удовлетворяет трем условиям:
1) , т.е. непрерывна на отрезке вместе со своими производными до второго порядка включительно;
2) на каждом частичном отрезке она является многочленом третьей степени, который мы обозначим через ;
3) .
Замечание. Условия 1)-2) определяют кубический сплайн; условие 3) есть условие интерполирования сплайном функции в узлах сетки .
Условие 2) означает, что при
.
Многочлен полностью определяется четырьмя коэффициентами . Поскольку общее число промежутков равно , то для полного определения сплайна нужно найти значения неизвестных .
Пользуясь определением сплайна, выпишем уравнения, из которых определяются эти коэффициенты.
Многочлен является функцией, бесконечное число раз непрерывно дифференцируемой, поэтому условие 1) сводится к требованию непрерывности самого сплайна и его первой и второй производных и во внутренних узлах , т.е. в точках, являющихся общими для соседних промежутков и :
, (1.15)
, (1.16)
. (1.17)
В силу условия 3) определения сплайна значения и не только совпадают между собой, но и равны значению . Поэтому, заменив соотношения (1.15) равенствами
, (1.18)
, (1.19)
мы учтем как условие непрерывности сплайна во внутренних узлах , так и условие интерполяционности сплайна в этих узлах.
К соотношениям (1.18) добавим условие интерполяционности сплайна в узле
,
а к соотношениям (1.19) – условие интерполяционности в узле
.
Соберем вместе все полученные уравнения:
, (1.20)
, (1.21)
, (1.22)
. (1.23)
Соотношения (1.20)-(1.23) представляют собой систему из линейных уравнений относительно неизвестных. Для замыкания этой системы требуется еще два уравнения. Для их получения привлекают некоторые дополнительные условия, обычно индуцируемые физической постановкой задачи и называемые краевыми условиями.
Рассмотрим наиболее употребительные краевые условия.
а) , где и – известные числа;
б) , где и – известные числа;
в) ,
где и – кубические многочлены, графики которых проходят соответственно через первые четыре и последние четыре узла сетки;
г) если известно, что приближаемая функция является периодической с периодом , то и сплайн также должен быть периодическим; в этом случае краевые условия задают в виде
.
В любом из перечисленных случаев мы получаем два дополнительных уравнения, необходимых для замыкания системы (1.20)-(1.23).
На практике нет необходимости в явном виде составлять и затем решать полученную систему. Сейчас мы покажем, что сплайн полностью определяется значениями своей первой или второй производной в узлах сетки, что позволяет вместо системы из уравнений относительно неизвестных решать систему только из уравнений с таким же числом неизвестных.
Введем обозначения:
,
.
Числа называются наклонами сплайна , а числа – его моментами.
Покажем, что сплайн целиком определяется, например, своими моментами . Для этого достаточно показать, что все коэффициенты можно выразить через моменты , а также через данные значения и .
На отрезке сплайн имеет вид
.
Коэффициенты найдем из условий:
, (1.24)
. (1.25)
Имеем:
.
Введем обозначение
.
Тогда из (1.24) находим
,
откуда
.
Из (1.25) находим
,
откуда
.
В результате получаем следующее выражение для через моменты:
.
Обычно выражение для преобразуют и используют в следующем виде:
.
В таком виде выражение для многочлена легче запоминается и может быть эффективно использовано для целей интерполирования функции , численного дифференцирования и интегрирования, если, конечно, предварительно определены моменты .
Для определения моментов воспользуемся пока еще не использованным условием непрерывности первой производной сплайна во внутренних узлах и краевыми условиями. Имеем:
.
Аналогично
.
Запишем теперь условия (1.20):
,
или
.
После умножения полученного равенства на оно принимает вид
.
Так как условия (1.20) должны быть выполнены для всех , то мы получили систему из уравнений для определения неизвестных . Недостающие два уравнения дает использование краевых условий.
Наиболее простые уравнения получаются при использовании краевых условий типа б). В этом случае они имеют вид
,
где числа и заданы.
В случае использования краевых условий типа а) дополнительные уравнения получаются немного сложнее. Имеем:
,
.
В результате элементарных преобразований получаем:
,
.
После умножения первого из этих равенств на , а второго – на они принимают вид
,
.
Аналогичным образом получают дополнительные уравнения и при использовании других краевых условий.
Итак, в случае использования краевых условий типа б) для определения моментов приходим к системе
(1.26)
а в случае краевых условий типа а) получаем систему
(1.27)
И система (1.26) и система (1.27) являются системами линейных алгебраических уравнений с трехдиагональной матрицей вида
,
причем ее элементы удовлетворяют условиям:
,
.
Последние условия означают, что матрица обладает свойством диагонального преобладания. Известно, что такая матрица невырождена. Следовательно, решения обеих систем (1.26) и (1.27) существуют и определяются единственным образом.
§ 7. Метод прогонки решения системы с трехдиагональной матрицей
Опишем эффективный метод решения системы с трехдиагональной матрицей, имеющей диагональное преобладание, – метод прогонки. Рассмотрим систему
(1.28)
Из первого уравнения выразим :
.
Введем обозначения:
.
Тогда запишется в виде
.
Предположим, что
,
и подставим это выражение в -ое уравнение системы:
.
Отсюда выразим :
.
Обозначим:
.
Тогда . Тем самым мы доказали, что для всех справедливо представление
,
и получили рекуррентные формулы для так называемых прогоночных коэффициентов , :
(1.29)
При имеем:
.
Подставим это выражение в последнее уравнение системы:
.
Отсюда найдем :
. (1.30)
Затем для находим по формулам
. (1.31)
Вычисление прогоночных коэффициентов , по формулам (1.29) называется прямым ходом метода прогонки; вычисление неизвестных системы (1.28) по формулам (1.30)-(1.31) называется обратным ходом метода прогонки.
Известно, что выполнение условий диагонального преобладания является не только достаточным для осуществления метода прогонки, но и обеспечивает вычислительную устойчивость этого метода, заключающуюся в том, что погрешности, допускаемые на отдельных шагах метода, не увеличиваются на следующих шагах.
Лекция 5
Глава 2. Численное дифференцирование
Простейшие формулы численного дифференцирования, т.е. формулы для приближенного вычисления значений производных таблично заданных функций, получаются дифференцированием интерполяционных формул.
Пусть функция задана таблично, т.е. известны ее значения в узлах , и – ее интерполяционный многочлен. Положим
. (2.1)
Получим представление остаточного члена построенной формулы (2.1). Имеем:
,
где – остаточный член интерполяционной формулы . По формуле Ньютона-Лейбница получаем:
.
Обозначим: . По определению разделенной разности с кратными узлами имеем:
.
Используя формулу связи между разделенной разностью и производной, находим:
,
где – некоторая точка на промежутке , , . Окончательно получаем следующее представление остаточного члена формулы численного дифференцирования (2.1):
. (2.2)
Отсюда легко следует оценка остаточного члена формулы численного дифференцирования (2.1):
. (2.3)
На практике особый интерес представляют формулы численного дифференцирования для функций, заданных таблично на равномерной сетке, т.е. на сетке, узлы которой являются равноотстоящими друг от друга.
В качестве примера получим приближенную формулу для вычисления значения первой производной с погрешностью на основе первой интерполяционной формулы Ньютона. Здесь – шаг сетки. Выпишем первую интерполяционную формулу Ньютона для узлов , где , :
.
Соответствующая ей формула численного дифференцирования есть
.
Остаточный член этой формулы согласно (2.2) имеет вид
.
Имеем:
,
.
Следовательно,
.
Поскольку мы поставили задачу построения формулы численного дифференцирования, имеющей погрешность , то следует принять . В таком случае мы должны дифференцировать многочлен
.
Имеем:
.
Так как , то , и, следовательно,
.
Это и есть искомая формула для приближенного вычисления первой производной с погрешностью .
Обычно нас интересуют приближенные значения производной таблично заданной функции в узлах таблицы. Например, пусть . Тогда , и мы получаем:
.
Итак, с погрешностью имеем формулу численного дифференцирования
.
Совершенно аналогичным образом на основе различных интерполяционных формул с конечными разностями можно получить другие формулы численного дифференцирования с различной величиной погрешности. Приведем без вывода несколько таких формул:
с погрешностью ;
с погрешностью ;
с погрешностью ;
с погрешностью .
Лекции 6-7
Глава 3. Численное интегрирование
§ 1. Квадратуры Ньютона-Котеса
Пусть требуется вычислить определенный интеграл
, (3.1)
где – некоторая интегрируемая функция.
Мы рассмотрим ряд методов приближенного вычисления интеграла , в которых интеграл заменяется линейной комбинацией значений функции в конечном числе точек отрезка интегрирования . Другими словами, мы будем рассматривать методы, в которых
. (3.2)
Формулы типа (3.2) называются квадратурными, числа – коэффициентами квадратуры (3.2), а точки – узлами квадратуры.
Простейшие квадратурные формулы получаются при помощи аппарата интерполирования.
На отрезке выберем узлов . С помощью замены переменной
(3.3)
отрезок переходит в отрезок , а точки – в некоторые точки . Таким образом, выбор узлов на отрезке индуцирует некоторый выбор узлов на отрезке . По узлам построим интерполяционный многочлен для функции . В общем случае некоторые из узлов будем считать кратными. Напомним, что узел называется кратным, если в нем требуется совпадение значений не только самой функции, но и ее производных до некоторого порядка . При этом число называется кратностью узла . Как и ранее, через обозначим суммарную кратность всех узлов, т.е.
.
Итак, в общем случае будем строить интерполяционный многочлен Эрмита степени . Положим теперь
.
Получаемые при этом формулы называются квадратурными формулами Ньютона-Котеса.
Рассмотрим остаточный член квадратуры Ньютона-Котеса:
.
Здесь – остаточный член интерполяционного многочлена Эрмита. Воспользуемся известным из теории интерполирования представлением
,
где
.
Так как все узлы принадлежат отрезку , то для всех справедлива оценка
,
где
.
Следовательно,
.
Воспользуемся заменой (3.3). Имеем:
,
.
Здесь
.
Тогда
.
Обозначим:
.
Тогда окончательно получаем следующую оценку:
. (3.4)
Рассмотрим теперь частный случай, когда все узлы являются простыми. В этом случае
;
.
Снова воспользуемся заменой переменной (3.3). Имеем:
,
.
Следовательно,
.
В результате получаем следующую квадратуру
, (3.5)
где
. (3.6)
Заметим, что поскольку при построении квадратуры (3.5) все узлы являлись простыми, то для нее оценка (3.4) остаточного члена принимает вид
.
Из этой оценки следует, что квадратура (3.5) точна в случае, когда функция является произвольным многочленом степени , ибо для такой функции производная порядка тождественно равна нулю.
§ 2. Симметричные квадратуры Ньютона-Котеса
На практике узлы квадратурной формулы обычно выбирают так, чтобы они располагались симметрично относительно середины отрезка интегрирования :
.
Для соответствующих им узлов на отрезке это означает, что
.
В дальнейшем будем рассматривать квадратуры именно с таким расположением узлов. Такие квадратуры называются симметричными.
По сравнению с общим случаем симметричные квадратуры Ньютона-Котеса (3.5), в которых коэффициенты определяются по формулам (3.6), обладают рядом дополнительных свойств. Рассмотрим эти свойства.
1. Коэффициенты квадратуры, соответствующие симметрично расположенным узлам, равны:
.
Доказательство. В интеграле
,
определяющем коэффициент , сделаем замену переменной . Тогда
,
,
.
2. Симметричная квадратура точна для любой нечетной относительно середины отрезка интегрирования функции.
Доказательство. Пусть функция нечетна относительно середины отрезка :
.
В этом случае
.
С другой стороны,
.
Если в этой сумме число членов четно, то все они взаимно уничтожаются. Если же нечетно, то центральное слагаемое остается, но оно равно нулю в силу нечетности функции относительно точки . Таким образом, в любом случае имеем:
.
Следовательно,
,
а это и означает, что квадратура точна.
3. Симметричная квадратура с нечетным числом узлов точна для всех многочленов степени .
Доказательство. Рассмотрим произвольный многочлен степени
.
Очевидно, что его можно представить в виде
,
где – некоторый многочлен степени . Покажем, что первое слагаемое – функция, нечетная относительно точки :
.
В силу свойства 2 симметричная квадратура точна для такой функции, т.е.
.
Так как – многочлен степени , то для него будет точна любая квадратура Ньютона-Котеса, построенная по узлам. Следовательно, и
.
Остается показать, что остаточный член квадратуры Ньютона-Котеса является аддитивной функцией, т.е. для любых функций и имеет место равенство
.
Имеем:
.
Отсюда следует, что
.
Тем самым свойство 3 доказано.
Естественно ожидать, что оценка остаточного члена симметричной квадратуры при нечетном должна быть более тонкой, чем в общем случае. Именно с целью получения более тонких оценок остаточного члена в данном и ему подобных случаях мы и допустили с самого начала, что некоторые узлы могут быть кратными. Конкретнее, в случае, когда квадратура симметрична, а число узлов нечетно, мы будем использовать центральный узел как двукратный, а все остальные узлы – как простые. При этом оказывается, что квадратура, получающаяся при замене функции многочленом Эрмита , совпадает с квадратурой, построенной только по простым узлам, т.е. получающейся при замене функции многочленом Лагранжа . Действительно, пусть узел – двукратный, а остальные узлы – простые. Построим многочлен Эрмита по совокупности узлов
:
.
Покажем, что функция в силу нечетности числа и симметричного расположения узлов нечетна относительно средней точки :
.
Отсюда следует, что
.
Поэтому
.
Таким образом, мы показали, что если число узлов нечетно и квадратура симметрична, то она может быть получена двумя способами: 1) заменой функции многочленом Лагранжа и 2) заменой функции многочленом Эрмита . Соответствующие этим способам получения квадратуры оценки остаточного члена имеют вид:
, (3.7)
. (3.8)
Естественно из двух оценок, справедливых для одной и той же формулы, пользоваться более тонкой оценкой (3.8).
§ 3. Употребительные квадратуры Ньютона-Котеса
1. Формула левосторонних прямоугольников.
Положим , , . Следовательно, . Имеем:
.
Поскольку
,
то формула левосторонних прямоугольников имеет вид
.
Согласно (3.7) оценка остаточного члена этой формулы имеет вид
.
2. Формула правосторонних прямоугольников.
В этой формуле , , , . Имеем:
.
Поскольку
,
то формула правосторонних прямоугольников имеет вид
.
Согласно (3.7) оценка остаточного члена этой формулы имеет вид
.
3. Формула центральных прямоугольников.
В этой формуле , , , . Имеем:
.
Так как
,
то формула принимает вид
.
В силу (3.8) оценка остаточного члена этой формулы имеет вид
.
4. Формула трапеций.
В этой формуле , , , , , . Имеем:
,
.
Так как
,
мы получаем формулу
,
оценка остаточного члена которой согласно (3.7) имеет вид:
.
5. Формула Симпсона.
В формуле Симпсона , , , , , , , . Тогда
,
,
.
Мы получаем формулу
с оценкой остаточного члена
.
Заметим, что при малой длине отрезка интегрирования эта оценка на два порядка лучше оценок формул центральных прямоугольников и трапеций.
§ 4. Составные квадратуры Ньютона-Котеса
До сих пор, сравнивая квадратурные формулы по порядку погрешности, мы предполагали, что длина отрезка интегрирования мала (меньше 1). Если же это не так, то рассмотренные формулы применяют к отдельным малым по длине частям этого отрезка. Получаемые при этом формулы называются составными квадратурными формулами.
Разобьем отрезок на частей точками
и для вычисления интеграла по каждой из частей применим одну из квадратурных формул Ньютона-Котеса:
.
Здесь
.
Пусть – кратность узла . Тогда оценка остаточного члена полученной составной формулы имеет вид:
.
Заметим, что, как и в простых формулах, здесь при четном все , а при нечетном и симметричном расположении узлов равны 1 все , кроме одного: .
В простейшем и наиболее употребительном случае точки выбираются равноотстоящими:
.
Тогда полученная формула принимает вид
,
а оценка остаточного члена – вид:
.
Приведем конкретные примеры составных квадратурных формул.
1. Составная формула левосторонних прямоугольников:
.
Оценка остаточного члена составной формулы левосторонних прямоугольников имеет вид
.
2. Составная формула правосторонних прямоугольников:
.
Оценка остаточного члена составной формулы правосторонних прямоугольников совпадает с оценкой погрешности составной формулы левосторонних прямоугольников:
.
3. Составная формула центральных прямоугольников:
.
Здесь
.
Оценка остаточного члена составной формулы центральных прямоугольников имеет вид
.
4. Составная формула трапеций:
.
Оценка ее остаточного члена:
.
5. Составная формула Симпсона:
.
Оценка ее остаточного члена:
.
Замечание. Обычно формулу Симпсона записывают в терминах расстояний между соседними узлами квадратуры
:
.
Оценка остаточного члена принимает в этом случае вид:
.
§ 5. Правило Рунге
Обозначим:
.
Зададим некоторое целое и будем последовательно вычислять при . После каждого вычисления найдем величину
,
где – суммарная кратность узлов используемой простой формулы, и проверим, выполнено ли неравенство
,
где – требуемая точность вычисления интеграла .
Если при некотором это неравенство выполняется, то в качестве приближенного значения с требуемой степенью точности принимаем или, что, как правило, лучше, .
Описанное правило достижения требуемой точности вычисления определенного интеграла называется правилом Рунге.
§ 5. Численное интегрирование с помощью кубических сплайнов
Рассмотрим другой способ приближенного вычисления определенного интеграла
.
На отрезке рассмотрим сетку
.
Пусть сетка является равномерной, т.е. ее узлы находятся на равном расстоянии друг от друга:
.
Для функции построим интерполяционный кубический сплайн относительно сетки . Заменяя подынтегральную функцию ее сплайном , получим приближенную формулу
. (3.9)
Имеем:
.
Здесь – кубический многочлен, представляющий сплайн на частичном отрезке .
Напомним, что интерполяционный кубический сплайн однозначно определяется своими наклонами или своими моментами . Подробно рассмотрим случай, когда сплайн задается моментами. В этом случае кубический многочлен , представляющий сплайн на частичном отрезке , имеет вид
.
Подставим это выражение в интеграл :
.
Тогда
.
Окончательно формула (3.9) принимает вид
.
Лекции 6-7
Глава 5. Прямые методы решения систем
линейных алгебраических уравнений
§ 1. Прямые методы решения СЛАУ. Метод Гаусса
Мы рассмотрим методы решения СЛАУ двух типов – прямые и итерационные. В прямых методах точное решение системы при отсутствии погрешностей округления находится за конечное число арифметических действий. В итерационных методах определяется некоторая последовательность приближений, которая в пределе стремится к точному решению системы.
Одним из основных прямых методов решения СЛАУ является метод Гаусса. Под методом Гаусса понимается целая совокупность методов, в основе которых лежит использование элементарных преобразований матрицы системы с целью приведения ее к треугольному виду. Мы рассмотрим одну из возможных схем этого метода, в результате применения которой матрица системы принимает верхнюю треугольную форму.
Рассмотрим систему линейных уравнений относительно неизвестных . Предположим, что все угловые миноры матрицы отличны от нуля:
.
На первом шаге матрица преобразуется к матрице , все элементы первого столбца которой, начиная со второго, равны нулю. Согласно предположению угловой минор первого отличен от нуля. Поэтому мы можем определить множители первого шага
.
Для всех вычтем из -ой строки матрицы первую, умноженную на . В результате получим матрицу требуемого вида:
.
При этом значения элементов находятся по формулам
.
Элемент , на который приходится делить при определении множителей первого шага, называется ведущим элементом первого шага.
Отметим для дальнейшего, что описанные преобразования матрицы не изменяют значений ее угловых миноров. Поэтому все угловые миноры полученной матрицы отличны от нуля. В частности, отличен от нуля угловой минор второго порядка этой матрицы. Но он равен произведению . Следовательно, отличен от нуля элемент , который мы назовем ведущим элементом второго шага.
Целью второго шага является получение нулей во втором столбце матрицы , начиная с третьего элемента этого столбца. Определим множители второго шага
,
и для всех из -ой строки матрицы вычтем вторую, умноженную на . Ясно, что описанные преобразования не изменяют первых двух строк и первого столбца матрицы . В результате приходим к матрице
,
где
Пусть в результате выполнения шагов мы пришли к матрице
.
Так как преобразования, выполняемые на каждом из шагов, не меняют значений угловых миноров, то все угловые миноры матрицы совпадают с соответствующими угловыми минорами исходной матрицы . В частности, отличен от нуля угловой минор -ого порядка матрицы . Но он равен произведению . Следовательно, . Назовем элемент ведущим элементом -ого шага и перейдем к описанию этого шага. Его цель – получение нулей в -ом столбце матрицы , начиная с -ого элемента. Определим множители -ого шага
,
и для всех из -ой строки матрицы вычтем ее -ую строку, умноженную на . Так как при этом первые строк и первые столбцов матрицы остаются неизменными, то мы приходим к матрице вида
,
где .
Ясно, что при мы придем к матрице , имеющей верхнюю треугольную форму. В силу невырожденности исходной матрицы матрица также невырождена.
Описанный процесс приведения матрицы системы к верхней треугольной форме называется прямым ходом метода Гаусса.
Отметим, что мы описали процесс преобразования только матрицы системы . Для того чтобы получающаяся в результате система оставалась эквивалентной исходной системе, следует, очевидно, с компонентами правой части выполнить те же самые преобразования, что и со строками матрицы . Так, на первом шаге для всех из следует вычесть первую компоненту , умноженную на соответствующий множитель . В результате мы придем к вектору
,
где .
Точно таким же образом поступаем на всех следующих шагах. Именно, пусть в результате выполнения шагов мы пришли к вектору
.
На -ом шаге для всех из компоненты вычтем компоненту , умноженную на . В результате получим вектор
,
где .
При придем к окончательно преобразованной правой части . При этом система преобразуется к эквивалентной ей системе .
В дальнейшем матрицу будем обозначать через , а вектор – через . Так как матрица системы является верхней треугольной, то решение полученной системы находится очень просто. Последнее уравнение системы имеет вид , причем . Отсюда находим
.
Пусть значения неизвестных уже определены. Подставляя их в -ое уравнение системы
и учитывая, что , находим
.
Описанный процесс решения системы с верхней треугольной матрицей называется обратным ходом метода Гаусса.
§ 2. -разложение матрицы
При выполнении прямого хода метода Гаусса мы выполняли элементарные преобразования со строками матрицы . Выполнение каждой такой операции равносильно умножению матрицы слева на некоторую невырожденную матрицу. Действительно, выполнение первого шага, т.е. вычитание из -ой строки матрицы ее первой строки, умноженной на , проделанное для всех , равносильно умножению матрицы слева на матрицу
.
Таким образом, .
Преобразования второго шага равносильны умножению матрицы слева на матрицу
.
Поэтому .
По аналогии легко выписать матрицу , осуществляющую преобразования, выполняемые на -ом шаге:
.
При этом .
В результате выполнения шагов мы приходим к невырожденной верхней треугольной матрице . Легко выписать выражение матрицы через исходную матрицу :
.
Обозначив произведение через , получим
.
В силу невырожденности множителей матрица невырождена. Обозначим через матрицу, обратную к . Оказывается, что матрица имеет вид
.
Доказательство этого утверждения состоит в непосредственной проверке равенства , где – единичная матрица.
Представление матрицы в виде произведения нижней треугольной матрицы с единицами на главной диагонали и невырожденной верхней треугольной матрицы называется ее -разложением.
-разложение играет важную роль в численных методах линейной алгебры. Так, если для матрицы системы получено -разложение , то решение этой системы сводится к последовательному решению двух систем с треугольными матрицами
.
При наличии -разложения матрицы особенно выгодным является решение большого количества систем с одной и той же матрицей и различными правыми частями. Дело в том, что по количеству необходимых арифметических действий вычисление -разложения является более трудоемкой операцией, чем решение систем с треугольными матрицами. Конкретнее, вычисление -разложения требует арифметических действий, а решение системы с треугольной матрицей – .
К решению большого числа систем с одной и той же матрицей сводится, например, задача обращения матрицы. Пусть требуется найти матрицу , обратную к , т.е. решить матричное уравнение . Обозначив через -ый столбец матрицы , матричное уравнение запишем в виде совокупности систем уравнений
.
Таким образом, для вычисления обратной матрицы требуется решить систем уравнений с одной и той же матрицей .
В заключение параграфа упомянем одно достоинство -разложения, связанное с возможностью экономичного использования машинной памяти в процессе его вычисления на ЭВМ. Вспомним, что на первом шаге мы получаем нули в первом столбце, начиная со второго элемента этого столбца. Этих нулей ровно столько, сколько имеется множителей первого шага, определяющих первый столбец матрицы . Поскольку хранить нули в памяти нет необходимости, то на их местах можно как раз разместить множители первого шага. Заметим также, что первая строка матрицы , которая вообще не изменяется, является в то же время и первой строкой результирующей матрицы . Таким образом, в результате выполнения первого шага первая строка и первый столбец двумерного массива , в котором первоначально размещалась матрица системы, содержат полную информацию о первой строке матрицы и первом столбце матрицы . Остальные элементы этого массива суть элементы матрицы . На втором шаге, оставляя первые две строки неизменными, мы преобразуем остальные строки, получая нули во втором столбце ниже главной диагонали. Этих нулей ровно столько, сколько множителей определяют второй столбец матрицы . Поэтому мы можем разместить множители на местах этих нулей. При этом вторая строка, начиная со второго элемента, полностью определяет вторую строку матрицы . Продолжая действовать подобным образом, мы сможем разместить полную информацию о сомножителях и в одном двумерном массиве , так что при вычислении -разложения не требуется никакой дополнительной памяти.
§ 3. Метод квадратного корня
Предположим, что все угловые миноры матрицы отличны от нуля. Такая матрица имеет -разложение , где – нижняя треугольная матрица с единицами на главной диагонали, а – невырожденная верхняя треугольная матрица. Невырожденность матрицы означает, что ее диагональные элементы отличны от нуля. В таком случае матрицу можно представить в виде , где
, ,
а , .
Нетрудно проверить, что данное представление определяется единственным образом. -разложение матрицы приобретает теперь вид
. (4.1)
Представление матрицы в виде произведения нижней треугольной матрицы с единицами на главной диагонали, невырожденной диагональной матрицы и верхней треугольной матрицы с единицами на главной диагонали называется ее -разложением. Таким образом, равенство (4.1) является -разложением матрицы . Легко видеть, что -разложение определяется единственным образом.
Рассмотрим -разложение симметричной матрицы . Из равенства следует , или . Так как – нижняя треугольная матрица с единицами на главной диагонали, а – верхняя треугольная матрица с единицами на главной диагонали, то представление также является -разложением матрицы . В силе единственности -разложения , так что
. (4.2)
Предположим теперь, что матрица не только симметрична, но и положительно определена. Согласно критерию Сильвестра все угловые миноры матрицы положительны. Поэтому для симметричной положительно определенной матрицы -разложение вида (4.2) заведомо существует. Покажем, что матрица в этом разложении также положительно определена. Воспользуемся невырожденностью матрицы и в выражении сделаем замену переменной :
.
Очевидно, что из условия вытекает и условие . Но матрица предполагается положительно определенной, следовательно, , а это и означает, что матрица также положительно определена.
В силу критерия Сильвестра все угловые миноры матрицы положительны, откуда следует, что при всех . Поэтому мы можем рассмотреть матрицу
,
которая, очевидно, удовлетворяет равенству .
Введем в рассмотрение еще одну матрицу , которая, очевидно, является невырожденной верхней треугольной матрицей. С использованием этого обозначения разложение (4.2) принимает вид
. (4.3)
Это разложение приводит к одному из наиболее употребительных методов решения систем линейных уравнений с симметричными положительно определенными матрицами – методу квадратного корня. Метод квадратного корня заключается в получении разложения (4.3) и последующем решении двух систем с треугольными матрицами
.
Лекция 8
Глава 5. Итерационные методы решения
систем линейных алгебраических уравнений
§ 1. Двухслойные итерационные методы
В отличие от прямых методов, в которых при отсутствии погрешностей округления точное решение системы находится за конечное число шагов, в итерационных методах всегда определяется некоторая последовательность векторов , которая только в пределе при дает точное решение системы:
,
где – точное решение рассматриваемой системы .
Сходимость последовательности векторов к вектору будем понимать в смысле сходимости по норме:
.
Из курса функционального анализа известно, что в конечномерном пространстве из сходимости последовательности в одной норме вытекает ее сходимость и в любой другой норме. Поэтому для доказательства сходимости достаточно использовать какую-нибудь одну, наиболее удобную, норму.
Из всего многообразия итерационных методов мы рассмотрим только так называемые двухслойные итерационные методы решения системы линейных уравнений , каждый из которых можно записать в следующем каноническом виде
,
где – некоторая невырожденная матрица, – итерационный параметр -ого шага, – произвольное начальное приближение. Если при всех , то метод называется стационарным. Далее мы будем рассматривать только стационарные двухслойные методы
. (5.1)
Предполагая матрицу невырожденной, покажем, что если метод (5.1) сходится, то он непременно сходится к точному решению системы . Действительно, пусть существует предел . переходя в равенстве (5.1) к пределу при , получаем . В силу единственности решения системы отсюда следует, что . Таким образом, выбирая тем или иным образом матрицу и параметр , нам следует позаботиться только о самом факте сходимости последовательности приближений . При этом, конечно, желательно обеспечить наиболее быструю сходимость. Это желание, как правило, вступает в противоречие с другим естественным желанием – простоты вычисления каждого очередного приближения по предыдущему приближению . В самом деле, для вычисления приближения на каждом шаге необходимо решить систему
. (5.2)
Трудоемкость решения этой системы определяется структурой матрицы . Ясно, что наиболее просто решаются системы с треугольными и диагональными матрицами. С другой стороны, можно обеспечить сходимость метода к точному решению за один шаг, выбрав в качестве матрицы матрицу . В этом случае метод фактически оказывается прямым. Но тогда теряет смысл сам итерационный метод – ведь мы отказываемся от применения прямого метода как раз из-за трудоемкости решения исходной системы. На практике в качестве матрицы , как правило, используют треугольную или диагональную матрицу.
Вектор называется вектором погрешности -ого шага, а вектор – вектором невязки -ого приближения.
В терминах векторов погрешностей метод (5.1) принимает вид
,
а условие сходимости:
.
Выразим вектор через вектор :
,
или
,
где . Матрица называется матрицей перехода итерационного метода (5.1). Как мы увидим, свойства сходимости метода полностью определяются его матрицей перехода.
Важной характеристикой итерационного метода является так называемая асимптотическая скорость сходимости. Пусть – наименьшее число шагов, достаточное для уменьшения нормы начальной погрешности в раз, т.е.
. (5.3)
Так как , то
.
Отсюда следует, что неравенство (5.3) будет заведомо выполнено, если потребовать, чтобы имело место неравенство
,
или
. (5.4)
Естественно, что говорить о скорости сходимости имеет смысл, только если метод сходится. В дальнейшем будет установлено, что достаточным условием сходимости итерационного метода является неравенство . Имея в виду это условие, из неравенства (5.4) получаем:
.
Величина и называется асимптотической скоростью сходимости итерационного метода (5.1).
Теорема 1 (критерий сходимости итерационного метода). Для того чтобы итерационный метод (5.1) сходился при любом начальном приближении, необходимо и достаточно, чтобы все собственные значения матрицы перехода этого метода были по модулю меньше 1.
Теорема 2 (достаточное условие сходимости итерационного метода). Для сходимости итерационного метода (5.1) при любом начальном приближении достаточно, чтобы хотя бы одна из норм матрицы перехода, согласованных с какой-либо векторной нормой, была меньше 1.
Доказательство. Пусть – произвольное собственное значение матрицы перехода и – соответствующий ему собственный вектор, т.е. . Если согласована с векторной нормой , то
.
Отсюда, учитывая, что , получаем:
,
т.е. модуль любого собственного значения матрицы перехода не превосходит любой ее согласованной нормы. Если , то в силу теоремы 1 итерационный метод сходится.
§ 2. Примеры двухслойных итерационных методов
1. Метод простой итерации.
В этом методе , а – любое положительное число. Следовательно, он записывается в следующем каноническом виде:
.
В методе простой итерации приближение находится по явным формулам
.
Поставим следующую оптимизационную задачу: подобрать параметр так, чтобы асимптотическая скорость сходимости была наибольшей (т.е. была наименьшей). Мы приведем решение этой задачи в простейшем случае, когда матрица системы является симметричной и положительно определенной. Пусть
.
Тогда решение поставленной задачи есть
.
При этом
.
2. Метод Якоби.
Представим матрицу в виде
,
где – нижняя треугольная матрица, диагональные элементы которой равны 0, а поддиагональные элементы совпадают с соответствующими элементами матрицы ; – верхняя треугольная матрица, диагональные элементы которой равны 0, а наддиагональные элементы совпадают с соответствующими элементами матрицы ; – диагональная матрица, на главной диагонали которой стоят соответствующие элементы матрицы .
Предположим, что матрица невырождена, т.е. . Полагая и , получаем метод Якоби:
,
или
.
Так как – диагональная матрица, то решение находится по явным формулам
.
Матрица перехода метода Якоби имеет вид
.
Приведем одно достаточное условие сходимости метода Якоби, которое формулируется в терминах самой матрицы системы, а не матрицы перехода . Для этого нам понадобится определение. Говорят, что матрица имеет строгое диагональное преобладание, если для всех выполняются неравенства
.
Теорема 3. Если матрица системы имеет строгое диагональное преобладание, то метод Якоби для такой системы сходится при любом начальном приближении.
3. Метод Зейделя.
Положим и . В результате получим метод Зейделя
,
или
.
Поскольку матрица является треугольной, то решение находится по явным формулам
.
Матрица перехода метода Зейделя имеет вид
.
Теорема 4. Если матрица системы имеет строгое диагональное преобладание, то метод Зейделя для такой системы сходится при любом начальном приближении.
Теорема 5. Если матрица системы является симметричной и положительно определенной, то метод Зейделя для такой системы сходится при любом начальном приближении.
Лекция 9
Глава 6. Численное интегрирование
обыкновенных дифференциальных уравнений
§ 1. Метод Эйлера
Рассмотрим задачу Коши для обыкновенного дифференциального уравнения первого порядка
(6.1)
Общее решение уравнения имеет вид
, (6.2)
где – произвольная постоянная. Геометрически общее решение (6.2) представляет собой семейство так называемых интегральных кривых, т.е. совокупность линий, соответствующих различным значениям постоянной . Интегральные кривые обладают тем свойством, что в каждой их точке угловой коэффициент касательной равен значению функции в этой точке. Если задать точку , через которую должна проходить интегральная кривая, то с помощью условия из бесконечного семейства интегральных кривых выделяется некоторая определенная интегральная кривая , соответствующая частному решению уравнения . Это частное решение и есть решение задачи Коши (6.1).
Задачу Коши не всегда удается решить аналитически, т.е. найти ее точное решение возможно не всегда. В таких случаях используют приближенные (численные) методы.
Рассмотрим один из наиболее употребительных численных методов решения задачи Коши (6.1) – метод Эйлера. Разделим отрезок на равных частей точками
,
где , а . В точке к интегральной кривой уравнения проведем касательную
.
Учитывая, что , уравнение касательной перепишем в виде
.
Будем считать, что значение искомого решения в точке приближенно равно значению ординаты точки на касательной, соответствующей абсциссе , т.е.
.
Обозначая это приближенное значение через , приходим к формуле
,
называемой формулой метода Эйлера.
Заметим, что точка лежит уже не на искомой интегральной кривой, а на некоторой другой, близкой к искомой, интегральной кривой уравнения. Исходя из этой точки и повторяя описанные рассуждения, получим приближенное значение решения задачи Коши в следующей точке :
.
В общем случае формула метода Эйлера имеет вид
.
Она позволяет вместо аналитической функции , являющейся точным решением задачи Коши (6.1), найти таблицу ее приближенных значений в точках сетки
.
Известно, что метод Эйлера относится к методам первого порядка точности, т.е. его погрешность аппроксимации имеет порядок .
Рассмотрим несколько более точных методов, являющихся модификациями метода Эйлера.
1. Первый улучшенный метод Эйлера описывается формулами
.
2. Второй улучшенный метод Эйлера описывается формулами
.
Первый и второй улучшенные методы Эйлера имеют погрешность .
§ 2. Методы Рунге-Кутты
Методы Эйлера являются частными случаями методов Рунге-Кутты, к описанию которых мы приступаем.
Явный -этапный метод Рунге-Кутты состоит в следующем. Пусть приближенное значение решения задачи (6.1) в точке уже известно. Задаются числовые коэффициенты
,
,
,
и последовательно вычисляются значения
,
,
.
Очередное приближенное значение решения задачи (6.1) в точке находится по формуле
.
Коэффициенты , , выбираются из соображений точности получаемой формулы.
При и получаем метод Эйлера, рассмотренный в предыдущем параграфе.
При получаем семейство методов
(6.3)
Исследуем погрешность аппроксимации методов (6.3) в зависимости от выбора параметров. Подставляя в последнее равенство выражения для и , получаем:
. (6.4)
По определению погрешностью аппроксимации называется выражение
(6.5)
получаемое заменой в (6.4) приближенных значений и точными значениями и соответственно. Найдем порядок погрешности в предположении, что функции и обладают достаточной степенью гладкости. Для этого разложим все величины, входящие в выражение (6.5) по формуле Тейлора в точке . Опуская достаточно громоздкие выкладки, окончательно получаем:
Отсюда видно, что если , то методы (6.3) имеют не менее чем первый порядок аппроксимации. Если же дополнительно потребовать, чтобы выполнялись условия и , то эти методы будут иметь уже второй порядок аппроксимации.
Таким образом, имеется однопараметрическое семейство двухэтапных методов Рунге-Кутты второго порядка аппроксимации. Это семейство можно записать в виде
, (6.6)
где , или .
В частности, при , получаем первый улучшенный метод Эйлера, рассмотренный в предыдущем параграфе:
При , получаем второй улучшенный метод Эйлера:
Приведем еще несколько примеров методов Рунге-Кутты более высокого порядка аппроксимации.
1. Метод третьего порядка:
,
,
.
2. Метод третьего порядка:
,
,
.
3. Метод четвертого порядка:
,
,
.
4. Метод четвертого порядка:
,
,
.

Методы численного дифференцирования

Формулы на основе разложения функций по формуле Тейлора

Рассмотрим решение задач 1 и 2 численного дифференцирования на различных шаблонах.

А. Двухточечный шаблон. Выберем шаблон [math]H_{2,i}= (x_i,x_{i+1})[/math] на неравномерной сетке [math]Omega_n[/math]. Предполагая, что [math]f(x)in C_2[a,b][/math], разложим функцию [math]f(x)[/math] по формуле Тейлора (В.20) при [math]k=1[/math] относительно точки [math]x_i[/math] с остаточным слагаемым в форме Лагранжа и найдем выражение для [math]f_{i+1}= f(x_{i+1})colon[/math]

[math]f_{i+1}= f_i+ h_{i+1}f’_i+ frac{h_{i+1}^2}{2}f”(xi),qquad scriptstyle{mathsf{(5.3)}}[/math]

где [math]xiin(x_i,x_{i+1}),~ h_i=x_{i+1}-x_i[/math]. Отсюда получаем [math]f’_i= frac{f_{i+1}-f_i}{h_{i+1}}-frac{h_{i+1}}{2}f”(xi)[/math]. Очевидно, справедлива оценка

[math]left|-frac{h_{i+1}}{2}f”(xi)right|leqslant frac{h_{i+1}}{2} max_{xin [x_i,x_{i+1}]} bigl|f”(x)bigr|= frac{h_{i+1}cdot M_{2,i}}{2},,[/math]

где [math]M_{2,i}= max_{xin [x_i,x_{i+1}]} bigl|f”(x)bigr|[/math]. Отсюда следует функциональная формула (функциональный оператор) для первой производной:

В скобках справа от аппроксимационных операторов здесь и далее указываются правые части оценок их погрешностей. Отметим, что формула (5.4) является несимметричной, односторонней (левосторонней). Если функцию [math]f(x)[/math] разложить по формуле Тейлора относительно точки [math]x_{i+1}[/math], то получим правостороннюю формулу

[math]widehat{f},’_{i+1,c}= frac{f_{i+1}-f_i}{h}qquad left(frac{h}{2} M_{2,i}right)!.[/math]

Б. Трехточечный шаблон. На неравномерной сетке [math]Omega_n[/math] выбираем трехточечный (двухшаговый) шаблон [math]H_{3,i}= (x_{i-1}, x_i,x_{i+1})[/math], характеризующийся шагами [math]h_{i+1}= x_{i+1}-x_i,~ h_i=x_i-x_{i-1}[/math] и параметром его нерегулярности [math]delta_{i+1}= h_{i+1}!!not{phantom{|}},h_i[/math] который в общем случае не равен единице.

Аппроксимационные функциональные (точечные) формулы второго порядка в левой крайней, центральной и правой крайней точках шаблона можно получить на основе разложения функции [math]f(x)[/math] по формуле Тейлора с остаточным слагаемым в форме Лагранжа. При этом предполагается, что [math]f(x)in C_3[a,b][/math]. Это позволяет получить разностные дифференциальные операторы [math]widehat{f},'(x_t)~ (t=i-1,i,i+1)[/math] и провести оценки их погрешностей.

Разложим функцию [math]f(x)[/math] при [math]x=x_i[/math] и [math]x=x_{i+1}[/math] по формуле Тейлора при [math]k=2[/math] относительно точки [math]x_{i-1}[/math]с остаточным слагаемым в форме Лагранжа. В результате находим соотношения, определяющие [math]f_i= f(x_i)[/math] и [math]f_{i+1}= f(x_{i+1})colon[/math]

где [math]H_{i}^{i+1}= h_i+h_{i+1},~~ xi_{-}=(x_{i-1},x_i),~~ xi_{+}= (x_{i-1}, x_{i+1}),~~ f_{i-1}^{(p)}(x_{i-1})= f^{(p)}(x_{i-1}),~ p=0,1,2[/math].

Исключая из (5.5), (5.6) слагаемое, содержащее вторую производную, и выражая из полученного соотношения [math]f’_{i-1}[/math], получаем следующую аппроксимацию первой производной в левой крайней точке (левостороннюю формулу или оператор)

При [math]h=text{const}~ (delta_{i+1}=1)[/math] формула (5.7) упрощается и приводится к известному виду:

Данную формулу можно записать через конечные разности: [math]widehat{f},’_{i-1, mathsf{c}}= frac{1}{2h}bigl(3 Delta f_{i-1}-Delta f_ibigr)[/math]. Здесь нижние индексы [math]mathsf{v}[/math] и [math]mathsf{c}[/math], относящиеся к аппроксимационным операторам (5.7) и (5.8), указывают на тип шаблона — нерегулярный [math](h_{i+1}= text{var})[/math] и регулярный [math](h_{i+1}= text{const})[/math]). Остаточное слагаемое для (5.7) получается равным [math]tfrac{1}{6}h_i^2(1+delta_{i+1})f”'(xi),[/math] [math]xiin (x_{i-1}, x_{i+1})[/math] и поэтому для этой аппроксимации справедлива следующая оценка погрешности:

[math]bigl|f’_{i-1}-widehat{f}_{i-1, mathsf{v}}bigr| leqslant frac{1}{6} h_i^2 (1+ delta_{i+1})M_{3,i}[/math], где [math]M_{3,i}= max_{xin H_{3,i}} f”'(x)[/math].

Приводимые здесь и ниже остаточные слагаемые для дифференциальных операторов также получаются путем дифференцирования остаточных слагаемых [math]R(x)[/math] интерполяционных многочленов соответствующей степени. Из (4.20) для трехточечных формул при [math]n=2[/math] следует соотношение

[math]R(x)= frac{}{}f”'(xi)cdot omega(x),qquad omega(x)= (x-x_{i-1})(x-x_i)(x-x_{i+1}).[/math]

Аналогично, разложив функцию [math]f(x)[/math] относительно точки [math]x_{i+1}[/math] и получив соотношения для [math]f_{i-1},,f_i[/math], найдем [math]widehat{f}_{i+1}[/math] — разностный оператор, аппроксимирующий первую производную [math]f’_{i+1}[/math] в правой крайней точке (правосторонняя формула):

Оператор [math]widehat{f}_{i+1, mathsf{v}}[/math] имеет остаточное слагаемое [math]frac{h_i^2}{6} delta_{i+1}(1+delta_{i+1})f”'(xi)[/math].

Разложение функции [math]f(x)[/math] относительно центральной точки [math]x_i[/math] шаблона, получение выражений для [math]f_{i-1},,f_{i+1}[/math] и исключение из них слагаемого со второй производной приводят к следующим разностным операторам функционального типа, аппроксимирующим первую производную в центральной точке (формула центрального вида)

Оператор [math]widehat{f}_{i, mathsf{v}}[/math] имеет остаточное слагаемое [math]frac{h_i^2}{6} delta_{i+1}f”'(xi)[/math].

Приведенные остаточные слагаемые разностных операторов обусловливают следующие оценки их погрешностей:

Замечания

1. Далее в тексте оценочная константа [math]M_{p,i}[/math] для краткости будет использоваться без дополнительного ее описания. В нижнем индексе этой константы всюду указывается [math]p[/math] — порядок производной.

2. Из оценок (5.11) вытекает, что разностные операторы [math]widehat{f}_{i-1},, widehat{f}_{i},, widehat{f}_{i+1}[/math] аппроксимируют при [math]h_{i+1}=text{var}[/math] соответствующие производные [math]f’_{i-1},, f’_i,, f’_{i+1}[/math] со вторым порядком, если шаблон произвольный (безусловная аппроксимация). Если же на шаблоне с [math]delta_{i+ 1}ll1[/math] наложить условие, например [math]delta_{i+1}leqslant h_i[/math], то есть [math]h_{i+1}leqslant h_i^2[/math], то порядок аппроксимации [math]widehat{f}_{i},, widehat{f}_{i+1}[/math] может быть повышен до третьего (условная аппроксимация). Такая возможность повышения порядка аппроксимации относительно [math]h_i[/math] без увеличения количества точек шаблона обеспечивается введением в мажоранты, соответствующие аппроксимациям, параметра [math]delta_{i+1}[/math], на который в случае необходимости можно наложить условие [math]delta_{i+1}leqslant h_i[/math]. При этом следует иметь в виду, что данный параметр входит в знаменатель некоторых слагаемых аппроксимационных формул (5.9), (5.10) и при его уменьшении увеличиваются погрешности арифметических операций.

3. Аппроксимации (5.8),(5.9) являются условными, так как для них справедливо условие [math]delta_{i+1}=1[/math]. Предположим, что [math]f(x)in C_4[a,b][/math], и разложим функцию [math]f(x)[/math] в точках [math]x_i,,x_{i+1}[/math] на трехточечном нерегулярном шаблоне [math]H_{3,i}= (x_{i-1},x_i,x_{i+1})[/math] до слагаемого четвертого порядка относительно шага. Складывая эти разложения и выражая из суммы вторую производную, получаем функционально-дифференциальную формулу для второй производной:

Подставляя в (5.12) аппроксимационную формулу (5.10) для первой производной, находим разностный аппроксимационный оператор [math]widehat{f},”_i[/math], выраженный через параметры [math]delta_{i+1},, h_i^2[/math] и аппроксимирующий (безусловно) вторую производную [math]f”_i[/math] на нерегулярном шаблоне с первым порядком:

Выражение (5.13) можно преобразовать к виду

[math]widehat{f},”_{i,mathsf{v}}= frac{2}{H_{i}^{i+1}}! left(frac{Delta f_i}{h_{i+1}}-frac{Delta f_{i-1}}{h_i}right)= frac{2}{H_{i}^{i+1}}! left[frac{f_{i-1}}{h_i}-left(frac{1}{h_i}+ frac{1}{h_{i+1}}right)!f_i+ frac{f_{i+1}}{h_{i+1}}right]!.[/math]

Если сетка равномерная [math](delta_{i+1}=1)[/math], то указанный порядок условной аппроксимации возрастает на единицу, так как третье слагаемое в (5.12) становится равным нулю. В этом случае из (5.13) получаем широко распространенный оператор, аппроксимирующий вторую производную на регулярном шаблоне:

Замечание. Из (5.12) следует, что порядок условной аппроксимации (5.13) можно повысить на единицу и на нерегулярном шаблоне, если принять [math]|delta_{i+1}-1|<h_i[/math], то есть [math]h_i-h_i^2 leqslant h_{i+1}leqslant h_i+h_i^2[/math] (квазиравномерная сетка).

Как следует из приведенных выше постановок задач, в вычислительной практике аппроксимационные формулы (операторы) для производных используются для вычисления значений производных либо для замены ими соответствующих дифференциальных операторов. В связи с этим приведем методику вычисления [math]Bigl.{f^{(p)}(x)}Bigr|_{x=x_j}[/math].

---

Методика вычисления значений производных

1. Выбрать конкретную аппроксимационную формулу (или несколько разных формул), в которой порядок аппроксимации должен соответствовать заданному в задаче порядку точности [math]t[/math].

2. Выбрать наборы точек (шаблоны [math]H_{k,i}[/math]), которым принадлежат точки [math]x_j~(x_jin H_{k,i})[/math], причем для каждого из наборов расположение точки [math]x_j[/math] должно быть зафиксировано относительно точек шаблона ([math]k[/math] — количество точек, определяющих шаблон). Эта фиксация определяется структурой формулы. Например, если формула имеет центральный тип (см. формулу (5.10)), то точка [math]x_j[/math] должна совпадать со средней точкой шаблона, а если формула левосторонняя (см. формулы (5.7) и (5.8)), то точка [math]x_j[/math] должна совпадать с левой крайней точкой шаблона и т.д.

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

4. Произвести требуемые вычисления с учетом того, что количество сохраняемых цифр должно приблизительно соответствовать величине остаточного слагаемого аппроксимационной формулы и порядку точности [math]t[/math].

---

Оценка погрешности аппроксимационного оператора

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

Данную процедуру можно выполнить с помощью различных подходов, один из которых основан на разложении функций, входящих в правую часть оператора, по формуле Тейлора относительно той точки [math]x_j[/math], для которой записан этот оператор. Другой подход использует анализ остаточного слагаемого, полученного для интерполяционного многочлена Лагранжа. При этом рассматривается соотношение (4.19), которое дифференцируется нужное количество раз:

[math]f^{(p)}(x)= L_{n}^{(p)}+ R_{n}^{(p)}(x).[/math]

Тогда если нужно найти погрешность численного дифференцирования в точке [math]x=x_j[/math], то осуществляется подстановка [math]x=x_j[/math]. В результате находится оценка погрешности в точке

[math]left|Bigl.{R_n^{(p)}(x)}Bigr|_{x=x_j}right|leqslant frac{M_{n+1}}{(n+1)!}cdot omega_n^{(p)}(x_j),[/math]

где [math]n[/math] — степень алгебраического многочлена, по которому получен аппроксимационный оператор. Если априори степень [math]n[/math] неизвестна, то она может быть определена путем подстановки в правую часть оператора произвольных узлов [math]x_i[/math] и значений [math]f(x_i)[/math] для многочленов [math]N_1(x)= ax+b;~ N_2(x)= ax^2+bx+c[/math] и т. д. Максимальная степень многочлена [math]N_n(x)[/math], для которого остаточное слагаемое равно нулю, является искомой.

Если необходимо оценить погрешность аппроксимационного оператора [math]widehat{f}^{(p)} (x)[/math] не в точке, а на всем отрезке [math][a,b][/math], то для этого используется неравенство (4.22), которое дифференцируется [math]p[/math] раз:

[math]max_{xin[a,b]} bigl|R_n^{(p)}(x)bigr|leqslant frac{M_{n+1}}{(n+1)!}cdot max_{xin [a,b]} bigl|omega_n^{(p)}(x)bigr|.[/math]

Изложим простейшую методику оценки погрешности в точке [math]x=x_j[/math], когда степень многочлена [math]L_{k-1}(x)[/math] и шаблон [math]H_{k,i}[/math], на котором этот многочлен получен, известны.

---

Методика оценки погрешности аппроксимационного оператора

1. На заданном шаблоне, который в общем случае имеет структуру [math]H_{k,i}= (x_{i-r}, ldots, x_i,ldots, x_{i+s}),~ k=r+s+1[/math], записать выражение для оценки остаточного слагаемого в точке [math]x=x_jin H_{k,i}colon[/math]

[math]left|Bigl.{R_{k-1}^{(p)}(x)}Bigr|_{x=x_j}right|leqslant frac{M_k}{k!} left|Bigl.{omega_{k-1}^{(p)}(x)}Bigr|_{x=x_j}right|.[/math]

Здесь [math]r[/math] и [math]s[/math] — количество точек, расположенных соответственно левее и правее точки [math]x_j[/math], которая фиксируется на этом шаблоне.

2. Найти производную [math]omega_{k-1}^{(p)}(x)= bigl[(x-x_{i-r})(x-x_{i-r+1})ldots (x-x_{i+s})bigr]^{(p)}[/math].

3. В полученную производную подставить значение [math]x_j[/math]. Далее преобразовать ее, выразив через [math]h[/math] (при [math]h=text{const}[/math]) или через параметр нерегулярности [math]delta_{i+1}= h_{i+1}!!not{phantom{|}},h_i[/math], и записать окончательно оценку погрешности.

---

Четырехточечный шаблон

Формулы третьего порядка для первых производных на регулярном шаблоне [math]H_{4,i}= (x_{i-2},x_{i-1},x_i,x_{i+1})[/math] имеют вид:

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

Формулы второго порядка на регулярном шаблоне для вторых производных имеют вид

---

Пятиточечный шаблон

Формулы четвертого порядка для первых производных на регулярном шаблоне [math]H_{5,i}= (x_{i-2}, x_{i-1},x_i,x_{i+1},x_{i+2})[/math] имеют вид

Формулы третьего порядка для вторых производных на указанном шаблоне имеют вид

Замечание. Формулы, записанные выше для аппроксимации производных через приращения функций, могут быть переписаны для восстановления функции [math]y=f(x)[/math] по значениям интегралов. С этой целью можно использовать изложенный ранее метод подобия, состоящий в соответствующем изменении порядка производной в левой и правой частях аппроксимационного выражения. Так, из оператора [math]widehat{f},’_k[/math] ([math]k[/math] — номер точки шаблона) можно получить оператор [math]widehat{f},’_k[/math] путем замены [math]Delta f_{k-1}= f_k-f_{k-1}[/math] на интеграл [math]I_{k-1}^k= F_k-F_{k-1}[/math]. Проделав это, вместо формулы, следующей за (5.8), вторых из формул (5.9),(5.10) и формул (5.15) получим операторы, восстанавливающие функцию [math]y=f(x)[/math] по значениям интегралов:

[math]begin{aligned}&widehat{f}_{i-1,c}= frac{1}{2h} (3I_{i-1}^{i}-I_{i}^{i+1}); &qquad & widehat{f}_{i,c}= frac{1}{2h} (I_{i}^{i+1}+ I_{i-1}^{i});\ &widehat{f}_{i-2,c}= frac{1}{6h} (11I_{i-2}^{i-1}-7I_{i-1}^{i}+ 2I_{i}^{i+1}); &qquad & widehat{f}_{i-1,c}= frac{1}{6h} (2I_{i-2}^{i-1}+ 5I_{i-1}^{i}-I_{i}^{i+1});\ &widehat{f}_{i,c}= frac{1}{6h} (-I_{i-2}^{i-1}+ 5I_{i-1}^{i}+ 2I_{i}^{i+1}); &qquad & widehat{f}_{i+1,c}= frac{1}{6h} (2I_{i-2}^{i-1}-7I_{i-1}^{i}+ 11I_{i}^{i+1}). end{aligned}[/math]

Первые две формулы имеют второй порядок аппроксимации по [math]h[/math] (они записаны на трехточечном шаблоне), а последние четыре — третий порядок (они записаны на четырехточечном шаблоне).

Вместо формулы (5.14) имеем следующую: [math]widehat{f}_{i,c}= frac{1}{2h}(I_{i}^{i+1}-I_{i-1}^{i})[/math]. Аналогичным путем можно найти и другие интегральные аппроксимации первой производной.

---

Формулы на основе разложения первообразных по формуле Тейлора

А. Двухточечный шаблон. На неравномерной сетке [math]Omega_n[/math] рассмотрим двухточечный шаблон [math]H_{2,i}= (x_{i},x_{i+1})[/math] с шагом [math]h_{i+1}= x_{i+1}-x_{i}[/math]. Предположив, что [math]f(x)in C_2[a,b][/math], разложим первообразную [math]F(x)[/math] относительно точки [math]x_i[/math] по формуле Тейлора при [math]k=2[/math] и найдем выражение для [math]Bigl.{F(x)}Bigr|_{x=x_{i+1}}= F_{i+1}colon[/math]

[math]F_{i+1}= F_{i}+ h_{i+1}cdot f_{i}+ frac{h_{i+1}^2}{2}cdot f’_{i}+ frac{h_{i+1}^3}{6}cdot f”(xi),[/math]

где [math]xiin (x_{i},x_{i+1})[/math]. Выражая из этого разложения сначала первую производную, а затем интеграл [math]I_{i}^{i+1}= F_{i+1}-F_{i}[/math], получаем

[math]f’_{i}= frac{2I_{i}^{i+1}}{h_{i+1}^2}-frac{2f_{i}}{h_{i+1}}-frac{h_{i+1}}{3}cdot f”(xi),qquad I_{i}^{i+1}= h_{i+1}cdot f_{i}+ frac{h_{i+1}^2}{2}cdot f’_{i}+ frac{h_{i+1}^3}{6}cdot f”(xi).[/math]

Первые два слагаемых в правых частях последних двух соотношений представляют собой интегрально-функциональную (интегрально-точечную) аппроксимацию производной [math]f’_{i}[/math] и функционально-дифференциальную аппроксимацию интеграла [math]I_{i}^{i+1}[/math] соответственно:

Порядки этих аппроксимаций устанавливают остаточные слагаемые в выражениях для [math]f’_{i}[/math] и [math]I_{i}^{i+1}[/math], из которых следуют оценки, правые части которых указаны в скобках рядом с формулами.

Данные оценки свидетельствуют о том, что порядок аппроксимации, обеспечиваемый (5.19), равен единице, а обеспечиваемый (5.20) равен трем при условии, что [math]I_{i}^{i+1}[/math] и [math]f_i[/math] для (5.19) известны с точностью не ниже [math]O(h_{i+1}^3)[/math] и [math]O(h_{i+1}^2)[/math], а [math]f_{i}[/math] и [math]f’_{i}[/math] для (5.20) с точностью не ниже [math]O(h_{i+1}^2)[/math] и [math]O(h_{i+1})[/math]. Этот же результат следует из рассмотренного ранее правила соответствия порядков аппроксимации математических моделей различного типа.

Замечания

1. Выразив из разложения для [math]F_{i+1}[/math] непосредственно функцию [math]f_{i}[/math], получим еще одно аппроксимационное выражение: [math]widehat{f}_{i}= frac{1}{h_{i+1}}cdot I_{i}^{i+1}-frac{h_{i+1}}{2}cdot f’_{i},[/math] позволяющее восстанавливать со вторым порядком аппроксимации значение функции [math]f_{i}[/math] в точке [math]x_{i}[/math] по значениям [math]I_[i}^{i+1}[/math] и [math]f’_{i}[/math], известным с точностью не ниже [math]O(h_{i+1}^3)[/math] и [math]O(h_{i+1})[/math] соответственно.

2. Если в точке [math]x_{i}[/math] известны производные более высоких порядков, то могут быть построены аппроксимации [math]widehat{f},’_{i},, widehat{I}_{i}^{i+1},, widehat{f}_{i}[/math] более высокого порядка, которые здесь не рассматриваются.

Б. Трехточечный шаблон. Возьмем на неравномерной сетке [math]Omega_{n}[/math] трехточечный шаблон [math]H_{3,i}= (x_{i-1},x_{i},x_{i+1})[/math], характеризующийся шагами [math]h_{i+1}= x_{i+1}-x_{i},~ h_{i}= x_{i}-x_{i-1}[/math] и параметром нерегулярности [math]Delta_{i+1}= frac{h_{i+1}}{h_{i}}[/math]. Предположив, что [math]f(x)in C_3[a,b][/math], найдем разложения для первообразной [math]F(x)[/math] при [math]x=x_{i+1}[/math] и [math]x=x_{i-1}[/math] относительно точки [math]x_{i}[/math] и выпишем выражения для [math]F(x_{i+1}),, F(x_{i+1})[/math] по формуле Тейлора (В.23) при [math]k=3colon[/math]

[math]begin{gathered}F_{i+1}= F_{i}+ h_{i+1}cdot f_{i}+ frac{h_{i+1}^2}{2}cdot f’_{i}+ frac{h_{i+1}^3}{6}cdot f”_{i}+ frac{h_{i+1}^4}{24}cdot f”'(xi_1),quad xi_1in (x_{i},x_{i+1}),\ F_{i-1}= F_{i}- h_{i}cdot f_{i}+ frac{h_{i}^2}{2}cdot f’_{i}-frac{h_{i}^3}{6}cdot f”_{i}+ frac{h_{i}^4}{24}cdot f”'(xi_2),quad xi_2in (x_{i-1},x_{i}).end{gathered}[/math]

Умножая первое соотношение на [math]h_{i}[/math] а второе на [math]h_{i+1}[/math], складывая их с учетом равенства [math]I_{i}^{i+1}= F_{i+1}-F_{i}[/math] и разрешая относительно [math]f’_{i}[/math], получаем

[math]f’_{i}= frac{2}{h_{i}(1+delta_{i+1})}! left(frac{I_{i}^{i+1}}{h_{i+1}}-frac{I_{i-1}^{i}}{h_{i}}right)-frac{h_{i}(delta_{i+1}-1)}{3}f”_{i}-frac{h_{i}^2}{12(1+ delta_{i+1})}bigl(delta_{i+1}^3 f”'(xi_{i})+ f”'(xi_2)bigr).[/math]

Отсюда следует интегрально-дифференциальная аппроксимационная формула для первой производной [math]f’_{i}[/math] на нерегулярном шаблоне:

Формула (5.21) при заданных интегралах [math]I_{i}^{i+1},, I_{i-1}^{i}[/math] (известных с точностью не ниже [math]O(h_{i+1}^3)[/math] и [math]O(h_{i}^3)[/math] и производной [math]f”_{i}[/math] (известной с точностью не ниже первого порядка) аппроксимирует первую производную [math]f’_{i}[/math] на нерегулярном шаблоне со вторым порядком (безусловная аппроксимация).

При условной аппроксимации, когда [math]delta_{i+1}=1[/math] (равномерная сетка), из (5.21) следует интегральная аппроксимационная формула для первой производной [math]f’_{i}[/math] на регулярном шаблоне:

Эта формула аппроксимирует первую производную [math]f’_{i}[/math] со вторым порядком.

Из сопоставления мажорант оценок аппроксимационных операторов [math]widehat{f}_{i,c}[/math] (см. (5.10)) и [math]widehat{f}_{i,c}[/math] (см. (5.22)) вытекает, что они имеют одинаковый (второй) порядок аппроксимации, однако, мажоранта или остаточное слагаемое оператора интегрального типа содержит константу [math](1!!not{phantom{|}},12)[/math], в два раза меньшую соответствующей константы в мажоранте оператора точечного (функционального) типа.

Замечание. При условии [math]f(x)in C_2[a,b][/math] на регулярном шаблоне из разложения первообразных получается аппроксимационная интегрально-функциональная формула для второй производной:

из которой можно выразить интеграл [math]I_{i-1}^{i+1}[/math] через значения функции [math]f_{i}[/math] и производной [math]f”_{i}colon, widehat{I}_{i-1}^{i+1}= 2hf_{i}+ frac{h^3}{3}f”_{i}[/math].

---

Оценка погрешности дифференциальных операторов

Изложенная ранее методика определения остаточного слагаемого основана на соответствующей теореме для интерполяционного многочлена Лагранжа. Полезно рассмотреть и другой способ получения оценки погрешности, использующей разложение функций по формуле Тейлора.

Методика оценки погрешности дифференциальных операторов на основе разложения функции по формуле Тейлора

1. Составить разность [math]f_{i}^{(p)}-widehat{f}_{i}^{(p)}[/math], где [math]widehat{f}_{i}^{(p)}[/math] — обозначение оператора, в качестве которого может быть принят оператор как функционального типа, так и интегрального или интегрально-функционального, а [math]f_{i}^{(p)}[/math] — точное значение производной.

2. Все функции, входящие в [math]widehat{f}_{i}^{(p)}[/math] разложить по формуле Тейлора с остаточным слагаемым в форме Лагранжа относительно точки [math]x_{i}[/math], в которой записан [math]widehat{f}_{i}^{(p)}[/math]. Если оператор [math]widehat{f}_{i}^{(p)}[/math] выражается через интегралы, то они записываются через разности первообразных [math]bigl(F(x)bigr)[/math]. При этом предполагается, что функции [math]f(x),, F(x)[/math] являются непрерывными, и для них существуют непрерывные производные соответствующего порядка.

3. Подставить разложения функции в разность (см. п. 1) и выполнить ее преобразование, в процессе которого получается остаточное слагаемое.

4. На основе остаточного слагаемого оценить погрешность аппроксимации.

---

Формулы, полученные на основе сплайнов

А. Двухточечный шаблон. Из рассмотрения интегрально-дифференциальных сплайн-многочленов второй степени в работе получены интегрально-функциональные аппроксимационные формулы для первых производных на двухточечном шаблоне [math]H_{2,i}= (x_{i},x_{i+1})colon[/math]

Формулы (5.24) в силу их одноинтервального характера могут использоваться при [math]h_{i}= text{const}[/math], если значение интеграла [math]I_{i-1}^{i}[/math] либо известно, либо заранее вычислено с точностью не ниже [math]O(h_{i}^4)[/math]. Для операторов [math]widehat{f},’_{i-1},, widehat{f},’_{i}[/math] справедливы оценки

[math]bigl|widehat{f},’_{k}-f’_{k}bigr|leqslant frac{h_{i}^2}{12}cdot M_{3,i},quad k=i-1,i.[/math]

Замечание. Из (5.24) можно выразить величину [math]I_{i-1}^{i}[/math], и тогда получаются следующие функционально-дифференциальные квадратурные формулы:

[math]begin{aligned}I_{i-1}^{i}&= frac{h_{i}}{3}cdot (2f_{i-1}+ f_{i})+ frac{h_{i}^2}{6}cdot f’_{i-1},\ I_{i-1}^{i}= frac{h_{i}}{3}cdot (f_{i-1}+ 2f_{i})-frac{h_{i}^2}{6}cdot f’_{i}. end{aligned}[/math]

Б. Трехточечный шаблон. Из рассмотрения интегрально-дифференциальных сплайн-многочленов третьей степени на трехточечном шаблоне [math]H_{3,i}= (x_{i-1},x_{i}, x_{i+1})[/math] получены следующие интегрально-функциональные формулы:

– для левой, центральной и правой точек нерегулярного шаблона аппрокси-мационные формулы для первой производной имеют вид

[math]begin{aligned}widehat{f},’_{i-1,nu}&= frac{1}{H_{i}^{i+1}}left[4! left(frac{h_{i}}{h_{i+1}^2}I_{i}^{i+1}+ frac{2H_{i}^{i+1}+h_{i}}{h_{i}^2}I_{i-1}^{i}right)-left(frac{h_{i}}{h_{i+1}}f_{i+1}+ frac{3(H_{i}^{i+1})^2}{h_{i}h_{i+1}}f_{i}+ frac{5H_{i}^{i+1}+h_{i}}{h_{i}}f_{i-1}right)right]!,\[4pt]
widehat{f},’_{i,nu}&= frac{1}{H_{i}^{i+1}}left[4! left(frac{h_{i}}{h_{i+1}^2}I_{i}^{i+1}-frac{h_{i+1}}{h_{i}^2}I_{i-1}^{i}right)-left(frac{h_{i}}{h_{i+1}}f_{i+1}+ frac{3(h_{i}^2-h_{i+1}^2)}{h_{i}h_{i+1}}f_{i}-frac{h_{i+1}}{h_{i}}f_{i-1}right)right]!,\[4pt] widehat{f},’_{i+1,nu}&= frac{1}{H_{i}^{i+1}}left[-4! left(frac{2H_{i}^{i+1}+ h_{i+1}}{h_{i+1}^2}I_{i}^{i+1}+ frac{h_{i+1}}{h_{i}^2}I_{i-1}^{i}right)+ left(frac{5H_{i}^{i+1}+ h_{i+1}}{h_{i+1}}f_{i+1}+ frac{3(H_{i}^{i+1})^2}{h_{i}h_{i+1}}f_{i}+ frac{h_{i+1}}{h_{i}}f_{i-1}right)right]!;end{aligned}[/math]

– для левой, центральной и правой точек регулярного шаблона аппрокси-мационные формулы для первой производной принимают форму

[math]begin{aligned}widehat{f},’_{i-1,c}&= frac{1}{2h}! left[frac{4}{h}(I_{i}^{i+1}+ 5I_{i-1}^{i})-(f_{i+1}+ 12f_{i}+ 11f_{i-1})right] &quad & left(frac{h^3}{60}M_{4,i}right)!,\[4pt] widehat{f},’_{i,c}&= frac{1}{2h}! left[frac{4}{h}(I_{i}^{i+1}-I_{i-1}^{i})-(f_{i+1}-f_{i-1})right]= frac{1}{2h}! left[frac{4}{h} Delta I_{i}^{i+1}-(Delta f_{i}+ Delta f_{i-1})right] &quad & left(frac{h^4}{360}M_{5,i}right)!,\[4pt] widehat{f},’_{i+1,c}&= frac{1}{2h}! left[-frac{4}{h}(5I_{i}^{i+1}+ I_{i-1}^{i})+ (11f_{i+1}+ 12f_{i}+ f_{i-1})right] &quad & left(frac{h^3}{60}M_{4,i}right)!; end{aligned}[/math]

– для левой, центральной и правой точек нерегулярного шаблона аппрокси-мационные формулы для второй производной имеют вид

[math]begin{aligned}widehat{f},”_{i-1,nu}&= frac{6}{H_{i}^{i+1}}! left[-4! left(frac{1}{h_{i+1}^2}I_{i}^{i+1}+ frac{H_{2i}^{i+1}}{h_{i}^3}I_{i-1}^{i}right)+ frac{1}{h_{i+1}}f_{i+1}+ frac{H_{i}^{i+1}(2H_{i}^{i+1}+ h_{i})}{h_{i}^2h_{i+1}}f_{i}+ frac{2H_{i}^{i+1}+ h_{i}}{h_{i}^2}f_{i-1}right]!,\[4pt] widehat{f},”_{i,nu}&= frac{6}{H_{i}^{i+1}}! left[4! left(frac{1}{h_{i+1}^2}I_{i}^{i+1}+ frac{1}{h_{i}^2}I_{i-1}^{i}right)-left(frac{1}{h_{i+1}}f_{i+1}+ frac{3H_{i}^{i+1}}{h_{i+1}h_{i}}f_{i}+ frac{1}{h_{i}}f_{i-1}right)right]!,\[4pt] widehat{f},”_{i+1,nu}&= frac{6}{H_{i}^{i+1}}! left[-4! left(frac{H_{i}^{2(i+1)}}{h_{i+1}^3}I_{i}^{i+1}+ frac{1}{h_{i}^2}I_{i-1}^{i}right)+ frac{2H_{i}^{i+1}+ h_{i+1}}{h_{i+1}^2}f_{i+1}+ frac{H_{i}^{i+1}(2H_{i}^{i+1}+ h_{i+1})}{h_{i}h_{i+1}^2}f_{i}+ frac{f_{i-1}}{h_{i}}right]!; end{aligned}[/math]

– для левой, центральной и правой точек регулярного шаблона аппроксима-ционные формулы для второй производной принимают форму

[math]begin{aligned}widehat{f},”_{i-1,c}&= frac{3}{h^2}! left[-frac{4}{h}(I_{i}^{i+1}+ 3I_{i-1}^{i})+ (f_{i+1}+ 10f_{i}+ 5f_{i-1})right] &quad & left(frac{7h^2}{20}M_{4,i}right)!,\[4pt] widehat{f},”_{i,c}&= frac{3}{h^2}! left[frac{4}{h}(I_{i}^{i+1}+ I_{i-1}^{i})+ (f_{i+1}+ 6f_{i}+ f_{i-1})right] &quad & left(frac{h^2}{20}M_{4,i}right)!,\[4pt] widehat{f},”_{i+1,c}&= frac{3}{h^2}! left[-frac{4}{h}(3I_{i}^{i+1}+ I_{i-1}^{i})+ (5f_{i+1}+ 10f_{i}+ f_{i-1})right] &quad & left(frac{7h^2}{20}M_{4,i}right)!. end{aligned}[/math]

Подчеркнем, что правые части данных аппроксимационных соотношений записаны через интегралы [math]I_{i-1}^{i},, I_{i}^{i+1}[/math] на двух смежных отрезках, составляющих шаблон, и через значения функций в точках этого шаблона. Если известны интегралы и значения самой функции, то по этим формулам можно вычислить первые производные с третьим порядком, а вторые производные — со вторым (остаточные слагаемые некоторых аппроксимационных формул указаны в скобках, расположенных рядом с этими формулами). При этом формула для [math]widehat{f},’_{i,c}[/math], имеющая симметричный вид относительно центральной точки х, и содержащая интегральные и функциональные разности, обеспечивает повышенный (четвертый) порядок аппроксимации. Если интегралы для исследуемой функции неизвестны, они должны быть предварительно рассчитаны с порядком, по крайней мере на два превышающим порядок аппроксимации дифференциальных операторов.

Замечание. В вычислительной практике могут оказаться полезными еще два аппроксимационных оператора [math]widehat{f},”_{L(Pi),nu},, widehat{f},”_{i,c}colon[/math]

где [math]widehat{f},”_{L(Pi),nu}[/math] — лево- или правосторонний оператор, а [math]widehat{f},”_{i,c}[/math] — центральный оператор, записанный на трехточечном шаблоне при [math]h=text{const}[/math].

В. Четырехточечный шаблон. Выше приведены формулы численного дифференцирования на трехточечном (или для некоторых формул на двухточечном) шаблоне, имеющие порядок аппроксимации [math]h^{3-(p-1)}[/math], где [math]p[/math] — порядок производных, для которых записаны эти формулы. В дополнение к изложенному материалу приведем формулы, аппроксимирующие производную [math]f^{(p)}~(p=1,,2)[/math] с порядками [math]h^{3-(p-1)}[/math] на четырехточечном шаблоне. Данные формулы получены путем анализа кубических дифференциальных и интегрально-дифференциальных сплайнов.

1. Формулы для первых производных третьего порядка аппроксимации на шаблоне [math]H_{4,i}= (x_{i-2},x_{i-1},x_{i},x_{i+1})colon[/math]

– для лево- и правосторонних внутренних точек [math]x_{i-1},,x_{i}[/math] нерегулярного шаблона справедливы функциональные формулы:

где [math]a= H_{i-1}^{i+1}H_{i-1}^{i}H_{i}^{i+1}Pi_{i-1}^{i+1};~~ K_{2i}^2= h_{i-1}(H_{i-1}^{i+1}+2h_{i})+ h_{i}(2H_{i}^{i+1}+ h_{i})[/math];
[math]K_{2i-1}^2= h_{i+1}(H_{i-1}^{i+1}+ 2h_{i})+ h_{i}(2H_{i-1}^{i}+h_{i});~~ Pi_{i-1}^{i+1}= h_{i-1}h_{i}h_{i+1};~~ H_{i-1}^{i+1}= h_{i-1}+ h_{i}+ h_{i+1};[/math]

– для левой и правой крайних точек [math]x_{i-2},,x_{i+1}[/math] нерегулярного шаблона справедливы функционально-дифференциальные формулы рекуррентного типа:

На регулярном четырехточечном шаблоне при [math]h=text{const}[/math] формулы (5.27) – (5.30) упрощаются и сводятся к формулам (5.15). Можно показать, что (5.27) – (5.30), так же как и формулы (5.15), имеют третий порядок аппроксимации.

Параболические и кубические интегрально-дифференциальные сплайны позволяют сконструировать интегральные формулы численного дифференцирования.

На четырехточечном нерегулярном шаблоне [math]H_{4,i+1}= (x_{i-1},x_{i},x_{i+1},x_{i+2})[/math] получаются следующие интегральные и рекуррентные формулы для первых производных:

где [math]A= h_{i+1}^2 bigl(2h_{i}+ h_{i+1}+ 2h_{i+2}bigr)+ h_{i+1} bigl(h_{i}^2+ h_{i+2}^2bigr)+ h_{i}h_{i+2} bigl(h_{i}+ 3h_{i+1}+ h_{i+2}bigr)[/math].

Формулы (5.31),(5.32) относятся к внутренним точкам шаблона и при [math]h=text{const}[/math] переходят в (5.22), а формулы (5.33) (рекуррентные) — к крайним точкам шаблона и являются подобными формулам (5.36), справедливыми для вторых производных.

На регулярном шаблоне из (5.33) при [math]h=text{const}[/math] легко получаются явные трехинтервальные аппроксимации интегрального типа:

[math]begin{aligned}widehat{f},’_{i-1,c}&= frac{1}{h^2}bigl(-2I_{i-1}^{i}+ 3I_{i}^{i+1}-I_{i+1}^{i+2}bigr) &quad & left(frac{11}{12}h^2M_{3,i}right)!,\[4pt] widehat{f},’_{i+2,c}&= frac{1}{h^2}bigl(I_{i-1}^{i}-3I_{i}^{i+1}+ 2I_{i+1}^{i+2}bigr) &quad & left(frac{11}{12}h^2M_{3,i}right)!. end{aligned}[/math]

Последние две формулы могут быть получены также методом подобия из первой и последней формул (5.16). Их можно записать через интегральные разности, где [math]Delta I_{k}^{k+1}= I_{k}^{k+1}-I_{k-1}^{k}[/math]:

[math]begin{aligned}widehat{f},’_{i-1,c}&- frac{1}{h^2}bigl(2 Delta I_{i}^{i+1}-Delta I_{i+1}^{i+2}bigr);\[4pt] widehat{f},’_{i+2,c}&- frac{1}{h^2}bigl(2 Delta I_{i+1}^{i+2}-Delta I_{i}^{i+1}bigr). end{aligned}[/math]

2. Формулы для вторых производных второго порядка апппроксимации на шаблоне [math]H_{4,i}= (x_{i-2},x_{i-1},x_{i},x_{i+1})[/math]

– для лево- и правосторонних внутренних точек [math]x_{i-1},,x_{i}[/math] нерегулярного шаблона справедливы следующие функциональные формулы.

где [math]K_{t}^{t+1}= Pi_{t}^{t+1}H_{t}^{t+1};~~ Pi_{t}^{t+1}= h_{t}h_{t+1};~~t=i-1,I;~~ H_{i-1}^{2i}= h_{i-1}+2h_{i}[/math];
[math]H_{2i}^{i+1}= 2h_{i}+ h_{i+1};~~ Delta h_{i+1}= h_{i+1}-h_{i};~~ a=H_{i-1}^{i+1} H_{i-1}^{i} H_{i}^{i+1} Pi_{i-1}^{i+1};[/math]

– для левой и правой крайних точек нерегулярного шаблона справедливы дифференциальные формулы рекуррентного типа:

Формулы (5.27), (5.28) и (5.34), (5.35) могут использоваться для расчета значений производных во внутренних точках сетки [math]{x_{i}},~ i=overline{1,n-1}[/math], а (5.29), (5.30) и (5.36) — для расчета производных в крайних точках [math]x_{0},,x_{n}[/math] сетки [math]Omega_{n}[/math].

На регулярном шаблоне при [math]h=text{const}[/math] последняя группа формул для аппроксимации вторых производных во внутренних точках шаблона преобразуется к традиционным:

[math]begin{aligned}widehat{f},”_{i,c}&= frac{1}{h^2}bigl(f_{i-1}-2f_{i}+ f_{i+1}bigr) quad left(frac{h^2}{12}M_{4,i}right)!;\ widehat{f},”_{i-1,c}&= frac{1}{h^2}bigl(f_{i-2}-f_{i-1}+ f_{i}bigr),end{aligned}[/math]

а (5.36) — к рекуррентным формулам: [math]widehat{f},”_{i-2,c}= 2widehat{f},”_{i-1,c}-widehat{f},”_{i,c};~ widehat{f},”_{i+1,c}= 2widehat{f},”_{i,c}-widehat{f},”_{i-1,c}[/math], а также к первой и последней формулам из (5.16), имеющим второй порядок аппроксимации.

---

Неявные алгоритмы численного дифференцирования

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

Первые производные по заданной сеточной функции [math]y_{i}= f(x_{i}),~ i=overline{0,n}[/math], можно вычислить несколькими способами:

а) следствием параметрических соотношений параболических сплайнов является система линейных алгебраических уравнений относительно первых производных, где [math]Delta f_{i}= f_{i+1}-f_{i};~ h_{i+1}= x_{i+1}-x_{i}colon[/math]

При заданных значениях производных на концах отрезка [math][x_{0},x_{n}][/math] эта система позволяет со вторым порядком аппроксимации вычислить значения первых производных [math]widehat{f},’_{i}[/math] во всех внутренних точках.

На регулярном шаблоне при [math]h=text{const}[/math] эта система упрощается:

[math]widehat{f},’_{i-1}+ 2widehat{f},’_{i}+ widehat{f},’_{i+1}= frac{2}{h}(f_{i+1}-f_{i-1}),quad i=overline{1,n-1};[/math]

б) на регулярном шаблоне значения первых производных [math]widehat{f},’_{i}[/math] со вторым порядком аппроксимации могут быть вычислены также из системы, удовлетворяющей свойству преобладания диагональных элементов:

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

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

[math]frac{widehat{f},’_{i-1}}{h_{i}}+ 2! left(frac{1}{h_{i+1}}+ frac{1}{h_{i}}right)! widehat{f},’_{i}+ frac{widehat{f},’_{i+1}}{h_{i+1}}= 3! left(frac{Delta f_{i}}{h_{i+1}^2}+ frac{Delta f_{i-1}}{h_{i}^2}right)!,quad i=overline{1,n-1}.[/math]

Данная система, удовлетворяющая свойству преобладания диагональных элементов, может быть замкнута либо значениями [math]widehat{f},’_{0},,widehat{f},’_{n}[/math], либо двумя функционально-дифференциальными граничными соотношениями (4.77);

г) первые производные со вторым порядком аппроксимации могут вычисляться также по значениям интегралов с использованием системы, подобной системе (4.72), в которой порядок производных понижен на единицу (в этом случае вместо [math]Delta f_{i}[/math] и [math]Delta f_{i-1}[/math] следует подставить [math]Delta F_{i}= I_{i}^{i+1}[/math] и [math]Delta F_{i-1}= I_{i-1}^{i}[/math]):

Замыкающие соотношения формируются аналогично предыдущему случаю. На регулярном шаблоне система (5.40) для [math]widehat{f},’_{i}[/math] записывается через интегральные приращения:

[math]widehat{f},’_{i-1}+ 4widehat{f},’_{i}+ widehat{f},’_{i+1}= frac{6}{h}(Delta I_{i}^{i+1}),quad i=overline{1,n-1}.[/math]

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

Данная система, удовлетворяющая условию преобладания диагональных элементов, может быть замкнута либо известными значениями [math]widehat{f},”_{0},, widehat{f},”_{n}[/math], либо двумя граничными соотношениями, следующими из равенства

[math]frac{Delta widehat{f},”_{i}}{h_{i+1}}= frac{Delta widehat{f},”_{i-1}}{h_{i}}quad left(frac{widehat{f},”_{i+1}-widehat{f},”_{i}}{h_{i+1}}= frac{widehat{f},”_{i}-widehat{f},”_{i-1}}{h_{i}}right)!,[/math]

которое записывается для трех крайних точек сетки [math]Omega_{n}[/math], т. е. для [math](x_0,x_1,x_2)[/math] и [math](x_{n-2},x_{n-1},x_{n})[/math].

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

---

Методика вычисления производных по неявным алгоритмам

1. С учетом характера задания сеточной функции (заданы значения функции [math]y_{i}= f(x_{i}),~ i= overline{0,n}[/math], или значения интегралов [math]I_{i}^{i+1},~ i= overline{0,n-1}[/math]) и порядка аппроксимации [math]t[/math], который необходимо обеспечить в алгоритме, выбрать систему алгебраических уравнений относительно значений производных во всех внутренних узлах сетки. При этом, если эта система является следствием параболических сплайнов, то для первых производных [math]t=2[/math], для вторых производных [math]t=1[/math]. В случае, когда система получена из кубических сплайнов, порядок [math]t[/math] возрастает на единицу. Данные системы являются незамкнутыми (число неизвестных превышает на два число уравнений).

2. Замкнуть выбранную систему двумя граничными условиями на концах сетки [math]Omega_{n}[/math]. Эти условия могут выбираться либо в виде явных аппроксимационных формул, либо в виде двух дополнительных алгебраических соотношений, включающих по два слагаемых с [math]f_{0}^{(p)},, f_{1}^{(p)}[/math] и [math]f_{n-1}^{(p)},, f_{n}^{(p)},, p=1,,2[/math]. При этом необходимо соблюсти соответствие порядков аппроксимации последних соотношений (или формул) и исходного порядка [math]t[/math].

3. Решить полученную замкнутую систему линейных алгебраических уравнений трехдиагонального вида методом прогонки.

Приведенная методика является составной частью построения кубического дифференциального сплайна. Она была применена при решении примера 4.13.

Математический форум (помощь с решением задач, обсуждение вопросов по математике).

Если заметили ошибку, опечатку или есть предложения, напишите в комментариях.