Что такое вычислительная математика

• Исходный вариант статьи (И. Б. Петров, Что такое вычислительная математика?) опубликован в журнале «Потенциал».

Аннотация

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

Введение

На первый взгляд кажется, что вычислительная математика — это та же математика, только с вычислениями конкретных величин в формулах. Иногда приходится слышать иную позицию: «Вычислительная математика? Это, кажется, что-то от информатики?» Однако и первая, и вторая точка зрения верны лишь отчасти. Да, у вычислительной математики много общего с математическими дисциплинами; порой говорят о том, что это одна из математических дисциплин.

Вычислительная математика — это наука о методах решения вычислительных задач на компьютере. Она появилась от необходимости решать практические задачи, такие, как управление сложными технологическими процессами, управление полётом ракет, моделирование физических процессов (процесса ядерного распада, химических реакций, роста кристаллов и др.), а не только создание всё более реалистичных и ярких компьютерных игрушек, как ныне думают многие школьники.

Задачами вычислительной математики занимались такие выдающиеся учёные, как Эйлер, Лагранж, Чебышёв, Якоби, Лежандр, фон Нейман и многие другие. Они, часто занимаясь сложными вычислениями вручную на бумаге, невольно заложили основы науки об эффективных безошибочных вычислениях на компьютерах.

Всем известно, что компьютеры имеют дело с числами с ограниченным количеством знаков после запятой^[1]. Казалось бы, какая мелочь! Однако именно эти «мелочи» могут сильно исказить результаты численных расчётов. Появился важнейший раздел вычислительной математики — теория устойчивости вычислительных методов, то есть таких методов, которые позволяют на компьютере с «неточной» арифметикой получать точные (правдивые) результаты.

Как вычислительная математика решает прикладные задачи

Во многих прикладных задачах искомым результатом является не одно число, а целая функция на отрезке. А как описать произвольную функцию на компьютере? Самый простой способ — это хранить значения функции в нескольких точках и «в уме» соединять их гладкой кривой. Таким образом, отрезок может рассматриваться как система нескольких выбранных точек $\left\{{x_{k}}\right\}_{k=0}^{K}\,\!$, а непрерывная функция — как набор $\left\{{f_{k}}\right\}_{k=0}^{K}\,\!$ значений функции в этих точках. Очень важный объект $f'(x)\,\!$ — производная функции в точке $x$ — характеризует угол наклона прямой, касающейся графика функции. Вместо него приходится использовать отношение ${{f_{k+1}-f_{k}} \over {x_{k+1}-x_{k}}}\,\!$, где $x_{k}\,\!$ и $x_{k+1}\,\!$ ближайшие к $x\,\!$ выбранные точки. При решении задач на компьютерах приходится приближать не только числа, но и сами задачи несколько «округлять» и вместо идеальных непрерывных объектов рассматривать их дискретные приближения. Этот вынужденный шаг называется аппроксимацией задачи.

Как узнать, насколько адекватные результаты даёт компьютерная модель, которую вы построили, или, как принято говорить у вычислительных математиков, имеется ли сходимость решения модели к решению исходной «непрерывной» задачи? Кроме ряда сложных и красивых теорем, есть несколько простых хитростей, которые позволяют определить адекватность вашей модели:

1. Точность вычислений. В машинных вычислениях всегда присутствует машинная погрешность, о чём уже упоминалось. И прежде чем верить результатам, полученным на компьютере, попробуйте увеличить точность всех расчетов на пару разрядов. Если это изменит результат в предыдущих знаках, значит, вычисления проводились с недостаточной точностью и нужно использовать более точные типы для представления чисел (например, в языке Си float заменить на double, а в языке Паскаль тип real заменить на extended).
2. Обусловленность. В вычислительной практике большое значение имеет чувствительность решения к малым изменениям входных данных. Задача называется плохо обусловленной, если малые изменения входных данных приводят к заметным изменениям решения. Измерить эту обусловленность на практике очень просто. Нужно просто попробовать чуть-чуть поменять входные данные и посмотреть, как меняется результат. Собственно, нужно выяснить, с какой погрешностью заданы входные данные, и экспериментально проверить в каких пределах меняется результат при варьировании входных данных в пределах их погрешности.
3. Зависимость от алгоритма и модели. Есть ещё один источник неадекватности численных результатов, полученных на компьютере. Это неточность выбранной модели и алгоритма вычисления. Здесь всё несколько сложнее, нежели в предыдущих пунктах. Но для богатых заказчиков можно предложить надёжное решение — следует дать одну и ту же задачу нескольким группам прикладных математиков. Они, скорее всего, построят разные модели и будут пользоваться разными алгоритмами. Если все выдадут один и тот же результат, значит задача хорошая — она устойчива к выбору модели и алгоритма, и не так важно каким способом её решать. Конечно, бывают и такие случаи, когда все делают одну и ту же ошибку. У начинающих вычислительных математиков есть несколько таких излюбленных ошибок.
   

Разумеется, если есть возможность сравнить результаты расчета с экспериментальными, то такие сравнения являются поводом для подтверждения верности расчетов. Что же касается исследования численных методов на аппроксимацию — то эта теория достаточно хорошо разработана, однако без вычислительных ошибок численных решений не бывает. Вопрос заключается в том, значительны ли они или незначительны, что, впрочем, также можно оценить.

Приведём некоторые примеры, демонстрирующие что «шутки с погрешностью плохи».

Пример 1: корни многочлена

Есть множество программ, позволяющие численно вычислять корни многочленов. Давайте рассмотрим уравнение $x^{4}-4x^{3}+8x^{2}-16x+15.\underbrace {9...9} _{8}=(x-2)^{4}-10^{-8}=0;\,\!$

Это уравнение просто решить аналитически:$(x-2)^{2}=10^{-4}\,\!$ , $x_{1,2}=2\pm 0.01;{\rm {}}x_{3,4}=2\pm 0,01i.\,\!$

Предположим, что мы проводим вычисления на компьютере, имеющем точность $\delta _{M}>10^{-8}\,\!$ . В этом случае последнее слагаемое в уравнение будет округлено до 16.0, и с позиции компьютера будет решаться уравнение $\left({x-2}\right)^{4}=0\,\!$ , у которого один корень 2. Таким образом, малые погрешности в задании числа порядка $10^{-8}\,\!$ привели к появлению погрешностей в определении корней уравнения порядка $10^{-2}\,\!$. Более того, вместо двух действительных и двух комплексных корней у нас получился лишь один действительный. Такие плохо обусловленные многочлены не редкость. Например, если в многочлене $\left({x-1}\right)\left({x-2}\right)\ldots \left({x-20}\right)=x^{20}-210x^{19}+\ldots \,\!$

коэффициент при $x^{19}\,\!$ изменить на малую величину порядка $10^{-7}\,\!$ , то вместо максимального корня $x_{20}=20\,\!$ получим корень $x_{20}\approx 20,8\,\!$

Пример 2: cистема уравнений

Другой пример плохо обусловленной задачи — это система двух линейных уравнений: $\left\{{\begin{matrix}u+10v=11,\\100u+1001v=1101.\end{matrix}}\right.\,\!$

Решением является пара чисел $\left\{{1;1}\right\}.\,\!$

«Возмутим» правую часть первого уравнения на 0,01 (вместо 11 напишем 11,01) и получим новую, «возмущенную» систему, решением которой является пара чисел {11,01; 0,00}, не имеющая ничего общего с решением невозмущенной системы. Здесь изменение значения одного параметра на $0,1\%\,\!$ привело к совсем другому решению.

Пример 3: вычисление sin x

В математическом анализе показывается, что функция $\sin x\,\!$ может быть представлена в виде ряда $\sin x=x-{{x^{3}} \over {3!}}+{{x^{5}} \over {5!}}-{{x^{7}} \over {7!}}+...\,\!$, причём этот ряд будет сходиться к значению $\sin x\,\!$ при любых значениях x. После вычисления значения нашей функции при $x{\rm {=0}}{\rm {,5236}}\,\!$ ($x\approx 30^{\circ }\,\!$) с четырьмя значащими цифрами и с учётом лишь членов ряда, больших ${10^{-4}}\,\!$, получим значение $\sin x\,\!$ с той же точностью: $\sin(0,5236)=0,5000.\,\!$

Теперь положим, что $x{\rm {=25}}{\rm {,66}}\,\!$ ($x\approx 1470^{\circ }\,\!$). Если мы проведём вычисления с восемью значащими цифрами с учётом членов ряда, больших ${10^{-8}}\,\!$, то получим ошеломляющий результат; $\sin(25,66)=24,...\,\!$

Впрочем, в данном случае выход из положения достаточно прост — необходимо использовать формулы приведения: $\sin 1470^{\circ }=\sin 30^{\circ }\,\!$. А что делать, если нам необходимо вычислить экспоненту, представленную в виде ряда $e^{x}=1+x+{{x^{2}} \over {2!}}+{{x^{3}} \over {3!}}+...\,\!$? Очевидно, что нас и в этом случае ждут подобные сюрпризы. Однако и здесь существует выход: будем вычислять экспоненту с отрицательным показателем $e^{-x}=1-x+{{x^{2}} \over {2!}}-{{x^{3}} \over {3!}}+...\,\!$ , а затем вычислим обратную величину: $e^{x}={1 \over {e^{-x}}}={1 \over {1-x+{{x^{2}} \over {2!}}-{{x^{3}} \over {3!}}+...}}\,\!$ . Дело в том, что $e^{-x}\,\!$ вычисляется как сумма знакопеременного ряда. После того, когда каждое следующее слагаемое в этой сумме станет меньше, чем предыдущее, по модулю, мы будем чётко представлять верхнюю и нижнюю оценку вычисляемого значения и будем знать, на каком слагаемом следует оборвать суммируемый ряд.

Пример 4: корень из двойки

Мы говорили о том, что на результат вычисления может оказать влияние выбор алгоритма решения задачи. Действительно, пусть нам нужно вычислить значение выражения $A=\left({{{\sqrt {2}}-1} \over {{\sqrt {2}}+1}}\right)^{3}.\,\!$ После избавления от иррациональности в знаменателе, получим: $A=\left({{\sqrt {2}}-1}\right)^{6}\,\!$ или $A=\left({3-2{\sqrt {2}}}\right)^{3}\,\!$ или $A=99-70{\sqrt {2}}.\,\!$ Эти три формулы и будут нашими тремя алгоритмами вычисления числа A. Выберем два разных приближения корня из двойки. ${\sqrt {2}}\approx {7 \over 5}=1,4\,\!$ и ${\sqrt {2}}={{17} \over {12}}\approx 1,416666...\,\!$

Проведём вычисления A по трём формулам для двух приближенных значений ${\sqrt {2}}\,\!$ . Результаты вычислений представим в виде таблицы:

	${\sqrt {2}}\,\!$	$\left({{\sqrt {2}}-1}\right)^{6}\,\!$	$\left({3-2{\sqrt {2}}}\right)^{3}\,\!$	$99-70{\sqrt {2}}\,\!$
1-й способ	7/5	0,004096	0,008000	1
2-й способ	17/12	0,005233	0,004630	-0,1(6)

Видно, что первый «алгоритм» $A=\left({{\sqrt {2}}-1}\right)^{6}\,\!$ является более устойчивым к погрешностям входных данных.

Задачи для самостоятельного решения

Читателю предлагается добраться до компьютера и компилятора какого-нибудь языка программирования (Си или Pascal) и провести несколько незатейливых численных экспериментов.

Задача 1 (частичная сумма гармонического ряда).

Напишите программу, которая суммирует первые миллион слагаемых гармонического ряда сначала с первого по последний элемент, а потом наоборот — с последнего по первый: $A={1 \over 1}+{1 \over 2}+{1 \over 3}+\ldots +{1 \over {10^{6}}}={1 \over {10^{6}}}+\ldots {1 \over 3}+{1 \over 2}+{1 \over 1}\,\!$. Убедитесь, что ассоциативный закон «$(a+b)+c=a+(b+c)\,\!$» при вычислении на компьютере с конечной точностью не выполняется, что может привести к довольно сильному расхождению между результатом вычислений и реальным значением выражения.

Задача 2 («страшный» предел).

Найдите с помощью компьютера значение выражения ${\frac {\sin(\tan(x))-\tan(\sin(x))}{\arcsin(\arctan(x))-\arctan(\arcsin(x))}}\,\!$

при x, равном 1, 1/2, 1/4, 1/8, 1/10, 1/100, … (x измеряется в радианах). Как вы думаете, к чему это выражение на самом деле стремится по мере того, как x стремится к 0?

Задача 3 («башня» из степеней).

Рассмотрим функцию $f(x)=x^{x^{x^{\cdot ^{\cdot ^{\cdot }}}}}\,\!$. Эта функция представляет собой бесконечную башню степеней. При некоторых $x$ такая бесконечная конструкция имеет смысл. Видно, что $f(x)=x^{f(x)}\,\!$. Это значит, что при $f(x)=a\,\!$ имеем $a=x^{a}\,\!$ и $x=a^{1/a}\,\!$ . Из этого можно ошибочно вывести, что при $x=3^{1/3}={\sqrt[{3}]{3}}\,\!$ имеем $f(x)=3\,\!$ . Это не так. Убедитесь в этом сами. Мы предлагаем вам изучить функцию $f(x)\,\!$ , приблизив её «башней» степеней конечной высоты. А именно, рассмотрите и запрограммируйте рекурсивную функцию $F_{n}(x)=x^{F_{n-1}(x)}\,\!$ , $F_{1}(x)=x\,\!$ . Напишите программу (или используйте GnuPlot, Mathematica, Maple или другие подобные инструменты), которая рисует график четырёх функций $F_{1000}(x)\,\!$, $F_{1001}(x)\,\!$, $F_{1002}(x)\,\!$, $F_{1003}(x)\,\!$ на интервале (0;3). Из вида этих графиков сделайте выводы. При каких $x$ определена функция $f(x)\,\!$, то есть при каких $x$ значения $F_{n}(x)\,\!$ при увеличении $n$ становятся все ближе и ближе друг к другу и к некоторому фиксированному числу? Когда числа из последовательности $a_{1},\;a_{2},\;\;\ldots ,\;a_{n},\;\ldots \,\!$ приближаются к некоторому числу $A$, то говорят, что последовательность $a_{i}\,\!$ имеет предел $A$ и пишут $a_{n}\to A\,\!$ при $n\to \infty \,\!$.

Уточним это понятие: какое бы маленькое положительное число $\varepsilon \,\!$ мы ни выбрали, всегда найдется такой элемент последовательности $a_{M}\,\!$, что он сам и все элементы последовательности после него оказываются удалены от числа A не более чем на $\varepsilon \,\!$

$|a_{M}-A|<\varepsilon ,\quad |a_{M+1}-A|<\varepsilon ,\quad |a_{M+2}-A|<\varepsilon ,\quad |a_{M+3}-A|<\varepsilon ,\quad \ldots \,\!$

Последний вопрос нашей задачи можно переформулировать так: при каких значениях $x$ у последовательности $a_{1}=F_{1}(x),\;a_{2}=F_{2}(x),\;a_{3}=F_{3}(x),\;\ldots \,\!$ есть предел, то есть существует некоторое число $A$, зависящее от $x$, к которому стремится последовательность $a_{i}=F_{i}(x)\,\!$?

Примечания

• ^  Это не совсем так. Есть системы символических вычислений, такие как Mathematica и MathCad, которые рассматривают рациональные числа, и числа вроде ${\sqrt {2}}\,\!$, $\pi \,\!$ с любой точностью. В последнем случае это означает, что система может по мере необходимости вычислить любое число знаков после запятой.