Как оценить производную функции в Matlab?

Question

Это должно быть очень просто. У меня есть функция f(x), и я хочу оценить f'(x) для данного x в MATLAB.

Все мои поиски были связаны с символической математикой, а это не то, что мне нужно, мне нужно численное дифференцирование.

например. если я определю: fx = inline('x.^2')

Я хочу найти, скажем f'(3), что будет 6, я не хочу найти 2x

Answer 1

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

f'(x) = (f(x + h) - f(x - h)) / 2h

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

f'(x) = (f(x + h) - f(x)) / h     (*)

это первый порядок в ч.

Это теория. На практике все довольно сложно. Я возьму вторую формулу (первый порядок), так как анализ проще. Выполните второй заказ как упражнение.

Самое первое наблюдение состоит в том, что вы должны убедиться, что (x + h) - x = h, иначе вы получите огромные ошибки. Действительно, f (x + h) и f (x) близки друг к другу (скажем, 2.0456 и 2.0467), и когда вы вычитаете их, вы теряете много значащих цифр (здесь это 0,0011, что на 3 значимых цифры меньше чем х). Поэтому любая ошибка в h может оказать огромное влияние на результат.

Итак, первый шаг, исправьте кандидата h (я через минуту покажу вам, как его выбрать), и возьмите в качестве значения h для вашего вычисления величину h '= (x + h) - x. Если вы используете такой язык, как C, вы должны позаботиться о том, чтобы определить h или x как изменчивые, чтобы эти вычисления не оптимизировались.

Далее выбор ч. Ошибка в (*) состоит из двух частей: ошибка усечения и ошибка округления. Ошибка усечения вызвана тем, что формула не является точной:

(f(x + h) - f(x)) / h = f'(x) + e1(h)

, где e1(h) = h / 2 * sup_{x in [0,h]} |f''(x)|.

Ошибка округления связана с тем, что f (x + h) и f (x) близки друг к другу. Это можно оценить примерно как

e2(h) ~ epsilon_f |f(x) / h|

, где epsilon_f - относительная точность вычисления f (x) (или f (x + h), что близко). Это должно быть оценено по вашей проблеме. Для простых функций epsilon_f можно принять за эпсилон машины. Для более сложных это может быть хуже, чем это на порядки.

Итак, вы хотите h, который минимизирует e1(h) + e2(h). Объединяя все вместе и оптимизируя в h, получаем

h ~ sqrt(2 * epsilon_f * f / f'')

, который должен оцениваться по вашей функции. Вы можете принять приблизительные оценки. Если сомневаетесь, возьмите h ~ sqrt (epsilon), где epsilon = точность машины. Для оптимального выбора h относительная точность, с которой известна производная, равна sqrt (epsilon_f), т.е. половина значащих цифр верна.

Вкратце: слишком маленькая ошибка h => округления, слишком большая ошибка усечения h =>.

Для формулы второго порядка те же вычисления дают

h ~ (6 * epsilon_f / f''')^(1/3)

и дробная точность (epsilon_f) ^ (2/3) для производной (которая, как правило, на одну или две значащие цифры лучше, чем формула первого порядка, при условии двойной точности).

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

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

Answer 2

Чтобы получить числовую разницу (симметричная разница), вы вычисляете (f(x+dx)-f(x-dx))/(2*dx)

fx = @(x)x.^2;
fPrimeAt3 = (fx(3.1)-fx(2.9))/0.2;

В качестве альтернативы вы можете создать вектор значений функции и применить DIFF , т.е.

xValues = 2:0.1:4;
fValues = fx(xValues);
df = diff(fValues)./0.1; 

Обратите внимание, что diff принимает прямую разницу и что предполагается, что dx равно 1.

Однако в вашем случае вам лучше определить fx как полином и оценить производную функции, а не значения функции.

Answer 3

Не имея символического набора инструментов, ничто не мешает вам использовать Derivest , инструмент для автоматического адаптивного числового дифференцирования.

derivest(@sin,pi)
ans =
           -1

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

fx = inline('x.^2');
[fp,errest] = derivest(fx,3)

fp =
            6
errest =
   3.6308e-14

Answer 4

пробовали ли вы diff (вычисляет разности и аппроксимирует производную), gradient или polyder (вычисляет производную полинома) функции?

Подробнее об этих функциях можно прочитать виспользуйте help <commandname> на консоли MATLAB или используйте браузер функций в меню Справка.

Answer 5

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

syms x
df = diff(x^2);
df3 = subs(df, 'x', 3);
fprintf('f''(3)=%f\n', df3);

Для чистых числовых производных используйте уже приведенные решения Йонаса и Посдефа.