Читаем Учебник по Haskell полностью

сходится, то производная определена. Для того чтобы решить эту задачу мы начнём с небольшого значе-

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

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

Этот процесс напоминает то, что мы делали при поиске корня уравнения методом неподвижной точки.

Мы можем взять из того решения функцию определения сходимости последовательности:

| 181

converge :: (Ord a, Num a) => a -> [a] -> a

converge eps (a:b:xs)

| abs (a - b) <= eps

= a

| otherwise

= converge eps (b:xs)

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

вычисляет промежуточные решения:

easydiff :: Fractional a => (a -> a) -> a -> a -> a

easydiff f x h = (f (x + h) - f x) / h

Мы возьмём начальное значение шага и будем последовательно уменьшать его вдвое:

halves = iterate (/2)

Соберём все части вместе:

diff :: (Ord a, Fractional a) => a -> a -> (a -> a) -> a -> a

diff h0 eps f x = converge eps $ map (easydiff f x) $ iterate (/2) h0

where easydiff f x h = (f (x + h) - f x) / h

Сохраним эти определения в отдельном модуле и найдём производную какой-нибудь функции. Проте-

стируем решение на экспоненте. Известно, что производная экспоненты равна самой себе:

*Numeric> let exp’ = diff 1 1e-5 exp

*Numeric> let test x = abs $ exp x - exp’ x

*Numeric> test 2

1.4093421286887065e-5

*Numeric> test 5

1.767240203776055e-5

Интегрирование

Теперь давайте поинтегрируем функции одного аргумента. Интеграл это площадь кривой под графиком

функции. Если бы кривая была прямой, то мы могли бы вычислить интеграл по формуле трапеций:

easyintegrate :: Fractional a => (a -> a) -> a -> a -> a

easyintegrate f a b = (f a + f b) * (b - a) / 2

Но мы хотим интегрировать не только прямые линии. Мы представим, что функция является ломаной

прямой линией. Мы посчитаем интеграл на каждом из участков и сложим ответы. При этом чем ближе точки

друг к другу, тем точнее можно представить функцию в виде ломаной прямой линии, тем точнее будет

значение интеграла.

Проблема в том, что мы не знаем заранее насколько близки должны быть точки друг к другу. Это зависит

от функции, которую мы хотим проинтегрировать. Но мы можем построить последовательность решений.

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

расти вдвое. Как только решение перестанет меняться мы вернём ответ.

Построим последовательность решений:

integrate :: Fractional a => (a -> a) -> a -> a -> [a]

integrate f a b = easyintegrate f a b :

zipWith (+) (integrate a mid) (integrate mid b)

where mid = (a + b)/2

Первое решение является площадью под прямой, которая соединяет концы отрезка. Потом мы делим от-

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

ции 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)