Какой лучший способ вычислить числовую производную в MATLAB?

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

1. У вас есть символическая функция или набор точек?

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

2. Ваша сетка расположена равномерно или неравномерно?

• Если ваша сетка расположена равномерно, вы, вероятно, захотите использовать схему конечных разностей (см. Любую из статей Википедии здесь или здесь), если вы не используете периодическую границу условия (см. ниже). Здесь хорошее введение в конечный разностные методы в контексте решения обыкновенных дифференциальных уравнений на сетке (особенно см. слайды 9-14). Эти методы, как правило, эффективны с точки зрения вычислений, просты в реализации, а ошибку метода можно просто оценить как ошибку усечения разложений Тейлора, использованных для ее получения.
• Если ваша сетка неравномерно разнесена, вы все равно можете использовать схему конечных разностей, но выражения сложнее, а точность очень сильно зависит от того, насколько однородна ваша сетка. Если ваша сетка очень неоднородна, вам, вероятно, потребуется использовать трафареты больших размеров (большее количество соседних точек) для вычисления производной в данной точке. Люди часто строят интерполирующий полином (часто полином Лагранжа) и дифференцируют этот полином для вычисления производной. См., Например, этот вопрос StackExchange. Часто бывает трудно оценить ошибку с помощью этих методов (хотя некоторые пытались это сделать: здесь и здесь). Метод Форнберга часто очень полезен в этих случаях ....
• Следует проявлять осторожность на границах вашего домена, потому что трафарет часто включает точки, которые находятся за пределами домена. Некоторые люди вводят «точки-призраки» или комбинируют граничные условия с производными разного порядка, чтобы устранить эти «точки-призраки» и упростить трафарет. Другой подход - использовать правосторонние или левосторонние методы конечных разностей.
• Вот отличная "шпаргалка" конечно-разностные методы, в том числе центрированные, правые и левосторонние схемы низких порядков. Я храню распечатку этого документа рядом с моей рабочей станцией, потому что считаю его очень полезным.

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

• Если ваш домен периодический, вы можете вычислять производные с очень высокой точностью, используя спектральные методы Фурье. Этот метод несколько жертвует производительностью, чтобы получить высокую точность. Фактически, если вы используете N точек, ваша оценка производной будет примерно N ^ -го порядка точности. Для получения дополнительной информации см. (Например) эту WikiBook.
• В методах Фурье часто используется алгоритм быстрого преобразования Фурье (БПФ) для достижения производительности примерно O (N log (N)), а не алгоритм O (N ^ 2), который может использовать наивно реализованное дискретное преобразование Фурье (ДПФ).
• Если ваша функция и область не периодичны, вам не следует использовать спектральный метод Фурье. Если вы попытаетесь использовать его с функцией, которая не периодическая, вы получите большие ошибки и нежелательные явления "звонка".
• Вычисление производных любого порядка требует 1) преобразования из сеточного пространства в спектральное (O (N log (N))), 2) умножения коэффициентов Фурье на их спектральные волновые числа (O (N)) и 2) обратное преобразование из спектрального пространства в сеточное (снова O (N log (N))).
• Следует соблюдать осторожность при умножении коэффициентов Фурье на их спектральные волновые числа. Кажется, что каждая реализация алгоритма БПФ имеет свой собственный порядок спектральных режимов и параметров нормализации. См., Например, ответ на этот вопрос на Math StackExchange, для заметок об этом в MATLAB.

4. Какой уровень точности вы ищете? Вам нужно вычислить производные с заданным допуском?

• Для многих целей может быть достаточно конечно-разностной схемы 1-го или 2-го порядка. Для большей точности вы можете использовать разложения Тейлора более высокого порядка, отбрасывая члены более высокого порядка.
• Если вам нужно вычислить производные в пределах заданного допуска, вы можете поискать схему высокого порядка, которая имеет нужную вам ошибку.
• Часто лучший способ уменьшить ошибку - это уменьшить шаг сетки в конечно-разностной схеме, но это не всегда возможно.
• Имейте в виду, что схемы конечных разностей более высокого порядка почти всегда требуют большего размера шаблона (большего количества соседних точек). Это может вызвать проблемы на границах. (См. Обсуждение точек-призраков выше.)

5. Имеет ли значение для вас, что ваша производная оценивается в тех же точках, что и ваша функция?

• MATLAB предоставляет функцию diff для вычисления различий между соседними элементами массива. Это можно использовать для вычисления приближенных производных с помощью схемы прямого дифференцирования (или прямой конечной разности), но оценки являются оценками низкого порядка. Как описано в документации MATLAB для diff (ссылка), если вы введете массив длины N, он вернет массив длины N-1. Когда вы оцениваете производные с помощью этого метода по N точкам, вы будете иметь оценки производной только в N-1 точках. (Обратите внимание, что это можно использовать для неравномерных сеток, если они отсортированы в порядке возрастания.)
• В большинстве случаев мы хотим, чтобы производная оценивалась во всех точках, а это означает, что мы хотим использовать что-то помимо метода diff.

6. Вам нужно рассчитать несколько заказов деривативов?

• Можно составить систему уравнений, в которой значения функции точки сетки и производные 1-го и 2-го порядка в этих точках зависят друг от друга. Это можно найти, комбинируя разложения Тейлора в соседних точках, как обычно, но сохраняя производные члены, а не сокращая их, и связывая их вместе с членами соседних точек. Эти уравнения можно решить с помощью линейной алгебры, чтобы получить не только первую производную, но и вторую (или более высокие порядки, если они настроены правильно). Я считаю, что это называется комбинированными конечно-разностными схемами, и они часто используются вместе с компактными конечно-разностными схемами, которые будут обсуждаться далее.
• Компактные конечно-разностные схемы (ссылка). В этих схемах создается матрица плана и производные вычисляются во всех точках одновременно с помощью матричного решения. Их называют «компактными», потому что они обычно проектируются так, чтобы требовать меньшего количества точек трафарета, чем обычные конечно-разностные схемы сопоставимой точности. Поскольку они включают матричное уравнение, связывающее все точки вместе, некоторые компактные конечно-разностные схемы имеют «спектральное разрешение» (например, Статья Леле 1992 года - отлично!), что означает, что они имитируют спектральные схемы, зависящие от всех узловых значений и, из-за этого , они сохраняют точность на всех масштабах длины. Напротив, типичные методы конечных разностей точны только локально (производная в точке №13, например, обычно не зависит от значения функции в точке №200).
• Текущая область исследований - как лучше всего найти несколько производных в компактном трафарете. Результаты таких исследований, комбинированные, компактные методы конечных разностей, являются мощными и широко применимыми, хотя многие исследователи склонны настраивать их для конкретных потребностей (производительность, точность, стабильность или конкретная область исследования, например как гидродинамика).

Готовые программы

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

• Подпрограмма MATLAB gradient (ссылка на документацию) - отличный вариант для многих целей. . Он реализует центральную разностную схему второго порядка. Он имеет преимущества вычисления производных в нескольких измерениях и поддержки произвольного шага сетки. (Спасибо @thewaywewalk за указание на это вопиющее упущение!)

• Я использовал метод Форнберга (см. Выше), чтобы разработать небольшую процедуру (nderiv_fornberg) для вычисления конечных разностей в одном измерении для произвольных интервалов сетки. Мне легко пользоваться. Он использует двусторонние трафареты с 6 точками по границам и центрированный 5-точечный трафарет в интерьере. Он доступен в разделе обмена файлами MATLAB здесь.

Заключение

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

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

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