Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3791 / 1086 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Специальности: Программист, Математик
Лекция 2:

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

1.6. Задача численного дифференцирования

В пункте "1.2" уже была введена простейшая формула численного дифференцирования. Рассмотрим задачу приближенного вычисления значения производной подробнее.

Пусть задана таблица значений xi. В дальнейшем совокупность точек на отрезке, на котором проводятся вычисления, иногда будут называться сеткой, каждое значение xiузлом сетки. Пусть сетка — равномерная, и расстояние между узлами равно hшагу сетки. Пусть узлы сетки пронумерованы в порядке возрастания, т.е.

x0 = a, 
xj = a + jh,, j = 0, 1, ... , N

Пусть f(xj) = fj — функция, определенная в узлах сетки. Такие функции будут называться табличными, или сеточными функциями. Считаем, кроме того, что рассматриваемая сеточная функция есть проекция (или ограничение) на сетку некоторой гладкой нужное число раз непрерывно дифференцируемой функции f(x). По определению производной

$ f^{\prime}(x) = \lim\limits_{\Delta x \to 0}\frac{f(x + \Delta x) - f(x)}{\Delta x}, $

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

$ f^{\prime}(x_j)  \approx  \frac{f(x_j  + h) - f(x_j)}{h}. $ ( 1.1)

Если параметр h достаточно мал, то можно считать полученное значение производной достаточно точным. Погрешность формулы (1.1) оценена в пункте "1.2" . Как показано выше, при уменьшении шага сетки h ошибка будет уменьшаться, но при некотором значении h ошибка может возрасти до бесконечности. При оценке погрешности метода обычно считается, что все вычисления были точными. Но существует ошибка округления. При оценке ее большую роль играет машинный \varepsilon — мера относительной погрешности машинного округления, возникающей из-за конечной разрядности мантиссы при работе с числами в формате с плавающей точкой. Напомним, что по определению машинным \varepsilon называют наибольшее из чисел, для которых в рамках используемой системы вычислений выполнено 1 + \varepsilon  = 1. Тогда абсолютная погрешность при вычислении значения функции (или представлении табличной функции) есть f(x_j) \cdot \varepsilon . Максимальный вклад погрешностей округления при вычислении производной по формуле (1.1) будет

\frac{2f(x_j) \cdot \varepsilon }{h},
тогда, когда члены в знаменателе (1.1) имеют ошибки разных знаков.

Пусть k = max |f(x)|, максимум ищется на отрезке, на котором вычисляются значения производных. Тогда суммарная ошибка, состоящая из погрешности метода и погрешности округления, есть

$ \Delta = \frac{2k\varepsilon }{h} + \frac{M_2 h}{2},\quad M_2  = \max \left|{f^{\prime\prime}(x)}\right|.$

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

$\frac{M_2}{2} - \frac{2k\varepsilon }{h^2} = 0 $
, откуда

$h_{opt} = 2\sqrt{\frac{k\varepsilon }{M_2}} . $

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

$ f^{\prime}(x) = \lim\limits_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x - \Delta x)}{2\Delta x}.$

По аналогии напишем конечно-разностную формулу

$ f^{\prime}(x_j)  \approx  \frac{f(x_j + h) - f(x_j - h)}{2h}. $ ( 1.2)

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

$ f(x_j + h) = f(x_j) + f^{\prime}(x_j)h + f^{\prime\prime}(x_j) \frac{h^2}{2} + f^{\prime\prime\prime}(x_j + \theta_1 h) \frac{h^3}{6}, \\ 
f(x_j - h) = f(x_j) - f^{\prime}(x_j)h + f^{\prime\prime}(x_j) \frac{h^2}{2} - f^{\prime\prime\prime}(x_j  - \theta_2 h) \frac{h^3}{6}.$

Погрешность метода определяется 3-й производной функции. Введем M_3  = \max\limits_{x \in [a,b]} \left|{f^{\prime\prime\prime}(x)}\right|, тогда суммарная погрешность при вычислении по формуле с центральной разностью есть

$  \Delta = \frac{M_3 h^2 }{3} + \frac{k\varepsilon}{h},$

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

$ \frac{2M_3 h}{6} - \frac{k\varepsilon }{h^2} = 0 $
, откуда
$ h_{opt}  = \sqrt[3]{\frac{3k\varepsilon }{M_3}}$.

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

Формула (1.1) — двухточечная, (1.2) — трехточечная: при вычислении производной используются точки (узлы) xj (узел входит с нулевым коэффициентом), xj + h, xj - h — совокупность узлов, участвующих в каждом вычислении производной, в дальнейшем будем иногда называть сеточным шаблоном.

Введем на рассматриваемом отрезке шаблон из нескольких точек.

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

$ f^{\prime}(x_j)  \approx  \frac{1}{h}\sum\limits_{k = - l}^{m} \alpha _k f(x_j + kh), $ ( 1.3)

шаблон включает l точек слева от рассматриваемой точки xj и m справа. Коэффициенты \alpha_kнеопределенные коэффициенты. Формула дифференцирования может быть и односторонней — либо l, либо m могут равняться нулю. В первом случае иногда называют (на наш взгляд, не слишком удачно) такую приближенную формулу формулой дифференцирования вперед, во втором — формулой дифференцирования назад. Потребуем, чтобы (1.3) приближала первую производную с точностью O(hl + m) . Используем разложения в ряд Тейлора в окрестности точки xj. Подставляя их в (1.3), получим

$ \frac{1}{h}\sum\limits_{k =  - l}^{m} \alpha_k f(x_j + kh)  = \frac{1}{h}f(x_j)\sum \alpha_k + f^{\prime}(x_j)\sum k\alpha_k + f^{\prime\prime}(x_j)\sum \frac{k^2}{2}\alpha_k h +  \\ 
+ f^{\prime\prime\prime}(x_j)\sum \frac{k^3}{6}\alpha_k h^2  +  \ldots  + f^{\left( n\right)}(x_j) \sum \alpha_k \frac{k^n}{n!}h^{n - 1} +  \ldots $

Потребуем выполнение условий:

$ \sum \alpha_k = 0, \sum k\alpha_k = 1, \sum \alpha_k \frac{k^2}{2} = 0, \ldots \sum \alpha_k \frac{k^n}{n!} = 0, \ldots $ ( 1.4)

Получаем систему линейных алгебраических уравнений для неопределенных коэффициентов \alpha (1.4). Матрица этой системы есть

\left( \begin{array}{cccc}
   1 & 1 &  \ldots  & 1  \\
   - l &  - l + 1 &  \ldots  & m  \\
   l^2 & {(l - 1)}^2 &  \ldots  & m^2   \\
   {(- l)}^3  &  {(- l + 1)}^3 & \ldots  &    m^3  \\
    \ldots  &  \ldots  &  \ldots  &  \ldots   \\
 \end{array} \right) .

Вектор правых частей (0, 1, 0, ..., 0)T.

Определитель данной матрицы — детерминант Вандермонда. Из курса линейной алгебры следует, что он не равен нулю. Тогда существует единственный набор коэффициентов \alpha, который позволяет найти на шаблоне из (1 + l + m) точек значение первой производной с точностью O(hl + m).

Для нахождения второй производной можно использовать ту же самую формулу (1.3) с небольшой модификацией

$ f^{\prime}(x_j)  \approx  \frac{1}{h^2}\sum\limits_{k = - l}^{m} \alpha_k f(x_j + kh) , $

только теперь

$ \sum k\alpha_k = 0, \sum \alpha_k k^2  = \frac{1}{2} $.

Очевидно, что и данная система уравнений для нахождения неопределенных коэффициентов имеет единственное решение. Для получения с той же точностью приближенных значений производных до порядка l + m включительно с точностью O(hl + m) модификации формулы (1.3) и условий (1.4) очевидны, набор неопределенных коэффициентов находится единственным образом.

Таким образом, доказано следующее утверждение. На сеточном шаблоне, включающем в себя N + 1 точку, с помощью метода неопределенных коэффициентов всегда можно построить единственную формулу для вычисления производной от первого до n порядка включительно с точностью O(hN).

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

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

Читателям предлагается оценить значение оптимального шага при вычислениях по формулам типа (1.3) самостоятельно.