C++ - вторая производная.

[Katmandu]

Необходимо посчитать значение второй производной в точке.

Я так понимаю что

$$ f''(x) =\lim_{\Delta x \rightarrow 0}{ \frac {f'(x+\Delta x) - f'(x)} {\Delta x} }, $$

Например исследовал для функции $ f(x) = x^2 + x + 1 $

Написал программу.
Первую производную считает всегда почти правильно.
А со второй проблемы.
При $ EPS = 0.01 $ был вывод на экран: $ 43 \quad 13.01\quad 2$
При $ EPS = 0.000001 $ был вывод на экран: $ 43\quad 13\quad 2.01084$
При $ EPS = 0.0000001$ вывод: $ 43\quad 13\quad 2.13163$
При $ EPS = 0.00000001$ вывод: $ 43\quad 13\quad -71.0543$ , -71 то откуда
При $ EPS = 0.000000001$ вывод: $ 43\quad 13\quad -7105.43$ ,
При $ EPS = 0.0000000001$ вывод: $43\quad 13\quad 0$ ,
При $ EPS = 0.00000000001$ вывод: $43\quad 13.0001\quad 0$
При $ EPS = 0.00000000000001$ вывод: $43\quad 12.7898\quad 0$

Я понимаю, конечно, вычислительная погрешность там, округления при вычислении.
Но почему при точности $0.01$ он посчитал точнее, чем при $0.00000001$



Код программы:

Используется синтаксис C++

 
#define EPS 0.00000001

long double diff(long double (*f)(long double), long double);
long double diff2(long double (*f)(long double),long double X);

long double MNS(long double (*f)(long double))
{
        std::cout << f(6) << " " << diff(f,6) << " " << diff2(f,6) << std::endl;
        return 0;
}

long double diff(long double (*f)(long double),long double X)
{
        return ((f(X+EPS)-f(X))/EPS);
}
        
long double diff2(long double (*f)(long double),long double X)
{
        return (diff(f,X+EPS) - diff(f,X))/EPS;
}

-- 27.04.2014, 14:51 --

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

Используется синтаксис C++

long double diff(long double (*f)(long double),long double X)
{
        return (f(X+EPS)-f(X-EPS))/(2*EPS);
}
        
long double diff2(long double (*f)(long double),long double X)
{
        return (diff(f,X+EPS) - diff(f,X-EPS))/(2*EPS);
}
 

[ewert]

Katmandu в сообщении #855679 писал(а):

Я понимаю, конечно, вычислительная погрешность там, округления при вычислении.
Но почему при точности $0.01$ он посчитал точнее, чем при $0.00000001$

Он и должен был посчитать точнее. У double запас -- где-то в 15-16 знаков. Грубо говоря, это означает, что абсолютная погрешность округления для числителя -- порядка десяти в минус четырнадцатой. И потом Вы её делите на десять в минус шестнадцатой. Вот и получается погрешность окончательного результата порядка сотни. А Вы чего ожидали?...

[EtCetera]

ewert в сообщении #855694 писал(а):

У double запас -- где-то в 15-16 знаков.

Справедливости ради, ТС использует long double, который в большинстве случаев реализован как 80-битное число с плавающей точкой, т.е. у мантиссы имеется порядка 19 значащих десятичных цифр.

Я бы сказал проще. Вычитание вещественных чисел одного знака и сложение чисел разных знаков $\text{---}$ это самые неприятные операции в компьютерной арифметике. Погрешности в случае их использования для операндов с близкими абсолютными значениями могут быть катастрофическими.

[ewert]

EtCetera в сообщении #855707 писал(а):

ТС использует long double, который в большинстве случаев реализован как 80-битное число с плавающей точкой, т.е. у мантиссы имеется порядка 19 значащих десятичных цифр.

Я на это не обратил внимания. Там действительно на два-три порядка лучше, и семидесяти быть действительно не должно -- в худшем случае порядка единицы. Единственно, что могу предположить -- что где-то в промежутке вклинивается обычный восьмибайтовый тип (ну, скажем, в эпсилоне).

EtCetera в сообщении #855707 писал(а):

Я бы сказал проще.

Я бы не стал так говорить. Вопрос о том, при каком шаге достигается наилучшая точность разностного дифференцирования -- вполне содержателен, причём с сугубо прикладной точки зрения.

[iifat]

80 бит — это размерность внутренних регистров сопроцессора. long double — таки 64, насколько мне известно.

[ewert]

iifat в сообщении #855714 писал(а):

long double — таки 64, насколько мне известно.

Ну я в сях не разбираюсь, но вообще-то 64 -- это везде просто double. Такой в определённом смысле стандарт де-факто.

[EtCetera]

ewert

ewert в сообщении #855711 писал(а):

где-то в промежутке вклинивается обычный восьмибайтовый тип (ну, скажем, в эпсилоне)

Да, 0.00000001 по стандарту C++ имеет тип обычного double, но вычисления, по-моему, все равно должны производиться в 80 битах $\text{---}$ формате более "мощного" операнда. Можно порекомендовать ТС исправить 0.00000001 на 0.00000001L (а еще лучше на 1e-8L, для большей читабельности). А также заглянуть в настройки компилятора. И перестать использовать макросы в качестве констант!

iifat

iifat в сообщении #855714 писал(а):

long double — таки 64, насколько мне известно.

В большинстве случаев это не так.