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

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

Я могу вычислить первую производную с помощью:

f = @(x) (x^3);
h = 1e-10;

df = @(x) (f(x+h) - f(x))/h;

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

f = @(x) (x^3);
h = 1e-10;

d2f = @(x) (f(x+h) - 2*f(x) + f(x-h))/(h^2);

Для d2f я должен получить функцию, аналогичную d2f = 6x, но если сюжет d2f, я получаю следующее: сюжет d2f

Что я делаю неправильно?


matlab numerical-methods differentiation
person Ulises    schedule 10.04.2019    source источник
comment
Было бы полезно, если бы вы предоставили информацию о том, что не работает. Что вы пытаетесь, какой результат вы ожидаете и что вы получаете взамен?   -  person Will    schedule 10.04.2019
comment
См. также: stackoverflow.com/questions/22160712/   -  person horchler    schedule 10.04.2019
comment
Проблема здесь с размером h и катастрофической отменой из-за числовой точности с плавающей запятой . Попробуйте h=1e-7 или больше. Однако, тратя на то, где вы оцениваете функции, вы все равно можете увидеть много шума.   -  person horchler    schedule 10.04.2019
comment
Вы также можете увеличить точность с плавающей запятой h, используя vpa: h = vpa(1e-10). Теперь у вас есть точность в 32 цифры, поэтому h должно быть меньше 1e-16.   -  person obchardon    schedule 10.04.2019
comment
@horchler Это сработало! Спасибо!!   -  person Ulises    schedule 10.04.2019
comment
Загляните на en.wikipedia.org/wiki/Numerical_дифференциация и ознакомьтесь с книгой Числовые рецепты. Ваш выбор h слишком мал, что приводит к ошибкам округления и катастрофической отмене.   -  person kvantour    schedule 10.04.2019
comment
Возможный дубликат Попытка реализовать формулу разности в MATLAB   -  person Cris Luengo    schedule 10.04.2019

Ответы (2)


arrow_upward
3
arrow_downward

Формула разделенной разности имеет теоретическую ошибку O (h ^ 2). Вычисление функций с плавающей запятой будет производить относительную ошибку, равную машинной точности mu. Затем это делится на h ^ 2. Наилучшая сумма обеих ошибок достигается там, где они находятся в равновесии, то есть где h^4=mu или h=1e-4.

Это, конечно, неверно, если коэффициент члена ошибки, который является 4-й производной от f, равен нулю, как это происходит с f(x)=x^3. Тогда единственный вклад в ошибку - это ошибки с плавающей запятой, которые наименьшие для больших h, даже h = 1 даст минимальную ошибку.

Для менее тривиальной функции, такой как f (x) = sin (x), ошибка для разных h ведет себя как на следующем графике (где переменная, помеченная x, представляет собой размер шага h)

введите здесь описание изображения

person Lutz Lehmann    schedule 10.04.2019

arrow_upward
0
arrow_downward

Я не уверен, что вы делаете неправильно, но код ниже работает

f=@(x) x.^3;

x = (0:1E-12:1E-6)' ;
d2y = secondDerivative(f,x(1),x(end),x(2)-x(1))';

fit(x,d2y,'poly1')

ans = 

     Linear model Poly1:
     ans(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =           6  (6, 6)
       p2 =   1.352e-15  (-4.734e-13, 4.761e-13)

Определение функции

function d2y = secondDerivative(f, x1, x2, dx)

y = f(x1:dx:x2);

d2y = nan(size(y));
d2y(2:end-1) = y(1:end-2) - 2*y(2:end-1) + y(3:end);

if length(d2y) == 3
    d2y(1) = y(1) - 2*y(2) + y(3);
    d2y(2) = y(end-2) - 2*y(end-1) + y(end);
elseif length(d2y) > 4
    d2y(1) = 2*y(1) - 5*y(2) + 4*y(3) - y(4);
    d2y(end) = -y(end-3) + 4*y(end-2) - 5*y(end-1) + 2*y(end);
end

d2y = d2y / dx^2 ;
end
person user1543042    schedule 11.04.2019