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

Question

Я хочу вычислить вторую производную анонимной функции в 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

Что я делаю не так?

Answer 1

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

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

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

Answer 2

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

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