17. Кое-что о вещественных вычислениях

В отличие от целочисленных выражений, которые всегда вычисляются точ­но, вещественные выражения дают приближенный результат, и вещественные переменные содержат приближенные значения. Например, выполним программу

VAR x : Real; BEGIN x:=1/3; WRITELN(x*3-1); END. Мы получим не 0 (точное значение выводимого выражения), а, например, 4.547E-13, а может быть, какое-нибудь другое маленькое число. Это обусловлено тем, что переменные типа Real хранят конечное число десятичных цифр (11-12 цифр), кроме того, эти цифры хранятся в двоичном коде, поэтому мы и не полу­чили 1E-12 или 1E-13. Таким образом, x/a*a далеко не всегда равно x. И наобо­рот, a+x может быть равно a, даже если x не равно нулю. Найдем такое положи­тельное число, которое удовлетворяет уравнению x+1=1 :

CONST x:Real=1; BEGIN WHILE x+1<>1 DO x:=x/2; WRITELN(x); END. Мы получим, например, 5.42E-20 (результат зависит от типа компьютера).

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

: вычислить сумму ряда XxVi"! , i=0,1,..., 8. Несмотря на то, что необходимо про­суммировать бесконечное число слагаемых, эта задача легко решается за конеч­ное время, так как общий член ряда быстро убывает и начиная с некоторого n прибавление очередного слагаемого уже не будет изменять сумму. Сначала на­пишем плохую программу:

FUNCTION Factorial(n:Word):Real;

VAR i:Word; f:Real;

BEGIN f:=1; FOR i:=1 TO n DO f:=f*i; Factorial:=f; END;

FUNCTION Power(x:Real; n:Woed):Real;

VAR i:Word; f:Real;

BEGIN IF n=0 THEN Power:=1

ELSE BEGIN f:=x; FOR i:=2 TO n DO f:=f*x; Power:=f; END;

END;

VAR x,S1,S2 : Real;

CONST i : Word = 0;

BEGIN WRITE('Bведите x '); READ(x); S2:=0;

REPEAT S1:=S2; S2:=S1+Power(x,i)/Factorial(i); INC(i); UNTIL S1=S2; WRITELN('Cумма ряда = ',S1);

END.

Запустим эту программу, задав x=1; мы получим верный результат 2.71828... (его легко проверить, поскольку сумма нашего ряда равна exp(x)). А теперь вве­дем x=10, ожидая получить exp(10)=22026.4... Однако ничего подобного не слу­чилось! Turbo Pascal выбросил нас обратно в среду, выдал сообщение "Error 205: Floating point overflow" и показал курсором строку

FOR i:=1 TO n DO f:=f*i; Это аварийное прерывание - переполнение, оно происходит всякий раз, когда ве­щественная величина превышает максимально допустимое значение 1.7E38. Следовательно, для некоторого i мы уже не можем вычислить i! .

Означает ли это, что невозможно вычислить exp(10)? Вовсе нет; конечно, мы не можем задать нашей программе очень большой x, но значение 10 вполне при­емлемо, дело здесь в качестве нашей программы. Действительно, посмотрим, как работает программа : для i=0 вычисляется x0 и 0!, затем для i=1 заново вычисля­ется x1 и 1! и т.д. до получения результата; но xl+1=x*xl и (i+1)!=(i+1)*i!, так что, зная предыдущие значения, достаточно выполнить всего одну операцию, чтобы получить последующие. Более того, нам вовсе не нужен факториал сам по себе, а только общий член ряда (в котором этот факториал находится в знаменателе). Нетрудно записать рекуррентную формулу для общего члена ряда : a0=1;

ai=x/i*ai-1 .Кроме того, что избавимся от переполнения, пользуясь этой формулой, мы во много раз увеличим скорость нашей программы (не всегда функции быва­ют полезны).

VAR x,S1,S2,a : Real; CONST i : Word = 0;

BEGIN WRITE('Введите x '); READ(x); a:=1; S2:=0;

REPEAT S1:=S2; S2:=S1+a; INC(i); a:=a*x/i; UNTIL S1=S2; WRITELN('Сумма ряда = ',S1);

END.

Программа сработала для x=10 и x=20 и x=50, но для x=100 снова произошло переполнение. Но здесь уже ничего сделать нельзя, exp(100)>1043 и никак не мо­жет быть представлена вещественным значением типа Real. Решим еще одну, более содержательную, задачу: найти корень уравнения f(x)=0 методом бисекции. Метод бисекции заключается в следующем: пусть уравнение f(x)=0 имеет един­ственный корень на отрезке [a,b] , это значит, что график функции один раз пе­ресекает ось X на этом отрезке. Определим знак функции в точке a и в точке x=(a+b)/2. Если эти знаки одинаковы, то, очевидно, корень лежит на отрезке [x,b] , в противном случае - на отрезке [a,x] . Таким образом, за один шаг метода мы ровно вдвое уменьшили наш отрезок; будем повторять эти операции до тех пор, пока отрезок не станет очень маленьким, и в качестве корня возьмем сере­дину этого маленького отрезка. Попробуем реализовать этот метод:

CONST a : Real=0; b : Real=10;

FUNCTION f(x:Real):Real; BEGIN f:=exp(x)-2; END;

CONST epsilon=1E-10; {"очень маленький" отрезок}

VAR x,Fa,Fx : Real;

BEGIN Fa:=f(a); {нет необходимости вычислять f(a) в цикле, т.к.ЗНАК

функции на левом конце отрезка не изменится! }

WHILE b-a>epsilon DO BEGIN

x:=(a+b)/2; Fx:=f(x);

IF Fx=0 THEN BEGIN WRITELN(x); Halt; END; IF Fa*Fx<0 THEN b:=x ELSE a:=x;

END;

WRITELN(x);

END.

Программа вычислила 0.693147. Легко проверить, что это правильный ре­зультат (корень уравнения равен ln2). Казалось бы, мы написали хорошую про­грамму, но давайте теперь вычислим корень уравнения ln(x)-50=0, a=1, b=1e30 -программа зациклится ! Выведем внутри цикла значения a и b : эти числа почти одинаковы (отличие в 12-й цифре) и не меняются, но поскольку их порядок 10 , b-a существенно превосходит наш "очень маленький" отрезок. Есть два способа, которыми мы можем исправить программу. Первый заключается в правильном подборе "очень маленького" отрезка, но надо понимать, что вам придется подби­рать этот отрезок для каждого нового уравнения, то есть фактически для каждо­го уравнения писать свою программу. Очевидно, что это бесперспективный путь. Выведем в нашей зацикленной программе не только a и b, но и x, может быть,это поможет нам придумать второй способ: значение x, оказывается, в точности равно b.

Мы могли бы прийти к выводу, что рано или поздно x станет равным или a или b, рассуждая чисто теоретически. Действительно, на каждом шаге цикла мы уменьшаем отрезок в два раза; если бы мы работали на вещественной оси, то ве­личины a и b стремились бы друг к другу бесконечно, но, поскольку множество вещественных чисел в компьютере дискретно (из-за конечного числа цифр), на­станет момент, когда между a и b больше не будет ни одного числа. После этого выражение (a+b)/2 будет давать либо a, либо b. Воспользуемся этим обстоятель­ством и напишем хорошую программу: CONST a : Real = 1; b: Real = 1E30; FUNCTION f(x:Real):Real; BEGIN f:=Ln(x)-50; END;

VAR x,Fa,Fx : Real;

BEGIN Fa:=f(a); x:=(a+b)/2;

WHILE (x<>a)AND(x<>b) DO BEGIN

Fx:=f(x);

IF Fx=0 THEN BEGIN WRITELN(x); Halt; END;

IF Fa*Fx<0 THEN b:=x ELSE a:=x; x:=(a+b)/2;

END;

WRITELN(x);

END.

Программа дала верный результат 5.184705...E21.

Решим еще одну задачу: вычислить значения функции

f(x)=ln(1+ln(1+exp(exp(x)))) на отрезке [0,1000] с шагом 5. CONST x0 = 0; x1 = 1000; h = 5;

FUNCTION f(x:Real):Real; BEGIN f:=ln(1+ln(1+exp(exp(x)))); END;

VAR i : Byte; x : Real;

BEGIN FOR i:=0 TO Round((x1-x0)/h) DO BEGIN

x:=x0+i*h; WRITELN('x=',x:4:0,' f(x)=',f(x)); END;

END.

При x=10 произошло переполнение. Означает ли это, что задача неразреши­ма? Нет, мы просто написали плохую программу, скопировав математическую формулу в оператор Паскаля. Посмотрим, в каком месте происходит переполне­ние - очевидно, при вычислении exp(exp(x)) , других возможностей просто не существует. Это значит, что полученное значение exp(exp(x)) превосходит 1E38. Посмотрим на аргумент внутреннего логарифма: прибавление единицы к очень большому числу никак не изменит это число, следовательно, этой единицей можно пренебречь.Таким образом, для x>5 наша формула упрощается: f(x)=ln(1+ln(1+exp(exp(x))))=ln(1+ln(exp(exp(x))))=ln(1+exp(x))

Исправим программу: FUNCTION f(x:Real):Real;

BEGIN IF x<5 THEN f:=ln(1+ln(1+exp(exp(x)))) ELSE f:=ln(1+exp(x));

Снова произошло переполнение, но теперь при x=90. Очевидно, что причина та же - не удалось вычислить exp(x). Еще раз упростим формулу для x>90 : f(x)=ln(1+exp(x))=ln(exp(x))=x - и запишем, наконец, верную программу: FUNCTION f(x:Real):Real;

BEGIN IF x<5 THEN f:=ln(1+ln(1+exp(exp(x)))) ELSE IF x<90 THEN f:=ln(1+exp(x)) ELSE f:=x;

END;

Сделаем некоторые выводы из вышесказанного. Компьютерные вычисления несколько отличаются от абстрактных математических вычислений. В математи­ке вещественная ось непрерывна (между двумя любыми вещественными числами находится бесконечное множество чисел) - компьютерное множество вещест­венных чисел дискретно. Математика оперирует с бесконечно большими и бес­конечно малыми величинами - компьютерные вещественные числа ограничены сверху и снизу. Математические вычисления точны - компьютерные вычисления приближенны. Вы должны учитывать это, когда программируете какую-либо вычислительную задачу.