7. Чисельне диференціювання

7.1. Побудова формул чисельного диференціювання

Задача чисельного диференціювання виникає у випадку коли необхідно обчислити похідну функції, значення якої задані таблицею. Нехай задано

\begin{equation} f_i = f(x), \quad i = \overline{0,n}, \quad x_i \in [a, b]. \end{equation}

Проінтерполюємо ці значення. Тоді

\begin{equation} \label{eq:7.1} f(x) = L_n(x) + r_n(x), \end{equation}

де залишковий член у формі Ньютона має вигляд:

\begin{equation} r_n(x) = f(x; x_0, \ldots, x_n) \omega_n(x), \end{equation}

де \begin{equation} \omega_n(x) = \Prod_{i = 0}^n (x - x_i). \end{equation}

Звідси

\begin{equation} \label{eq:7.2} f^{(k)}(x) = L_n^{(k)}(x) + r_n^{(k)}(x). \end{equation}

За наближене значення похідної в точці $x$ беремо $f^{(k)}(x) \approx L_n^{(k)}(x)$, $x \in [a, b]$.

Оцінимо похибку наближення — $r^{(k)}(x)$. За формулою Лейбніца:

\begin{equation} \label{eq:7.3} r_n^{(k)}(x) = \Sum_{j = 0}^k C_k^j f^{(j)}(x; x_0, \ldots, x_n) \omega_n^{(k - j)}(x). \end{equation}

З властивості розділених різниць маємо для $f(x) \in C^{n + k + 1}[a,b]$:

\begin{equation} f^{(j)}(x; x_0, \ldots, x_n) = j! \cdot f ( \underset{j + 1}{\underbrace{x, \ldots, x}}; x_0, \ldots, x_n ) = \frac{j!}{(n + j + 1)!} \cdot f^{(n + j + 1)}(\xi_j), \end{equation}

де всі $x_j \in [a, b]$.

Остаточно вираз для похибки наближення похідної має вигляд:

\begin{equation} \label{eq:7.4} r_n^{(k)}(x) = \Sum_{j = 0}^k \frac{k!}{(k - j)! (n + j + 1)!} \cdot f^{(n + j + 1)})(\xi_j) \omega_n^{(k - j)}(x). \end{equation}

Оцінка похибки матиме вигляд:

$$\begin{equation} \left| f^{(k)}(x) - L_n^{(k)}(x) \right| \le M \Sum_{j = 0}^k \frac{k!}{(k - j)! (n + j + 1)!} \cdot \left| \omega_n^{(k - j)}(x) \right|, \end{equation}$$

де $M = \Max_{0 \le j \le k} \Max_{x \in [a, b]} \left\vert f^{(n + j + 1)}(x) \right\vert$.

Нагадаємо, що процес інтерполювання розбіжний. Крім того, якщо $k > n$, то $L_n^{(k)}(x) \equiv 0$. Тому не можна брати великими значення $n$ та $k$. Як правило $k = 1, 2$, іноді $k = 3, 4$. Відповідно, $n = k$, або $n = k + 1$, або $n = k + 2$.

Подивимося як залежить порядок збіжності процесу чисельного диференціювання від кроку. Нехай $x_i = x_0 + i h$, $h > 0$ — крок. Тоді за умови $x_n - x_0 = O(h)$:

\begin{equation} \omega_n(x) = (x - x_0) \cdot \ldots \cdot (x - x_n) = O(h^{n + 1}), \end{equation}

де $x \in [x_0, x_n]$.

Перша похідна від $\omega_n(x)$ має порядок на одиницю менше, тобто

\begin{equation} \omega_n’(x) = O(h^n). \end{equation}

Далі

\begin{equation} r_n^{(k)}(x) = O(h^{n + 1 - k}), \end{equation}

тому

\begin{equation} f^{(k)}(x) - L_n^{(k)} = O(h^{n + 1 - k}). \end{equation}

При умові $n \ge k$ останній вираз збігається до нуля, тобто

\begin{equation} \label{eq:7.5} f^{(k)}(x) - L_n^{(k)}(x) \xrightarrow[h \to 0]{} 0. \end{equation}

Далі

\begin{equation} r_n^{(k)}(x) = \underset{O(h^{n + 1 - k})}{\underbrace{f(x; x_0, \ldots, x_n) \omega_n^{(k)(x)}}} + \underset{O(h^{n + 2 - k})}{\underbrace{\Sum_{j = 1}^k C_k^j f^{(j)}(x; x_0, \ldots, x_n) \omega_n^{(k - j)}(x)}}. \end{equation}

Якщо

\begin{equation} \omega_n^{(k)}\left(\overline{x}\right) = 0, \end{equation}

то

\begin{equation} \label{eq:7.6} r_n^{(k)}\left(\overline{x}\right) = O(h^{n + 2 - k}). \end{equation}

Точки $x = \overline{x}$ називаються точками підвищеної точності формул чисельного диференціювання.

Приклад 1: Виведемо формули чисельного диференціювання для $k = 1$, $n = 1$.

Виберемо точки $x_0$, $x_1 = x_0 + h і$ інтерполяційний багаточлен має вигляд:

\begin{equation} L_1(x) = f_0 + (x - x_0) \cdot \frac{f_1 - f_0}{h}. \end{equation}

Для похідної отримаємо вираз:

\begin{equation} f’(x) \approx L_1’(x) = \frac{f_1 - f_0}{h}, \quad x \in [x_0, x_1]. \end{equation}

Розписавши за формулою Тейлора, отримаємо вираз для похибки:

\begin{equation} r_1’(x) = \frac{f^{(3)}(\xi_1)}{3!} \cdot (x - x_0) (x - x_1) + \frac{f^{(2)}(\xi_0)}{2!} \cdot (2 x - x_1 - x_0) = O(h). \end{equation}

Якщо $2 \overline{x} - x_1 - x_0 = 0$, то $r_1'\left(\overline{x}\right) = O(h^2)$. Тобто $\overline{x} = \frac{x_1 + x_0}{2}$ — точка підвищеної точності. Більш точно (див. приклад 3):

\begin{equation} \left| r_1’\left(\overline{x}\right) \right| \le \frac{h^2 M_3}{24}, \end{equation}

де $M_3 = \Max_{x \in [a, b]} \left\vert f^{(3)}(x) \right\vert$.

Приклад 2: Аналогічно виведемо формули чисельного диференціювання для $k = 1$, $n = 2$.

Виберемо точки $x_0$, $x_1 = x_0 + h$, $x_2 = x_0 + 2 h$. Інтерполяційний поліном має вигляд:

\begin{equation} L_2(x) = f_0 + (x - x_0) \cdot \frac{f_1 - f_0}{h} + (x - x_0) (x - x_1) \cdot \frac{f_2 - 2 f_1 + f_0}{2 h^2}. \end{equation}

Тоді замінимо $f'(x) \approx L_2;(x) = \frac{f_1 - f_0}{h} + (2 x - x_0 - x_1) \cdot \frac{f_2 - 2 f_1 + f_0}{2 h^2}$, $x \in [x_0, x_2]$.

Якщо сюди підставити $x = x_0$, то отримаємо $f'(x_0) \approx \frac{-f_2 + 4 f_1 - 3 f_0}{2 h}$. Для точки $x = x_1$ маємо $f'(x_1) \approx \frac{f_2 - f_0}{2 h} = f_{x,1}^0$. Для точки $x = x_2$ маємо $f'(x_2) \approx \frac{f_0 - 4 f_1 + 3 f_2}{2 h}$, $x \in [x_0, x_2]$. Для похибки маємо оцінку $r_2'(x) = O(h^2)$.

Позначимо

• для $x \in [x_i, x_{i + 1}]$, $f_{x, i} = \frac{f_{i + 1} - f_i}{h} \approx f'(x)$, (різницева похідна вперед);

• для $x \in [x_{i - 1}, x_i]$ — $f_{\overline{x}, i} = \frac{f_i - f_{i - 1}}{h} \approx f'(x)$, (різницева похідна назад);

• для $x \in [x_{i - 1}, x_{i + 1}]$ — $f_{x,i}^0 \approx f'(x)$ (центральна різницева похідна).

Замість $f'(x_i)$ можна взяти будь-яке із значень: $f_{x,i}$, $f_{\overline{x}, i}$ або $f_{x,i}^0$.

Задача 21: Знайти точки підвищеної точності формул чисельного диференціювання для $k = 1$, $n = 2$ і оцінити похибку в цих точках.

Приклад 3: При $n = 1$, $k = 1$ оцінимо точність формул чисельного диференціювання за формулою Тейлора.

1. Нехай $f(x) \in C^2([a, b])$. Тоді

   \begin{equation} f’(x_0) - \frac{f(x_0 + h) - f(x_0)}{h} = f’(x_0) - \frac{1}{h} \left( f_0 + h f_0’ + \frac{h^2}{2} f’’(\xi) - f_0 \right) = - \frac{h}{2} \cdot f’’(\xi), \end{equation}

   \begin{equation} \left| f’(x_0) - \frac{f_1 - f_0}{h} \right| \le \frac{M_2 h}{2}, \end{equation}

   де $M_2 = \Max_{[x_0, x_1]} \vert f''(\xi) \vert$.

2. Нехай $f(x) \in C^3([a,b])$. Тоді, розписавши розклад по формулі Тейлора до третьої похідної, маємо оцінку:

   $$\begin{equation} \begin{aligned} & f'\left(\overline{x}\right) - \frac{f(x_0 + h) - f(x_0)}{h} = \newline & \quad = f'\left(\overline{x}\right) - \frac{1}{h} \left( f'\left(\overline{x}\right) + \frac{h}{2} \cdot f'\left(\overline{x}\right) + \frac{h^2}{8} \cdot f''\left(\overline{x}\right) + \frac{h^3}{48} \cdot f'''(\xi)) - \right. \newline & \qquad \left. - f\left(\overline{x}\right) + \frac{f}{2} \cdot f'\left(\overline{x}\right) - \frac{h^2}{8} \cdot f''\left(\overline{x}\right) + \frac{h^3}{48} \cdot f'''(\eta) \right) = -\frac{h^2}{24} \cdot f'''(\zeta). \end{aligned} \end{equation}$$

   \begin{equation} \left| f’\left(\overline{x}\right) - \frac{f_1 - f_0}{h} \right| \le \frac{h^2 M_3}{24}, \end{equation}

   де $\overline{x} = \frac{x_1 + x_0}{2}$.

Задача 22: Показати, що якщо $f(x) \in C^3([a,b])$, то $\left\vert f'(x_1) - \frac{f_2 - f_0}{2 h} \right\vert \le \frac{M_3 h^2}{6}$.

Приклад 4: При $n = 2$, $k = 2$ маємо:

\begin{equation} \begin{aligned} L_2(x) &= f_{i - 1} + \frac{f_i - f_{i - 1}}{h} \cdot (x - x_{i - 1}) + \newline &\quad + \frac{f_{i + 1} - 2 f_i + f_{i - 1}}{h^2} \cdot (x - x_{i - 1}) \cdot (x - x_i) \end{aligned} \end{equation}

$$\begin{equation} L_2''(x) = \frac{f_{i + 1} - 2 f_i + f_{i - 1}}{h^2}. \end{equation}$$

Для $f (x) \in C^4([a,b])$ оцінимо точність формул чисельного диференціювання за формулою Тейлора:

$$\begin{equation} \begin{aligned} & f''(x_1) - \frac{f_2 - 2 f_1 + f_0}{h^2} = \newline & \quad = f''(x_1) - \frac{f(x_1 + h) - 2 f(x_1) + f(x_1 - h)}{h^2} = \newline & \quad = f''(x_1) - \frac{1}{h^2} \left( f_1 + h f_1'' + \frac{h^2}{2} \cdot f_1'' + \frac{h^3}{6} \cdot f_1''' + \frac{h^4}{24} \cdot f^{(4)}(\xi_1) - 2 f_1 + \right. \newline & \qquad \left. + f_1 - h f_1'' + \frac{h^2}{2} \cdot f_1'' - \frac{h^3}{6} \cdot f_1''' + \frac{h^4}{24 \cdot f^{(4)}}(\xi_2) \right) = \frac{h^2}{12} \cdot f^{(4)}(\xi), \end{aligned} \end{equation}$$

де $\xi_1, \xi_2, \xi \in [x_0, x_2]$.

Отже,

\begin{equation} \left\vert f_1’’ - \frac{f_2 - 2 f_1 + f_0}{h^2} \right\vert \le \frac{M_4 h^2}{12}. \end{equation}

Задача 23: Побудувати формулу чисельного диференціювання $k = 2$, $n = 2$ у випадку нерівновіддалених вузлів: $x_0$, $x_1 = x_0 + h_1$, $x_2 = x_1 + h_2$. Оцінити точність формули. Знайти точки підвищеної точності оцінити похибку.

Крім інтерполяційних формул для чисельного диференціювання можна застосовувати сплайни. Нехай $f_i = f(x_i)$. Побудуємо інтерполяційний сплайн першого степеня $s_1(x)$, для якого має місце оцінка $\left\vert f^{(k)}(x) - s_1^{(k)}(x) \right\vert = O(h^{2 - k})$, $k = 0, 1$. Звідси при $k = 1$ маємо $f'(x) - s_1'(x) = O(h)$.

Для кубічного інтерполяційного сплайну $s_3(x)$ маємо для першої та другої похідних:

\begin{equation} \left| f^{(k)}(x) - s_3^{(k)}(x) \right| - O(h^{4 - k}), \quad k = 1, 2. \end{equation}

7.2. Про обчислювальну похибку чисельного диференціювання

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

Перед цим зробимо зауваження про вплив збурення функції на значення звичайних похідних.

Нехай $f(x) \in C^1([a,b])$ і її збурення має вигляд:

\begin{equation} \tilde f(x) = f(x) + \frac{\sin (\omega x)}{n}. \end{equation}

При $n \to \infty$ маємо $\left\| f(x) - \tilde f(x)\right\|\void_{C([a,b])} = \frac{1}{n} \to 0$, звідси $\tilde f(x) \xrightarrow[n \to \infty]{} f(x)$. Таким чином це малі збурення. Маємо $\tilde f'(x) = f'(x) + \frac{\omega}{n} \cdot \cos(\omega x)$. Нехай $\omega = n^2$, тоді

\begin{equation} \left| f(x) - \tilde f(x)\right|\void_{C([a,b])} = \frac{|\omega|}{n} = n \xrightarrow[n \to \infty]{} \infty. \end{equation}

Цей приклад ілюструє нестійкість оператора диференціювання. Є сподівання, що ця нестійкість має місце і для чисельного диференціювання.

Нехай $\tilde f_i = f_i + \delta_i$, $f_i = f(x_i)$, $i = \overline{0, n}$, $\vert \delta_i \vert \le \delta$. Розглянемо вплив похибок $\delta_i$ на конкретних формулах чисельного диференціювання.

Приклад 1: Оцінімо вплив збурень на похибку обчислення першої похідної $n = 1$, $k = 1$.

$$\begin{equation} f_i' = \frac{\tilde f_i - \tilde f_{i - 1}}{h} = f_i' - \frac{f_i - f_{i - 1}}{h} - \frac{\delta_i - \delta_{i - 1}}{h}, \end{equation}$$ $$\begin{equation} \label{eq:7.2.1} \left| f_i' - \frac{\tilde f_i - \tilde f_{i - 1}}{h} \right| \le \left| f_i' - \frac{f_i - f_{i - 1}}{h} \right| + \left| \frac{\delta_i - \delta_{i - 1}}{h} \right| \le \frac{M_2 h}{2} + \frac{2 \delta}{h} \xrightarrow[h \to 0]{} \infty \end{equation}$$

Таким чином, як і для аналітичного диференціювання, маємо некоректність: при малих збуреннях $\vert \delta_i \vert le \delta_i$ можуть бути як завгодно великі похибки, якщо $\frac{\delta}{h} \to \infty$ при $h \to 0$.

Мінімізуємо вплив цих збурень. Позначимо

\begin{equation} \varphi(h) = \frac{M_2 h}{2} + \frac{2 \delta}{h}. \end{equation}

Тоді мінімум цієї функції досягається для таких $h$:

\begin{equation} \varphi’(h) = \frac{M_2}{2} - \frac{2 \delta}{h^2} = 0, \end{equation}

звідки $h_0 = 2 \sqrt{\frac{\delta}{M_2}}$. При такому значенні $h$ оцінка похибки $(35)$ така:

\begin{equation} \varphi(h_0) = 2 \sqrt{M_2 \delta} = O\left(\sqrt{\delta}\right) \xrightarrow[\delta \to 0]{} 0. \end{equation}

Приклад 2: Подивимося на вплив збурень на похибку обчислення першої похідної при використанні центральної різницевої похідної.

$$\begin{equation} \left| f_i' - \frac{\tilde f_{i + 1} - \tilde f_{i - 1}}{2 h} \right| \le \left| f_i' - \frac{f_{i + 1} - f_{i - 1}}{2 h} \right| + \left|\frac{\delta_{i + 1} - \delta_{i - 1}}{2 h} \right| \le \frac{M_3 h^2}{6} + \frac{\delta}{h} = \varphi(h). \end{equation}$$

З рівняння $\varphi'(h) = \frac{M_3 h}{3} - \frac{\delta}{h^2} = 0$ маємо: $h_0^3 = \frac{3\delta}{M_3}$, $h_0 = \sqrt[3]{\frac{3 \delta}{M_3}}$. Отже,

$$\begin{equation} \varphi(h_0) = \frac{M_3}{6} \sqrt[3]{\frac{9 \delta^2}{M_3^2}} + \frac{\delta}{\sqrt[3]{\frac{3\delta}{M_3}}} = \frac{1}{2} \sqrt[3]{\frac{M_3 \delta^2}{3}} + \sqrt[3]{\frac{M_3 \delta^2}{3}} = \frac{3}{2} \sqrt[3]{\frac{M_3 \delta^2}{3}} = O\left(\sqrt[3]{\delta^2}\right). \end{equation}$$

Таким чином швидкість збіжності при $\delta \to 0$ похибки формули чисельного диференціювання центральною похідною вища ніж для формули з прикладу 1 (похідна вперед або назад).

Задача 24: Дослідити похибку чисельного диференціювання для $n = 2$, $k = 2$, вибрати оптимальний крок $h_0$, дати оцінку $\varphi(h_0)$.