Численное дифференцирование

Постановка задачи для 1-ой производной

Пусть таблично задана некоторая функция \(u(x)\):

\[\begin{split} \begin{array}{|l|l|l|l|l|l|l|l|} \hline u(x) & u\left(x_0\right) & u_1 & u_2 & \cdots & u_j & \cdots & u_J \\ \hline x & x_0 & x_1 & x_2 & \cdots & x_j & \cdots & x_J \\ \hline \end{array} \end{split}\]

Как раз в такой форме задана любая функция в компьютере (массив «иксов» и массив «игреков»). Отметим, что разбиение по «иксам» может быть как равномерным, так и неравномерным - всё зависит от конкретной табдично заданной функции.

Хотим вычислить её производную в каждой точке \(x_j\). Вспомним определение из матанализа:

\[u'\left(x_j\right)=\lim _{\Delta x \rightarrow 0} \frac{u\left(x_j+\Delta x\right)-u\left(x_j\right)}{\Delta x}\]

Но ведь на компьютере \(\Delta x\) не может стремится к нулю, т.к. функция по сути дискретная, - придётся считать приближённо, с помощью аппроксимации производной путём использования конечно-разностных схем. Некоторые примеры:

a) \(u_j^{(1)}=\cfrac{u_j-u_{j-1}}{x_j-x_{j-1}}=\cfrac{u_j-u_{j-1}}{h}\)

b) \(u_j^{(1)}=\cfrac{u_{j+1}-u_j}{x_{j+1}-x_j}=\cfrac{u_{j+1}-u_j}{h}\)

c) \(u_j^{(1)}=\cfrac{u_{j+1} - u_{j-1}}{x_{j+1}-x_{j-1}}=\cfrac{u_{j+1}-u_{j-1}}{2 h}\)

Подчеркнём разницу между обозначениями \(u'\) и \(u^{(1)}\). \(u'\) - это истинная (математическая) производная, а \(u^{(1)}\) - это лишь аппроксимация (наше приблежение).

Порядок аппроксимации

Чем же отличаются данные схемы? Оказывается, точностью. Но сначала вспомним про ряд Тейлора с остаточным членом в форме Лагранжа и Пеано:

\[ \begin{aligned} &u\left(x_j+\Delta x\right)=u\left(x_j\right)+u^{\prime}\left(x_j\right) \cdot \Delta x+u^{\prime \prime}\left(x_j\right) \frac{\Delta x^2}{2 !} +u^{\prime \prime \prime}\left(x_j\right) \frac{\Delta x^3}{3 !}+\cdots+u^{(n)}\left(x_i\right) \frac{\Delta x^n}{n !}+\frac{u^{(n+1)}(\xi)}{(n+1) !} \Delta x^{n+1}=\sum_{k=0}^n \frac{u^{(k)}\left(x_j\right)}{k !} \Delta x^k+O\left(\Delta x^{n+1}\right) \end{aligned} \]

Подставляя это разложение в наши примеры, посмотрим на их «точность»:

a) Левая разность

a) \(u_j^{(1)}=\cfrac{u_j-u_{j-1}}{x_j-x_{j-1}}=\cfrac{u_j-u_{j-1}}{h}\)

\( u_{j-1}=u\left(x_{j-1}\right)=u(x_j - h)=u_j-h u_j^{\prime}+\frac{h^2}{2} u_j^{\prime \prime}-\frac{h^3}{6} u_j^{\prime \prime \prime}+O\left(h^4\right) \)

\(u_j^{(1)}=\cfrac{u_j-u_{j-1}}{h}=\cfrac{u_j-(u_j-h u_j^{\prime}+\frac{h^2}{2} u_j^{\prime \prime}-\frac{h^3}{6} u_j^{\prime \prime \prime}+O\left(h^4\right))}{h} = u_j^{\prime}-\frac{h u_j^{\prime \prime}}{2}+\frac{h^2}{6} u_j^{\prime \prime \prime} +o(h^3) \approx u_j^{\prime}+O(h)\)

Видим, что при \(h \rightarrow 0\) наша аппроксимация стремится к истинной производной. В данном случае эта разностная схема имеет 1-ый порядок аппроксимации, т.к. ошибка пропорциональна первой степени \(h\).

b) Правая разность

b) \(u_j^{(1)}=\cfrac{u_{j+1}-u_j}{x_{j+1}-x_j}=\cfrac{u_{j+1}-u_j}{h}\)

\( u_{j+1}=u\left(x_{j+1}\right)=u(x_j + h)=u_j+h u_j^{\prime}+\frac{h^2}{2} u_j^{\prime \prime}+\frac{h^3}{6} u_j^{\prime \prime \prime}+O\left(h^4\right) \)

\(u_j^{(1)}=\cfrac{u_{j+1}-u_j}{h}=\cfrac{(u_j+h u_j^{\prime}+\frac{h^2}{2} u_j^{\prime \prime}+\frac{h^3}{6} u_j^{\prime \prime \prime}+O\left(h^4\right)) - u_j}{h} = u_j^{\prime}+\frac{h u_j^{\prime \prime}}{2}+\frac{h^2}{6} u_j^{\prime \prime \prime} +O(h^3) \approx u_j^{\prime}+O(h)\)

Опять получился 1-ый порядок аппроксимации.

c) Центральная разность

c) \(u_j^{(1)}=\cfrac{u_{j+1} - u_{j-1}}{x_{j+1}-x_{j-1}}=\cfrac{u_{j+1}-u_{j-1}}{2 h}\)

\( u_{j+1}=u\left(x_{j+1}\right)=u(x_j + h)=u_j+h u_j^{\prime}+\frac{h^2}{2} u_j^{\prime \prime}+\frac{h^3}{6} u_j^{\prime \prime \prime}+O\left(h^4\right) \)

\( u_{j-1}=u\left(x_{j-1}\right)=u(x_j - h)=u_j-h u_j^{\prime}+\frac{h^2}{2} u_j^{\prime \prime}-\frac{h^3}{6} u_j^{\prime \prime \prime}+O\left(h^4\right) \)

\(u_j^{(1)}=\frac{u_{j+1}-u_{j-1}}{2h}=u_j^{\prime}+\frac{h^2}{6} u_j^{\prime \prime \prime}+o(h^2)=u_j^{\prime}+O(h^2)\)

Получился 2-ой порядок аппроксимации - то есть эта разностная аппроксимация на порядок точнее двух предыдущих. К примеру, без учёта констант перед \(h\), если в a) и b) ошибка была бы \(10^{-3}\), то в с) на том же наборе данных мы бы имели ошибку \(10^{-6}\).

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

Постановка задачи для k-ой производной

Пусть таблично задана некоторая функция \(u(x)\):

\[\begin{split} \begin{array}{|l|l|l|l|l|l|l|l|} \hline u(x) & u\left(x_0\right) & u_1 & u_2 & \cdots & u_* & \cdots & u_J \\ \hline x & x_0 & x_1 & x_2 & \cdots & x_* & \cdots & x_J \\ \hline \end{array} \end{split}\]

Хотим приближенно найти производную \(k\)-го порядка в точке \(x_*\). Будем искать решение в виде конечно-разностной схемы:

\[u^{(k)}\left(x_*\right)=u_*^{(k)}=\frac{\overset{m}{\underset{j \,=-l}{\sum}} \alpha_j \cdot u\left(x_*+\Delta x_j\right)}{h^k}\]

Здесь \(\alpha_j\) - некоторые числовые коэффициенты, которые определяют вид конечно-разностной схемы. Как же их найти?

Алгоритм составления разностной схемы

1. Фиксировать \(l\) и \(m\) исходя из п.3 ниже.

2. Разложить \(u\left(x_*+\Delta x_j\right)\) в ряд Тейлора в точке \(x_*\). Привести множители перед производными.

3. Составить систему линейных уравнений на \(\alpha_j\) исходя из

   • зануления скобок перед производными порядка меньше k

   • равенства скобки единице перед искомой k-ой производной

   • занулению скобок перед некоторым количеством производных порядка выше k для увеличения точности.

4. Решить эту систему аналитически \(\rightarrow\) найти \(\alpha_j\) \(\rightarrow\) схема задана.

Если же в нашей функции нет фиксированного шага по «искам» \(h\), то меняется лишь изначальная формула, откуда \(h\) пропадает из знаменателя, и остаётся только числитель. В остальном алгоритм такой-же. Более подробно см. лекцию С1.

Вообще, все эти задачи для различных \(k, l, m\) давно решены за нас. Достаточно загуглить «конечно-разностные схемы». На википедии есть большой список схем на все случаи жизни.

k-ая производная функции нескольких переменных

В целом «дух», алгоритм и понятие порядка аппроксимации абсолютно такие же как и в предыдущих пунктах, только теперь необходимо рассматривать различные шаги \(h_1, h_2, ...\) по разным осям. Также в рядах Тейлора появятся перекрестные слагаемые. За более подробным описанием см. лекцию С1.