6. Исследование погрешностей вычислений

Работая над решением задачи, всегда полезно знать правильный

ответ

Закон Мерфи

Цель работы – исследование вычислительных погрешностей при использовании различных вещественных типов данных и получение практических навыков при программировании сугубо вычислительных алгоритмов. (4 час.)

Методические указания.

Как известно, любые вычисления, в нетривиальных случаях, можно выполнить только с некоторой точностью. Например, результатом операции 1.0/3.0 является 0 целых и бесконечное число "троек". В то же время выполнение операторов

cout.precision(20);

cout<<1./3.;

даст следующий результат:

0.33333333333333331482

Отметим, что для констант 1. и 3. в приведенном выше выражении компилятором выбран тип double по умолчанию, в связи с чем значащих цифр в полученном результате 16, а последние 4 цифры (1482) являются случайными (или псевдослучайными, строго говоря).

Выполнение же оператора

cout<<(float)1./(float)3.;

дает такой результат

0.33333334326744079589

В этом случае явно заданы типы констант 1. и 3. посредством приведения типа, вследствие чего число значащих цифр и сокращается до 7. Еще один поучительный (или умопомрачительный?!) пример - это такой фрагмент программы:

cout.precision(20); // выводить 20 значащих цифр

cout<<0.5; // выводится 0.5

cout<<0.6; // выводится 0.59999999999999998

cout<<0.6f; // выводится 0.60000002384185791

cout<<0.5+1e-20; // выводится 0.5

Или другой пример – вывод значений от -1 до +1 с шагом 0.1 следующим образом (проверьте, что получится, но не забудьте оператор cout.precision(20);):

for(double d=-1.0; d<=1.0; d+=0.1) cout<<d<<endl;

Не правда ли, интересные результаты? Может быть, прикладная математика действительно имеет право на существование, как и специалисты в этой области?

Задание.

В данной работе требуется реализовать вычисление числа Эйлера Е в степени X по формуле:

n

E^X = 1+ ∑ (X^k / k!) (1)

k=1

с максимально возможной точностью для двух вещественных типов: float и double. Как нетрудно видеть, искомый результат – E в степени X – получается путем суммирования членов ряда, причем с ростом k абсолютные значения членов ряда убывают, стремясь к нулю. Так как число значащих цифр для любого типа данных ограничено, то, начиная с некоторого k, прибавление очередного члена ряда уже не будет влиять на результат – E^X. Если мы найдем такой член ряда, то по его величине мы и сможем судить о погрешности вычисления E^X. Нахождение такого члена ряда для каждого из названных вещественных типов данных (float и double) и является Вашей благородной задачей, решив которую Вы прославите себя на века. Для того чтобы убедиться в корректности программы, необходимо вывести на монитор для каждого типа данных:

• значение E^X;

• погрешность вычисления E^X;

• число членов ряда n, при котором достигается максимальная точность;

• значение E^X, полученное с помощью библиотечной функции exp().

Проверьте корректность работы программы при нулевых, положительных, отрицательных и нецелых значениях X. Вычисление E^X для каждого из вещественных типов данных удобно вынести в отдельную функцию. Программу можно оформить в виде консольного приложения.

Для того чтобы с помощью потока cout выводить заданное число значащих цифр выражения вещественного типа, можно воспользоваться компонентной функцией precision() класса iostream

precision(число_значащих_цифр)

следующим образом:

cout.precision(20);

Для первоначальной проверки корректности вычислений напишите программу, которая будет вычислять E^X при заданном n. Посмотрите, будет ли уточняться значение E^X с ростом n.

При X==1 значение E^X должно быть равно приблизительно 2.718281. Можно воспользоваться также библиотечной функцией exp(), объявленной в заголовочном файле Math.h, для сравнения результатов своих вычислений с библиотечными.

Если при вычислении какой-либо величины (вещественного типа) происходит выход за границу диапазона для данного типа, то состояние ошибки не возникает, а значение величины при выводе представляется как 1.#INF (сокращение от infinite – бесконечный). Это характерное значение можно (и нужно!) проконтролировать с помощью функции infinity(), которая в этом случае возвращает значение true. Например:

float fact;

// вычисления

if(fact==numeric_limits<float>::infinity()) cout<<"Караул, переполнение!!!"<<endl;

Для типа double надо использовать функцию таким образом:

if(fact==numeric_limits<double>::infinity())cout<<"Ужас, переполнение!!!"<<endl;

Для использования функции infinity() надо подключить заголовочный файл:

#include <limits>

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

X=709.75

Exp = 1.7398368732641691e+308 eps=7.8270482774736306e+291 n=938

Math::exp = 1.7398368732641605e+308