Содержание

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

Математическая постановка задачи

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

Математически производная функции ff определяется как

f(x)=limΔx0f(x+h)f(x)Δx.f'(x) = \underset{\Delta x \rightarrow 0}{\lim} \frac{f(x+h) - f(x)}{\Delta x}.

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

Для производной используют разные обозначения.

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

f(x,y)x=fx(x,y)=limΔx0f(x+Δx,y)f(x,y)Δx,f(x,y)y=fy(x,y)=limΔy0f(x,y+Δy)f(x,y)Δy.\frac{\partial f(x,y)}{\partial x} = f_x(x,y) = \underset{\Delta x \rightarrow 0}{\lim} \frac{f(x+\Delta x,y) - f(x,y)}{\Delta x},\\ \frac{\partial f(x,y)}{\partial y} = f_y(x,y) = \underset{\Delta y \rightarrow 0}{\lim} \frac{f(x,y+\Delta y) - f(x,y)}{\Delta y}.

Наконец, существуют смешанные производные — это производные от нескольких аргументов сразу

f(x,y)yx=fyx(x,y)=limΔx0fy(x+Δx,y)fy(x,y)Δx.\frac{\partial f(x,y)}{\partial y \partial x} = f_{yx}(x,y) = \underset{\Delta x \rightarrow 0}{\lim} \frac{f_y(x+\Delta x,y) - f_y(x,y)}{\Delta x}.

Для функции, у которой существуют производные fx,fy,fxy,fyxf_x,f_y,f_{xy},f_{yx}, частные производные равны между собой: fxy=fyxf_{xy}=f_{yx}. В общем случае это не так.

В многомерном случае для обозначения производной используют оператор набла

f(x,y)=[xf(x,y)yf(x,y)].\nabla f(x,y) = \left[\begin{array}{l}\frac{\partial}{\partial x} f(x,y)\\\frac{\partial}{\partial y} f(x,y) \end{array}\right].

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

2f(x,y)=f(x,y)=22xf(x,y)+22yf(x,y).\nabla^2 f(x,y) = \nabla\cdot\nabla f(x,y) = \frac{\partial^2}{\partial^2 x}f(x,y)+\frac{\partial^2}{\partial^2 y}f(x,y).

Здесь точкой обозначается скалярное произведение векторов. Вторую производную можно записать и в векторном виде, однако такая запись редка, и отдельного символа для нее нет:

f(x,y)=[22xf(x,y)22yf(x,y)].\nabla\circ\nabla f(x,y) = \left[\begin{array}{l}\frac{\partial^2}{\partial^2 x}f(x,y)\\\frac{\partial^2}{\partial^2 y}f(x,y)\end{array}\right].

Здесь кругом обозначается покомпонентное произведение векторов (произведение Адамара).

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

Конечные разности

Очевидно, взяв формулу производной из определения и некоторый малый шаг Δx\Delta x, можно получить простейшую аппроксимацию

fxf(x+Δx)f(x)Δx.f'_x \approx \frac{f(x + \Delta x) - f(x)}{\Delta x}.

Эта формула называется разностью вперед. Также существует разность назад

fxf(x)f(xΔx)Δx,f'_x \approx \frac{f(x) - f(x - \Delta x)}{\Delta x},

которая отличается лишь направлением шага, и центральная разность

fxf(x+Δx)f(xΔx)2Δx.f'_x \approx \frac{f(x + \Delta x) - f(x - \Delta x)}{2\Delta x}.

Все эти разности при Δx0\Delta x \rightarrow 0 дают точное определение производной. Конечные разности интуитивны, однако их обобщение для аппроксимаций более высокой точности неочевидно, и может происходить несколькими способами. Сейчас мы рассмотрим эти способы.

Аппроксимация с помощью ряда Тейлора

Ряд Тейлора ставит в соответствие значению функции в определенной точке линейную комбинацию значений ее производных. Наша идея состоит в том, чтобы разложить функцию в ряд Тейлора в нескольких точках, отличающихся на Δx\Delta x, и из получивишхся уравнений получить выражение для производной.

f(x)=n=0f(n)(a)n!(xa)n.f(x) = \sum\limits_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^n.

Для удобства будем использовать альтернативную запись этого ряда для x+Δxx + \Delta x

f(x+Δx)=n=0f(n)(x)n!(Δx)n.f(x + \Delta x) = \sum\limits_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!} (\Delta x)^n.

Теперь разложим функцию в точках x+Δxx + \Delta x и xΔxx - \Delta x до второго слагаемого и получим

f(x+Δx)=f(x)+f(x)Δx+O((Δx)2),f(xΔx)=f(x)f(x)Δx+O((Δx)2).f(x + \Delta x) = f(x) + f'(x) \Delta x + O((\Delta x)^2),\\ f(x - \Delta x) = f(x) - f'(x) \Delta x + O((\Delta x)^2).

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

f(x)=f(x+Δx)f(xΔx)2Δx+O((Δx)2).f'(x) = \frac{f(x+\Delta x) - f(x-\Delta x)}{2 \Delta x} + O((\Delta x)^2).

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

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

Полиномиальная аппроксимация

Более автоматизированным способом аппроксимации производных является представление функции в виде полинома. Мы записываем значение функции в точках x,x+Δx,x, x + \Delta x, \ldots с помощью полинома, решаем получившуюся систему уравнений, чтобы узнать коэффициенты полинома, а затем получаем формулы для производных путем дифференцирования полинома.

Запишем систему уравнений для функции в двух точках

f(x)=a1+a2x,f(x+Δx)=a1+a2(x+Δx).f(x) = a_1 + a_2 x,\\ f(x + \Delta x) = a_1 + a_2 (x + \Delta x).

Решая эту систему, получим

a1=(x+Δx)f(x)xf(Δx+x)Δx,a2=f(x+Δx)f(x)Δx.a_1 = \frac{(x+\Delta x) f(x)-x f\left(\Delta x+x\right)}{\Delta x}, \quad a_2 = \frac{f\left(x+\Delta x\right)-f(x)}{\Delta x}.

Первая производная от полинома равна a2a_2, поэтому это и есть наш финальный ответ:

f(x)=a2=f(x+Δx)f(x)Δx.f'(x) = a_2 = \frac{f\left(x+\Delta x\right)-f(x)}{\Delta x}.

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

Solve[{a1 + a2 x == f[x],a1 + a2 (x+dx) == f[x+dx]}, {a1,a2}] // TeXForm

Этот код без изменений работает в программе Mathematica. Аналогичный код можно написать и для Matlab. При должном усердии можно написать генератор формул конечных разностей для полиномов произвольных порядков.

Аппроксимация с помощью преобразования Фурье

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

F{f(x)}(t)=2πitF{f(x)}(t)f(x)=F1{2πitF{f(x)}(t)}(x)\mathcal{F}\left\{f'(x)\right\}(t) = 2 \pi i t \,\mathcal{F}\left\{f(x)\right\}(t) \Rightarrow f'(x) = \mathcal{F}^{-1}\left\{2 \pi i t \,\mathcal{F}\left\{f(x)\right\}(t)\right\}(x)

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

Производная по направлению

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

f(x)=limx2x10f(x2)f(x1)x2x1.f'(\vec{x}) = \underset{\left|\vec{x}_2-\vec{x}_1\right| \rightarrow 0}{\lim} \frac{f(\vec{x}_2) - f(\vec{x}_1)}{\left|\vec{x}_2-\vec{x}_1\right|}.

🌊🌊🌊

Мы узнали, что производные можно аппроксимировать различными способами и с различной точностью. Теперь мы рассмотрим задачи, в которых применяются эти аппроксимации, а именно дифференциальные уравнения в частных производных (ДУЧП).

Задания

Уравнение второго порядка v21 балл

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

d2f(t)dt2+αf(t)=0.\frac{d^2 f(t)}{dt^2} + \alpha f(t) = 0.

Значение функции ff и ее первой производной при t=0t=0 произвольное.

Затопление отсеков v32 балла

Решите задачу затопления отсеков с помощью метода Ньютона. Для вычисления первой и второй производной воспользуйтесь центральными разностями. Величину шага по zz выберите самостоятельно.

Производная по направлению3 балла

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

Ti+1=Ti+αΔtj=1nTi,jTixjx2.T_{i+1} = T_{i} + \alpha \Delta t \sum\limits_{j=1}^n \frac{T_{i,j} - T_{i}}{\left|\vec{x}_j-\vec{x}\right|^2}.

Здесь nn — количество соседних тетраэдров, Ti,jT_{i,j} — температура тетраэдра jj в момент времени ii, TiT_{i} — температура текущего тетраэдра в момент времени ii, x\vec{x} — координата центра тетраэдра, α\alpha — коэффициент теплопроводности.

Для выполнения этого задания вам понадобится таблица смежности тетраэдров. Тетраэдры являются смежными, если три любых вершины одного тетраэдра совпадают с тремя любыми вершинами другого.

Видео

2023

  • 00:00 Введение.
  • 01:09 Обозначение производных.
  • 07:18 Вывод с помощью ряда Тейлора.
  • 10:50 Вывод с помощью полиномов.
  • 16:53 Вычисление с помощью преобразования Фурье.
  • 20:49 Производная по направлению.
  • 23:32 Задания.
  • 35:03 Вопросы.
Запись лекции 14.04.2023.

2022

  • 00:00 Введение.
  • 00:40 Обозначение производных.
  • 05:58 Конечные разности.
  • 11:20 Вывод с помощью ряда Тейлора.
  • 18:36 Вывод с помощью полиномов.
  • 39:15 Вычисление с помощью преобразования Фурье.
  • 50:54 Производная по направлению.
  • 56:26 Задания.
Запись лекции 14.04.2022.

Доска

2023

numerical-differentiation-2023 thumbnail.
Запись доски 14.04.2023.

2022

numerical-differentiation thumbnail.
Запись доски 14.04.2022.