Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3791 / 1086 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Специальности: Программист, Математик
Лекция 2:

Предмет вычислительной математики. Обусловленность задачи, устойчивость алгоритма, погрешности вычислений. Задача численного дифференцирования

1.2. Влияние выбора вычислительного алгоритма на результаты вычислений

Пример 1.6. Пусть необходимо вычислить значение выражения

$ \left({\frac{{\sqrt{2} - 1}}{{\sqrt{2} + 1}}}\right)^3 . $

Избавившись от знаменателя, получаем (\sqrt{2} - 1)^6  = (3 - 2\sqrt{2})^3  = 99 - 70\sqrt{2}.

Полагая

а)

$ \sqrt{2}  \approx  \frac{7}{5} = 1,4 , $

в)

$ \sqrt{2}  \approx  \frac{17}{12} = 1,41(6) $

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

\sqrt{2} {(\sqrt{2} - 1)}^6 {(3 - 2\sqrt{2})}^3 99 - 70\sqrt{2}
7/5 0,004096 0,008000 1
17/12 0,005233 0,004630 - 0,1(6)

Очевидно, что столь значительное различие в результатах вызвано влиянием ошибки округления в задании \sqrt{2}.

Пример 1.7. Вычисление функции sin x с помощью ряда Тейлора.

Из курса математического анализа известно, что функция синус представляется своим рядом Тейлора

$ \sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +  \ldots , 
$

причем радиус сходимости ряда равен бесконечности — ряд сходится при любых значениях x.

Вычислим значения синуса при двух значениях аргумента. Пусть сначала x = 0,5236 (30^\circ). Будем учитывать лишь члены ряда, большие, чем 10- 4. Выполнив вычисления с четырьмя значащими цифрами, получим sin (0.5236) = 0.5000, что соответствует принятой точности.

Пусть теперь x = 25,66 (1470^\circ). Если вычисления по данной формуле проводить с восемью значащими цифрами, то получим абсурдный результат: sin (25,66) \approx  24 (учитывались члены ряда, большие, чем 10-8 ).

Разумеется, выходом из создавшейся ситуации может быть использование формул приведения.

Пример 1.8. Вычисление функции ex с помощью ряда Тейлора.

Из курса математического анализа известно, что экспонента представляется своим рядом Тейлора

$ e^x  = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} +  \ldots , $

радиус сходимости этого ряда также равен бесконечности.

Приведем некоторые результаты расчетов ( e^x_M — значения экспоненты, вычисленные на компьютере).

x ex exM
1 2,718282 2,718282
20 4,8516520 \cdot 10^8 4,8516531 \cdot 10^8
-1 0,3678795 0,3678794
-10 4,5399930 \cdot10 ^{-5} -1,6408609 \cdot 10^{-4}
-20 2,0611537 \cdot 10^{-9} 1,202966

Выходом из этой ситуации может быть использование для отрицательных аргументов экспоненты формулы

$ e^{- x}  = \frac{1}{e^x } = \frac{1}{1 + x + \frac{x^2}{2!} +  \ldots } .
$

Естественно ожидать рост ошибок округления при вычислении рассматриваемой функции при больших значениях аргумента x. В этом случае можно использовать формулу ex = en + a = enea, где n = [x].

Пример 1.9. Рассмотрим следующий метод вычисления интеграла

I_n = \int\limits_0^1 {x^n}{e^{1-x}dx}, n=0,1,2\ldots

Интегрирование по частям дает

I_n = \int\limits_0^1 {x^n}{d(-e^{1-x})} = - x^{n-1}\cdot e^{1-x} |\limits_0^1 + \int\limits_0^1  e^{1-x}\cdot d(x^n) = -1 + \int\limits_0^1 nx^{n-1}e^{1-x}dx,

откуда следует

I_0  = \int\limits_0^1 e^{1 - x}dx = e - 1  \approx  1,71828.
I_n = nI_{n-1}-1, n \ge 1.

Тогда

\begin{gather*}
I_1 = 1 \cdot I_0 - 1  \approx  0.71828, \quad I_2 = 2I_1 - 1  \approx  0.43656, \quad I_3 = 3I_2 - 1  \approx  0.30968, \\ 
I_4 = 4I_3 - 1  \approx  0.23872, \quad I_5 = 5I_4 - 1  \approx  0.1936, \quad I_6 = 6  \cdot I_5 - 1  \approx  0.16160, \\ 
I_7 = 7I_6 - 1  \approx  0.13120, \quad I_8 = 8I_7 - 1  \approx  0.00496, \quad I_9 = 9I_8 - 1  \approx  - 0.55360, \\ 
{I}_{10} = 10\cdot I_9 - 1  \approx  -6.5360 
\end{gather*}

Очевидно, что отрицательные значения при n = 9,10 не имеют смысла. Дело в том, что ошибка, сделанная при округлении I0 до 6-ти значащих цифр сохранилась при вычислении I1, умножилась на 2! при вычислении I2, на 3! — при вычислении I3, и так далее, т.е. ошибка растет очень быстро, пропорционально n!.

Пример 1.10. Рассмотрим методический пример вычислений на модельном компьютере, обеспечивающем точность \delta_M = 0,0005. Проанализируем причину происхождения ошибки, например, при вычитании двух чисел, взятых с точностью до третьей цифры после десятичной точки u = 1,001, v = 1,002, разность которых составляет \Delta   = |v_{M} - u_{M}| = 0,001.

В памяти машины эти же числа представляются в виде

u_M  = u(1 + \delta_M^u), v_M = v(1 + \delta_M^v), \mbox{ причем } \mid \delta_M^u\mid \le \delta_M \mbox{ и } \mid \delta_M^v\mid \le \delta_M.

Тогда u_M - u  \approx  u\delta_M^u, v_M - v  \approx  v\delta_M^v.

Относительная ошибка при вычислении разности uM - vM будет равна

$ \delta  = \frac{(u_M - v_M) - (u - v)}{(u - v)} = \frac{(u_M  - u) - (v_M  - v)}{(u - v)} = \frac{u\delta_M^u - v\delta_M^v}{(u - v)}. $

Очевидно, что $ \delta  = \left|{\frac{\delta_M^u  - \delta_M^v}{\Delta }} \right| \le \frac{2\delta_M}{0,001}  \approx  2000\delta_M  = 1, $ т.е. все значащие цифры могут оказаться неверными.

Пример 1.11. Рассмотрим рекуррентное соотношение ui+1 = qui, i \ge 0, u0 = a, q > 0, ui > 0.

Пусть при выполнении реальных вычислений с конечной длиной мантиссы на i -м шаге возникла погрешность округления, и вычисления проводятся с возмущенным значением u_i^M  = u_i  + \delta_i, тогда вместо ui+1 получим u_{i + 1}^M  = q(u_i  + \delta_i) = u_{i + 1}  + q\delta_i, т.е. \delta_{i + 1} = q\delta_i,\quad i = 0,1,\ldots.

Следовательно, если |q| > 1, то в процессе вычислений погрешность, связанная с возникшей ошибкой округления, будет возрастать ( алгоритм неустойчив ). В случае \mid q\mid \le 1 погрешность не возрастает и численный алгоритм устойчив.

1.3. Экономичность вычислительного метода

Пример 1.12. Пусть требуется вычислить сумму S = 1 + x + x2 + … + x1023 при 0 < x < 1. Для последовательного вычисления x, x^2 = x \cdot x, \ldots, x^{1023} = x^{1022} \cdot x необходимо проделать 1022 умножения, а затем столько же сложений.

Однако если заметить, что

$ S = \frac{1 - x^{1024}}{1 - x} $
, то количество арифметических действий значительно уменьшается; в частности, для вычисления x1024 требуется всего 10 умножений: x^2 = x \cdot x; x^4 = (x)^2 (x)^2,  \ldots , x^{1024} = (x)^{512}(x)^{512}.

Пример 1.13. Вычисления значений многочленов. Если вычислять значение многочлена P(x) = a0 + a1x + a2x2 + … + anxn "в лоб", т.е. вычислять значения каждого члена и суммировать, то окажется, что необходимо выполнить (n2 + [n/2]) умножений и n сложений. Кроме того, такой способ вычислений может привести к накоплению ошибок округления при вычислениях с плавающей точкой.

Его очевидным улучшением является вычисление каждого члена последовательным умножением на x. Такой алгоритм требует (2n - 1) умножение и n сложений.

Еще более экономичным алгоритмом является хорошо известная в алгебре схема Горнера:

P(x) = ((...((anx + an - 1)x + an - 2)x + ... + a0),

требующая n операций сложения и n операций умножения. Этот метод был известен в средние века в Китае под названием Тянь-Юань и был заново открыт в Европе в начале XIX века англичанином Горнером и итальянцем Руффини.

Пример 1.14. Рассмотрим систему линейных алгебраических уравнений (СЛАУ) вида {\mathbf{Au}} = {\mathbf{f}}, {\mathbf{u}} = \{ u_1 , \ldots , 
u_n \}^T, {\mathbf{f}} = \{ f_1 , \ldots ,f_n \}^T, с трехдиагональной матрицей

\mathbf{A} =
\left( \begin {array}{cccccc}
b_1 & c_1 & 0 & \ldots & 0 \\
a_2 & b_2 & c_2 & \ldots & 0 \\
0 & a_3 & b_3 & c_3 & \ldots & 0 \\
0 & \ldots & 0 & a_n & b_n \\
\end {array} \right).

Если проводить решение такой системы без учета специфической структуры матрицы (например, с помощью метода Гаусса), то количество арифметических действий будет порядка n3, если же учесть эту структуру, то количество операций можно уменьшить до n.