Открыть на GitHub

Введение в численное моделирование

Методы вычисления элементарных функций

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

Разложение в ряд Тейлора

Если вы только задумались над тем, как вычисляются значения функций, то вероятно первым пришедшим вам в голову способом будет использовать ряд Тейлора, не зря же его проходили на первом курсе. Например, мы знаем, что экспоненту числа можно найти как сумму ряда: \(e^x=\sum_{n=0}^\infty \frac{x^k}{k!}=:S.\) По определению, если ряд сходится, то последовательность его частичных сумм (сумм первых нескольких слагаемых) сходится к сумме ряда с ростом числа слагаемых. Так как $N$-ая частичная сумма \(S_N=\sum_{n=0}^N \frac{x^k}{k!}=\sum_{n=0}^N a_k x^k\) содержит только умножения и сложения, а коэффициенты $a_k$ можно расчитать заранее или найти по рекуррентной формуле: \(a_0=1, a_k=k\cdot a_{k-1}\text{ для }k>0,\) то вычисления легко провести, имея только реализацию арифметических операций.

Так $N$-ая частичная сумма степенного ряда представляет собой многочлен степени $N$, то найте его значение в одной точке можно за $N$ flops (операций с плавающей запятой, одна операция включает в себя одно сложение и одно умножение) по схеме Горнера: \(S_N=\sum_{n=0}^N a_n x^n=a_0+x(a_1+x(a_2+\ldots x(a_N)\ldots)).\)

При вычислении функций мы обращаем главным образом на две характеристики используемого метода:

Как мы видели, сложность вычислений суммы ряда пропорциональна числу использованных слагаемых. Точность вычислений также будет зависеть от числа слагаемых, так как отбрасывание остатка ряда естественно изменяет результат. Выясним, как найти число слагаемых, которое необходимо взять, чтобы вычислить значение экспоненты с наперед заданной точностью $\epsilon$. Будем пока учитывать только погрешность, возникующую из-за отбрасывания остатка ряда. Тогда точное значение есть сумма $S$ ряда, а приближенное значение $S_N$, откуда относительная погрешность равна: \(\delta_N=\frac{|S-S_N|}{|S|} =\frac{\bigg|\sum_{n=0}^\infty \frac{x^n}{n!}-\sum_{n=0}^N \frac{x^n}{n!}\bigg|}{|\sum_{n=0}^\infty \frac{x^n}{n!}|} =e^{-x}\bigg|\sum_{n=N+1}^\infty \frac{x^n}{n!}\bigg|.\) Мы хотим найти такое $N$, что найденная выше погрешность была меньше $\epsilon$. Заменим сумму большей величиной, которую мы сможем вычислить, в результате мы найдем несколько большее $N$, чем необходимо, однако получим простую формулу для оценки: \(\delta_N\leq e^{-x}\sum_{n=N+1}^\infty \frac{|x|^n}{n!} \leq e^{-x}x^{N+1}\sum_{k=0}^\infty \frac{|x|^k}{(k+N+1)!} \leq e^{-x}\frac{x^{N+1}}{(N+1)!}\sum_{k=0}^\infty \frac{|x|^k}{k!} \leq e^{-x}\frac{x^{N+1}}{(N+1)!}e^{|x|} =\frac{x^{N+1}}{(N+1)!}\begin{cases}1,& x>0\\ e^{2|x|},& x<0 \end{cases}.\) Таким образом относительная погрешность меньше первого отброшенного слагаемого, по крайней мере для положительных значений аргумента. Следовательно суммирование достаточно производить до тех пор, пока первое отброшенное (можно последнее вычисленное) слагаемое больше желаемой относительной погрешности.

Возникает соблазн заключить, что для любого ряда относительная погрешность отбрасывания остатка ряда равна первому отброшенному слагаемому, однако это не так и в каждом конкретном случае нужно начинать анализ с самомого начала, с вычисления относительной погрешности по определению. Для оценки остатка ряда Тейлора можно воспользоваться формулой остаточного члена ряда Тейлора, например, в форме Лагранжа: если $f(x)=\sum_{n=0}^\infty a_n(x-x_0)^n$, то \(R_{N}=\sum_{n=N+1}^\infty a_n x^n =\frac{(x-x_0)^{N+1}}{(N+1)!}f^{(N+1)}(x_0+\theta(x-x_0)),\) для некоторого $\theta\in[0,1]$, где $f^{(n)}$ обозначает $n$-ую производную. Так для экспоненты остаток оценивается следующим образом: \(|R_N|=\bigg|\frac{x^{N+1}}{(N+1)!}e^{\theta x}\bigg| \leq \frac{|x|^{N+1}}{(N+1)!} \begin{cases}e^{|x|}, & x>0,\\ 1, & x<0,\end{cases},\) что даже несколько лучше полученной нами раньше оценки для отрицательных значений аргумента.

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

В действительности, найденная нами погрешность суммирования является черезмерно оптимистичной. Вспомним для начала, что исходные данные также имели погрешность, поэтому даже самый лучший алгоритм вычисления значения функции будет возвращать значение с погрешностью, равной произведению погрешности исходных данных и числа обусловленности для экспоненты. Как мы знаем из прошлой лекции, число обусловленности можно вычислить следующим образом: \(\kappa(f(x))=\bigg|\frac{xf'(x)}{f(x)}\bigg|,\) т.е. для экспоненты число обусловленности равно: \(\kappa(e^x)=\bigg|\frac{xe^x}{e^x}\bigg|=|x|.\) Таким образом, для больших значений $x$ экспоненту нельзя вычислить с машинной точностью даже теоретически, если $x$ имело точность меньше или равную машинной.

Для отрицательных значений аргумента ситуация обстоит еще хуже. Действительно, $n$-ое слагаемое в сумме вычисляется с относительной погрешностью $n\epsilon$, где $\epsilon$ – относительная погрешность для $x$. Тогда абсолютная погрешность каждого слагаемого $n\epsilon |x|^n/n!$, причем погрешности нескольких первых членов ряда могут очень велики для больших отрицательных $x$, и как минимум не меньше $\epsilon |x|$. Так как при суммировании абсолютные погрешности складываются, то вычисление частичной суммы имеет абсолютную погрешность не меньше $\epsilon|x|$. Так как точное значение экспоненты быстро убывает при $x\to-\infty$, то относительная погрешность результата $\epsilon|x|/e^x$ быстро растет. В результате формула суммирования через ряд Тейлора практически не может использоваться для отрицательных $x$.

Ряд Тейлора для экспоненты очень быстро сходися, т.е. слагаемые быстро стремяться к нулю (как $1/n!$), поэтому суммирование этого ряда относительно просто, даже с учетом приведенных выше замечаний. Многие другие ряды сходятся гораздо медленнее, например, члены ряда для логарифма убывают только как $1/n$: \(\ln(1+x)=\sum_{n=1}^N \frac{(-1)^{n-1}}{n}x^n,\) это означает, что для получения той же точности, потребуется значительно большее количество слагаемых, более того, максимально достижимая точность будет много меньше.

В отличии от ряда для экспоненты, ряд для логарифма сходится не при всех значениях $x$, а только для $|x|<1$. Т.е. суммой ряда принципиально нельзя пользоваться для аргументов $x\notin(0,2)$. Более того, остаток ряда быстро растет при стремлении $x$ к границам интервала, при этом количество слагаемых, необходимых для достижения точности стремится к бесконечности, а максимальная достижимая точность падает. Тем самым формула суммирования применима только для $x$ близких к $1$. Подытоживая можно сказать, что формулу суммирования через ряд Тейлора имеет смысл использовать только для быстросходящихся рядов и только для значений аргумента, близких к точке разложения в ряд.

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

Интерполяция многочленами

Пусть даны значения $f_k$ функции $f$ в конечном наборе точек $x_k$, $k=0..K$. Задача интерполяции заключается в восстановлении значений функции в промежуточных точках. Так как существует бесконечное число функций, принимающих данные значения, то восстановить значения однозначно невозможно. Однако если предположить, что функция имеет некий специальный вид, то оказывается возможным однозначно восстановить в промежуточных точках, т.е. построить интерполяционную функцию. Из функционального анализа известно, что любую непрерывную функцию на отрезке можно приблизить многочленом, поэтому в качестве интерполяционной функции естественно взять многочлен (однако это не обязательно лучший выбор). В следующей лекции мы увидим, как построить такой многочлен.

Использование комплексных чисел

Для вычисления синусов и косинусов не обязательно использовать специальная алгоритмы, можно воспользоваться уже известным алгоритмом для экспоненты. Это возможно благодаря формуле Эйлера: \(\exp i\theta=\cos\theta+i\sin\theta,\) согласно которой косинус и синус находятся как вещественная и мнимая часть экспоненты комплексного аргумента.

Сложение двух комлпексных чисел требует только два сложения: \((x+iy)+(a+ib)=(x+a)+i(y+b).\) Умножение комплексных чисел определено так: \((x+iy)(a+ib)=(xa-yb)+i(xb+ya),\) однако его можно вычислить, выполняя только три умножения, воспользовавшись наблюдением Карацубы: \((x+iy)(a+ib)=(xa-yb)+i((x+y)(a+b)-xa-yb),\) естественно, для вычисления мнимой части мы переиспользуем ранее вычисленные значения $xa$, и $yb$.

Параллельные вычисления

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

Литература для дальнейшего изучения