резок пополам, строим последовательность приближений и складываем частичные суммы с помощью функ-
ции zipWith.
Эта версия функции хоть и наглядная, но не эффективная. Функция f вычисляется заново при каждом ре-
курсивном вызове. Было бы хорошо вычислять её только для новых значений. Для этого мы будем передавать
значения с предыдущего шага:
integrate :: Fractional a => (a -> a) -> a -> a -> [a]
integrate f a b = integ f a b (f a) (f b)
where integ f a b fa fb = (fa+fb)*(b—a)/2 :
zipWith (+) (integ f a m fa fm)
(integ f m b fm fb)
where m
= (a + b)/2
fm = f m
182 | Глава 11: Ленивые чудеса
В этой версии мы вычисляем значения в функции f лишь один раз для каждой точки. Запишем итоговое
решение:
int :: (Ord a, Fractional a) => a -> (a -> a) -> a -> a -> a
int eps f a b = converge eps $ integrate f a b
Мы опять воспользовались функцией converge, нам не нужно было её переписывать. Проверим решение.
Для проверки также воспользуемся экспонентой. В прошлой главе мы узнали, что
? x
ex = 1 +
etdt
0
Посмотрим, так ли это для нашего алгоритма:
*Numeric> let exp’ = int 1e-5 exp 0
*Numeric> let test x = abs $ exp x — 1 —
exp’ x
*Numeric> test 2
8.124102876649886e-6
*Numeric> test 5
4.576306736225888e-6
*Numeric> test 10
1.0683757864171639e-5
Алгоритм работает. В статье ещё рассмотрены методы повышения точности этих алгоритмов. Что инте-
ресно для улучшения точности не надо менять существующий код. Функция принимает последовательность
промежуточных решений и преобразует её.
11.2 Степенные ряды
Напишем модуль для вычисления степенных рядов. Этот пример взят из статьи Дугласа МакИлроя
(Douglas McIlroy) “Power Series, Power Serious”. Степенной ряд представляет собой функцию, которая опре-
деляется списком коэффициентов:
F ( x) = f 0 + f 1 x + f 2 x 2 + f 3 x 3 + f 4 x 4 + …
Степенной ряд содержит бесконечное число слагаемых. Для вычисления нам потребуются функции сло-
жения и умножения. Ряд F ( x) можно записать и по-другому:
F ( x) = F 0( x)
= f 0 + xF 1( x)
= f 0 + x( f 1 + xF 2( x))
Это определение очень похоже на определение списка. Ряд есть коэффициент f 0 и другой ряд F 1( x)
умноженный на x. Поэтому для представления рядов мы выберем конструкцию похожую на список:
data Ps a = a :+: Ps a
deriving (Show, Eq)
Но в нашем случае списки бесконечны, поэтому у нас лишь один конструктор. Далее мы будем писать
просто f + xF 1, без скобок для аргумента.
Определим вспомогательные функции для создания рядов:
p0 :: Num a => a -> Ps a
p0 x = x :+: p0 0
ps :: Num a => [a] -> Ps a
ps []
= p0 0
ps (a:as) = a :+: ps as
Обратите внимание на то, как мы дописываем бесконечный хвост нулей в конец ряда. Теперь давайте
определим функцию вычисления ряда. Мы будем вычислять лишь конечное число степеней.
eval :: Num a => Int -> Ps a -> a -> a
eval 0 _
_ = 0
eval n (a :+: p) x = a + x * eval (n—1) p x
В первом случае мы хотим вычислить ноль степеней ряда, поэтому мы возвращаем ноль, а во втором
случае значение ряда a+ xP складывается из числа a и значения ряда P умноженного на заданное значение.
Степенные ряды | 183
Арифметика рядов
В результате сложения и умножения рядов также получается ряд. Также мы можем создать ряд из числа.
Эти операции говорят о том, что мы можем сделать степенной ряд экземпляром класса Num.
Сложение
Рекурсивное представление ряда f + xF позволяет нам очень кратко выражать операции, которые мы
хотим определить. Теперь у нас нет бесконечного набора коэффициентов, у нас всего лишь одно число и ещё
один ряд. Операции существенно упрощаются. Так сложение двух бесконечных рядов имеет вид:
F + G = ( f + xF 1) + ( g + xG 1) = ( f + g) + x( F 1 + G 1)
Переведём на Haskell:
(f :+: fs) + (g :+: gs) = (f + g) :+: (fs + gs)
Умножение
Умножим два ряда:
F ? G = ( f + xF 1) ? ( g + xG 1) = f g + x( f G 1 + F 1 ? G)
Переведём:
(.*) :: Num a => a -> Ps a -> Ps a
k .* (f :+: fs) = (k * f) :+: (k .* fs)
(f :+: fs) * (g :+: gs) = (f * g) :+: (f .* gs + fs * (g :+: gs))
Дополнительная операция (.*) выполняет умножение всех коэффициентов ряда на число.
Класс Num
Соберём определения для методов класса Num вместе:
instance Num a => Num (Ps a) where
(f :+: fs) + (g :+: gs) = (f + g) :+: (fs + gs)
(f :+: fs) * (g :+: gs) = (f * g) :+: (f .* gs + fs * (g :+: gs))
negate (f :+: fs) = negate f :+: negate fs
fromInteger n = p0 (fromInteger n)
(.*) :: Num a => a -> Ps a -> Ps a
k .* (f :+: fs) = (k * f) :+: (k .* fs)
Методы abs и signum не определены для рядов. Обратите внимание на то, как рекурсивное определение
рядов приводит к рекурсивным определениям функций для рядов. Этот приём очень характерен для Haskell.
Поскольку наш ряд это число и ещё один ряд за счёт рекурсии мы можем воспользоваться операцией, которую
мы определяем, на “хвостовом” ряде.
Деление
Результат деления Q удовлетворяет соотношению:
F = Q ? G
Переписав F , G и Q в нашем представлении, получим
f + xF 1 = ( q + xQ 1) ? G = qG + xQ 1 ? G = q( g + xG 1) + xQ 1 ? G
= qg + x( qG 1 + Q 1 ? G)
Следовательно
q
= f / g
Q 1 = ( F 1 ? qG 1)/ G
Если g = 0 деление имеет смысл только в том случае, если и f = 0. Переведём на Haskell:
184 | Глава 11: Ленивые чудеса
class Fractional a => Fractional (Ps a) where
(0 :+: fs) / (0 :+: gs) = fs / gs
(f :+: fs) / (g :+: gs) = q :+: ((fs — q .* gs)/(g :+: gs))
where q = f/g
fromRational x = p0 (fromRational x)
Производная и интеграл
Производная одного члена ряда вычисляется так:
d xn = nxn? 1
dx
Из этого выражения по свойствам производной
d
d
d
( f ( x) + g( x)) =
f ( x) +
g( x)
dx
dx
dx
d ( k ? f( x)) = k ? d f( x)
dx
dx
мы можем получить формулу для всего ряда:
d F( x) = f 1 + 2 f 2 x + 3 f 3 x 2 + 4 f 4 x 3 + …
dx
Для реализации нам понадобится вспомогательная функция, которая будет обновлять значение допол-
нительного множителя n в выражении nxn? 1:
diff :: Num a => Ps a -> Ps a
diff (f :+: fs) = diff’ 1 fs
where diff’ n (g :+: gs) = (n * g) :+: (diff’ (n+1) gs)
Также мы можем вычислить и интеграл степенного ряда:
int :: Fractional a => Ps a -> Ps a
int (f :+: fs) = 0 :+: (int’ 1 fs)
where int’ n (g :+: gs) = (g / n) :+: (int’ (n+1) gs)
Элементарные функции
Мы можем выразить элементарные функции через операции взятия производной и интегрирования. К
примеру уравнение для ex выглядит так:
dy = y
dx
Проинтегрируем с начальным условием y(0) = 1:
? x
y( x) = 1 +
y( t) dt
0
Теперь переведём на Haskell:
expx = 1 + int expx
Кажется невероятным, но это и есть определение экспоненты. Так же мы можем определить и функции
для синуса и косинуса:
d sin x = cos x,
sin(0) = 0 ,
dx
d cos x = ? sin x, cos(0) = 1
dx
Что приводит нас к:
sinx = int cosx
cosx = 1 — int sinx
И это работает! Вычисление этих функций возможно за счёт того, что вне зависимости от аргумента
функция int вернёт ряд, у которого первый коэффициент равен нулю. Это значение подхватывается и ис-
пользуется на следующем шаге рекурсивных вычислений.
Через синус и косинус мы можем определить тангенс:
tanx = sinx / cosx
Степенные ряды | 185
11.3 Водосборы
В этом примере мы рассмотрим одну интересную технику рекурсивных вычислений, которая называется
мемоизацией (memoization). Она заключается в том, что мы запоминаем все значения, с которыми вызывалась
