Большинство задач, решаемых вычислительной математикой, имеют одну из двух форм:
Ax = b (b Î Y) (1)
y = Ab (b Î X) , (2)
где A – некоторый абстрактный оператор, действующий из полного нормированного пространства X в пространство Y ( A: X ® Y).
Задачам первого типа соответствуют (в зависимости от оператора A) уравнения различного рода.
Примеры: 1) Если A – вещественная функция вещественной переменной (A: R ® R), то речь идет о решении обычного уравнения (линейного, квадратичного, нелинейного); 2) если A – дифференциальный оператор, то решается дифференциальное уравнение (обыкновенное или с частными производными); 3) если A – матрица, то решается система линейных алгебраических уравнений.
Задачи второго типа соответствуют задачам, аналогичным вычислению значения функции (значение оператора на элементе из пространства X). Например, если A – определенный интеграл, то речь идет о нахождении значения определенного интеграла от некоторой интегрируемой функции (пространство X– множество интегрируемых на заданном интервале функций).
Основным методом решения задач (1), (2) является дискретизация, то есть замена (аппроксимация) элементов исходных пространств X и Y конечными векторами из эвклидовых пространств En и Em, а оператор A заменяется (аппроксимируется) матрицей Am,n размерности m*n. Таким образом, задачи (1), (2) аппроксимируются задачами (1/) и (2/):
Am,n xn = bm (xn Î En, bm Î Em) (1/)
ym = Am,n bn (ym Î Em, bn Î En) (2/)
Как правило, при аппроксимации задач (1), (2) задачами (1/), (2/) качество аппроксимации (ее точность) характеризуется некоторым числовым параметром h, связанным с n и m. Поэтому можно обозначить Am,n = A(h). Для корректно поставленных задач (и исходных и аппроксимирующих) всегда справедливо утверждение о том, что если при h ® 0 оператор A(h) в некотором смысле сходится к оператору A (по некоторой норме), то и приближенные решения xn и ym сходятся к точным решениям (h ® 0 – эквивалентно n ® ¥ и m ® ¥). Иными словами, при приближенном решении задач (1), (2) мы всегда имеем дело с процессом предельного перехода (точное решение – это предел последовательности решений аппроксимационных задач). Тогда можно говорить о скорости сходимости этого процесса и об его ускорении.
Говорят, что процесс предельного перехода имеет скорость сходимости a, если можно установить справедливость неравенства
ôx* – xnô = C∙ha + O(ha+1) или ôy* – ymô = C∙ha + O(ha+1) ,
где x*, y* – решения соответствующих точных задач, xn, ym – решения аппроксимационных задач, C – константа (не зависит от h), ô∙ô – норма в соответствующем пространстве, a > 0, h < 1.
В нашем контексте под ускорением сходимости будем понимать нахождение линейных комбинаций
(h) = a1x(h1) + a2x(h2) + … + apx(hp),
(h) = a1y(h1) + a2y(h2) + … + apy (hp),
причем hi £ h < 1 при i = 1, …p. Пусть каждое приближенное решение x(hi) и y(hi) имеют скорость сходимости a, то есть
ôx* – x(hi)ô= C∙hia + O(hia+1), ôy* – y(hi)ô= C∙hia + O(hia+1).
Если построенные решения будет удовлетворять равенствам
ôx* – (h)ô = C1∙hb + O(hb+1), ôy* – (h)ô= C∙hb + O(hb+1) при b > a,
то это решение дает лучшую скорость сходимости по сравнению с решениями x(hi), из которых оно составлено (произошло улучшение сходимости).
Рассмотрим три примера ускорения сходимости в различных методах приближенных вычислений.
1. Метод Эйлера численного решения задачи Коши для обыкновенного дифференциального уравнения первого порядка.
Пусть требуется найти решение задачи Коши
x/ (t) = f(t, x), t Î [0; T], при x(0) = g. (3)
Метод Эйлера заключается в следующем. Отрезок [0; T] разбивается равноотстоящими точками (узлами) tk = k∙h, k = 0, 1, …, n, h = T/n. Это –дискретизация промежутка [0; T], на котором ищется решение, – этот отрезок заменили дискретным множеством точек {tk}. Затем в каждой узловой точке производная x/(tk) заменяется (аппроксимируется) разностным отношением
x/(tk) » .
Тогда исходная задача (3) аппроксимируется дискретной задачей
xk+1 = xk + h∙f(tk, xk), x0 = g, при k = 0, 1, …, n–1. (3/)
Алгоритмически получилась очень простая схема: по известному значению x0 = g вычисляется x1 и т.д. Схема нахождения приближенного решения задачи Коши в узловых точках относится к рекуррентным формулам. В (3/) мы воспользовались обозначением xk » x*(tk) (xk – приближенное решение, x*(tk) – точное решение в узловой точке).
В качестве пространства X в задаче (3) можно взять пространство непрерывно дифференцируемых на [0; T] функций, Y – множество непрерывных на [0; T] функций, A – дифференциальный оператор (Ax ≡ x/(t)). В качестве нормы и в пространстве X, и в пространстве Y возьмем
ôxô = maxïx(t)ï, при t Î [0; T].
Пространство Xn – множество векторов x = (x0, x1, …, xn), а пространство Ym – множество векторов y = (y0, y1, …, yn–1). В качестве нормы в этих пространствах возьмем первую норму вектора, согласованную с нормами в пространствах X и Y: ôxôI = maxïxkï, по всем k = 1, …, n. Вид дискретного оператора Am,n определяется тождеством Am,n x ≡ , при k = 0, 1, …, n–1.
Выясним, какова погрешность метода Эйлера на одном шаге. Будем считать, что на k-м шаге известно точное решение, то есть xk = x*(tk). Тогда абсолютная погрешность e(tk+1) на k + 1 шаге определяется разностью
e(tk+1) = x*(tk+1) – xk+1.
В предположении существования у точного решения x*(t) непрерывной производной третьего порядка можно x*(tk+1) представить в виде отрезка ряда Тейлора с остаточным членом
x*(tk+1) = xk + h∙x*/(tk) + ∙ + ∙ , q Î [tk; tk+1].
Отсюда x*(tk+1) = xk + h∙f(tk, xk) + ∙ + ∙ =
= xk+1 + ∙ + ∙ .
Обозначим ï ï = C. Тогда имеем
ïe(tk+1)ï = ïx*(tk+1) – xk+1ï = С∙h2 + O(h3).
Показали, что скорость сходимости метода Эйлера на одном шаге имеет порядок 2 (a = 2).
Выберем h1=h/2, h2 = h. Опираясь на значение в точке tk = k∙h, найдем три приближенных решения:
xk+1 = xk + h∙f(tk, xk),
= xk + ∙f(tk, xk), = + ∙f(tk + , ).
Таким образом, в точке tk+1 имеем два приближенных решения: первое – вычислено при шаге разбиения h, второе – . По доказанному ранее свойству имеем
x*(tk+1) = xk+1 + C∙h2 + O(h3). (4)
Для имеем
x*(tk + ) = + x//(h1) = + x//(tk) + O(h3),
В точке tk+1 имеем
x*(tk + h) = + x//(h1) + f(tk + , ) + x//(h2) =
= + x//(h1)+ x//(h2) – × x//(h1)=
= + x//(h)– × x//(h1)= + x//(tk) + O(h3),
то есть
x*(tk + h) = + x//(tk) + O(h3), или
x*(tk+1) = + 2C∙( )2 + O(h3). (5)
Умножим равенство (4) на a1, а (5) – на a2 и сложим получающиеся равенства
a1×x*(tk+1) + a2 ×x*(tk+1) = a1×xk+1 + a2× + (a1×h2 + a2× )×С + O(h3).
Пренебрегая слагаемыми порядка h3, выбирая a1 и a2 так, чтобы выполнялись два равенства
a1 + a2 = 1,
a1 + = 0 ,
получаем a1 = –1, a2 =2. С помощью линейной комбинации –xk+1 + 2∙ , получаем новое решение, скорость сходимости которого на единицу лучше, чем у решений полученных при значениях шага h и . Конкретнее, имеем
x*(tk+1) – [–xk+1 + 2∙ ] = O(h3).
Итак, вдвое увеличивая количество вычислений, улучшаем точность (повышаем скорость сходимости) на порядок. Увеличение точности на порядок без приема ускорения скорости сходимости потребовало бы увеличение количества вычислений приблизительно в 1/h раз.
В предположении большей гладкости решения x*(t) (а, значит, и функции f(t, x)) аналогичным образом можно повысить скорость сходимости еще на единицу и т.д.
Пример. Рассмотрим задачу Коши для обыкновенного дифференциального уравнения
x/ (t)= –x2 + t× x + ; x(0) = 1. (6)
T
0.1
0.05
x**
0.025
0.0125
x**
x*(t)
0.1
0.95491
0.95710
0.95931
0.95836
0.95877
0.95918
0.95917
0.2
0.92704
0.93059
0.93414
0.93263
0.93329
0.93395
0.93395
0.3
0.91265
0.91704
0.92144
0.91958
0.92041
0.92124
0.92123
0.4
0.90924
0.91417
0.91911
0.91703
0.91797
0.91891
0.91890
0.5
0.91508
0.92036
0.92566
0.92345
0.92446
0.92547
0.92546
0.6
0.92900
0.93452
0.94004
0.93775
0.93881
0.93987
0.93986
0.7
0.95017
0.95585
0.96154
0.95918
0.96028
0.96138
0.96137
0.8
0.97801
0.98381
0.98962
0.98722
0.98835
0.98947
0.98946
0.9
1.01213
1.01803
1.02392
1.02149
1.02264
1.02378
1.02377
1.0
1.05225
1.05822
1.06419
1.06173
1.06289
1.06405
1.06405
В таблице приведены результаты решения задачи (6) методом Эйлера при различных значениях шага интегрирования h = 0.1; 0.05; 0.025; 0.0125, причем решения распечатаны не все, а только в узлах кратных 0.1. Через x** обозначено решение, полученное методом ускорения сходимости. В четвертом столбце это решение было получено линейной комбинацией решений в методе Эйлера для h = 0.1 и h = 0.05, а в седьмом – h = 0.025 и h = 0.0125. В последнем столбце приведены точные решения. Выводы может сделать сам читатель.
2. Приближенное вычисление определенного интеграла.
Пусть стоит задача вычисления определенного интеграла
y = .
В качестве метода выберем квадратурную формулу трапеций (составную).
y(h) = yn = , (7)
где xk = k×h, k = 0, 1, …, n, h = .
Найдем представление остаточного члена квадратурной формулы трапеций на отрезке [xk; xk+1], то есть
R1(f; xk) = – ×[f(xk) + f(xk+1)].
В предположении существования третьей производной у подинтегральной функции f(x), разлагая ее в ряд Тейлора с остаточным членом в точке xk, используя обобщенную теорему о среднем значении интеграла, можно получить следующее представление:
R1(f; xk) = – ∙f //(xk) + ∙[f ///(hk) – f ///(xk)],
где hk, xk Î [xk; xk+1]. Отсюда ясно, что скорость сходимости формулы трапеций на отрезке [xk; xk+1] имеет порядок О(h3).
Составная формула трапеций (7) получается, как сумма простых формул трапеций при интегрировании на частичных отрезках [xk; xk+1], поэтому остаточный член интегрирования по составной формуле трапеций является суммой остаточных членов на [xk; xk+1], то есть
Rтр(f, h) = –yn = = – f //(xk) + [f ///(hk) – 2 f ///(xk)].
Уменьшая шаг в два раза, то есть h1 = , получим новую составную формулу трапеций с остаточным членом
Rтр(f, h1) = –y(h1) = – f( ) – f( ) +
+ [f ///( ) – 2 f ///( )].
В последнем представлении = k∙ , при k = 0, 1, …, N, N = 2n. Таким образом, первые суммы (без коэффициентов) в представлениях Rтр(f, h) и Rтр(f, h1) совпадают. Обозначим эту сумму через М. Вторую сумму заменим на следующую
f( ) = f(xk) + f /(γk) = M + f /(γk),
где γk Î [xk; xk+1].
Легко заметить, что
n∙ min f(x) ≤ f(xk) ≤ n∙ max f(x),
поэтому существует такая точка η Î [0; 1] , что f(xk)= n∙f (η)= ∙ f (η). Рассуждая аналогично, можно получить
[f ///(hk) – 2 f ///(xk)] = [f ///(α) – 2f ///(β)],
[f ///( ) – 2 f ///( )] = [f ///(α1) – 2f ///(β1)],
f /(γk) = f /(γ),
где α, β, α1, β1, γ – некоторые числа из [ 0; 1 ]. Тогда можно записать
Rтр(f, h) = – ∙ f (η ) + O(h3), Rтр(f, h1) = – ∙ f (η ) + O(h3).
Итак, имеем = y(h) – ∙ f (η ) + O(h3),
= y( ) – ∙ f (η ) + O(h3).
Умножая первое равенство на C1, а второе – на С2 и складывая полученные равенства, приходим к соотношению
(С1 + С2)∙ = С1 ∙y(h) + C2 ∙y( ) – ∙[C1 + C2]∙ f (η ) + O(h3).
Полагая C1 = – , а C2 = , приходим к равенству
= – ∙ y(h) + ∙ y( ) + O( h3).
Получили новую формулу численного интегрирования, которая на порядок точнее двух исходных
» – ∙ y(h) + ∙ y( ). (8)
Таким образом, рассмотренный прием улучшения сходимости квадратурной формулы трапеций можно рассматривать как один из способов получения новых неинтерполяционных квадратурных формул большей точности.
Пример. Вычислить интеграл
I = .
Решения сведены в таблицу
n=20
n=40
n=80
20 – 40
40 – 80
I (точное)
0.6353102
0.6362925
0.6365380
0.6366199
0.6366198
0.6366198
В трех первых столбцах приведены приближенные значения интеграла I, полученные по составной формуле трапеций при различных значениях n, в четвертом и пятом – по формуле (8), в последнем – точное значение интеграла I.
3. Предельный переход и машинная арифметика.
В качестве модельной задачи рассмотрим задачу нахождения точных знаков числа e, опираясь непосредственно на определение числа e как предела
(1 + )n e.
Построим последовательность чисел {en} таких, что
en = (1 + )n . (9)
Алгоритмически процесс построения этой последовательности прост. Если бы вычисления по формуле (9) осуществлялись точно (без округлений), то в соответствии с определением предела, построив достаточно длинную последовательность {en}, можно было бы получить сколь угодно точное представление числа e (с любым, заданным, количеством верных знаков). Однако, рассматривая фактически построенную (вычисленную на компьютере) последовательность {en}, замечаем, что характер изменения членов последовательности не совпадает с теоретически ожидаемым изменением. При увеличении n до некоторого достаточно большого N действительно происходит стабилизация верных знаков, то есть последовательность ведет себя в соответствии с теорией. Дальнейшее увеличение n не добавляет верных знаков. Наконец, наступает момент, когда верные знаки начинают замещаться сомнительными знаками, вплоть до полного искажения результата (ни одного верного знака). Эта ситуация характерна при осуществлении предельного перехода при реализации его на вычислительных устройствах.
Рассмотрим структуру абсолютной погрешности a = e – en. Она состоит из двух слагаемых, первое из которых определяется заменой предела последовательности e некоторым ее членом en, а второе – погрешностью округления и количеством цифр в представлении чисел данного вычислительного устройства (разрядностью чисел).
Рассмотрим функцию f(n) = (1 + )n. Вводя новую переменную x = и доопределяя ее по непрерывности в точке x = 0, получаем хорошо изученную функцию y = f(x), которая непрерывна в точке x = 0 вместе со всеми производными. Легко вычислить, что f /(0) = . Таким образом, первое слагаемое в представлении абсолютной погрешности может быть записано в виде ×x ( × ). Учитывая, что мантиссы числа типа real в Pascal записываются в пяти байтах, причем приближенно (с округлением), погрешность числа равна 2–39, погрешность числа (1 + )n равна 2–38. Тогда вычислительная погрешность степени (1+ )n будет равна n×(1 + )n–1×2–38 » n×e×2–38. При больших значениях n эта погрешность будет влиять не только на знаки, в которых осуществляется округление, но и на информативные знаки. Так, например, при n » 238 величина погрешности будет около e (2,7…), так что у соответствующего значения en все знаки сомнительные. В связи с этим возникает проблема, каким выбирать значение n и как видоизменить вычислительный процесс, чтобы получить максимальное количество верных знаков, (например, 11– стандартная выдача числа типа real на Pascal). Ясно, что при очень большом n вычислительная погрешность не позволит найти даже меньше, чем одиннадцать верных знаков. С другой стороны, при небольших значениях n в принципе невозможно получить много верных знаков.
В создавшейся ситуации выходом из положения является применение приема ускорения сходимости. Вычислим при сравнительно небольших значениях n несколько (n = n1, n2,, …, np) чисел en. Затем с помощью линейной комбинации этих значений найдем более точное приближение. При этом, если коэффициенты линейной комбинации выбраны так, что точность ее выше точности каждого из составляющих, то получим большее количество верных знаков. Выбор небольших значений чисел n сделают вычислительную погрешность сравнительно малой.
Найдем коэффициенты упомянутых линейных комбинаций, установим точность, с которой линейная комбинация приближает число e, и проведем численный эксперимент.
Ранее введенная функция f(x) может быть разложена в ряд Тейлора в точке x = 0 с остаточным членом
f(x) = f(0) + f /(0)×x + + …+ + .
Возвращаясь к переменной n, учитывая, что f(0) = e, а f( )= en, перепишем это равенство в следующем виде
e =en – f /(0)× – … – ∙ – ∙ .
Составим линейную комбинацию E = C1×E1 + C2×E2 + … + Cp×Ep , выбрав в качестве Ci решения системы уравнений
C1 + C2 + … + Cp = 1,
C1× h1 + C2×h2 + … + Cp×hp = 0,
........................................................,
C1× h1p + C2×h2p + … + Cp×hp p = 0,
где через Ei обозначены en при соответствующих значениях ni, а hi = , при
i = 1, …, p. Написанная система линейных уравнений разрешима единственным образом, так как ее определитель – это транспонированный определитель Вандермонда. Удобнее всего величины hi выбрать равными при некотором фиксированном h. Приведем системы линейных алгебраических уравнений для определения коэффициентов Ci и их решения при различных значениях параметра p.
1. p = 1: C1 = 1. Точность – порядка O(h) = O( ).
2. p = 2: C1 + C2 = 1, C1 = –1 , C2 = 2.
C1 + C2 = 0.
3. p=3: C1 + C2 + C3 = 1, C1 = , C2 = –2, C3 = .
C1 + C2 + C3 =0,
C1 + C2 + C3 = 0.
4. p = 4: C1 + C2 + C3 + C4 = 1, C1 = – 0,0476190476,
C1 + C2 + C3 + C4 = 0, C2 = 0,6666666667,
C1 + C2 + C3 + C4 = 0, C3 = – 0,2666666667
C1 + C2 + C3 + C4 = 0. C4 = 3,0476190476.
5. p = 5: C1 + C2 + C3 + C4 + C5 = 1,
C1 + C2 + C3 + C4 + C5 = 0,
C1 + C2 + C3 + C4 + C5 = 0,
C1 + C2 + C3 + C4 + C5 = 0,
C1 + C2 + C3 + C4 + C5 = 0,
С1 = 0,00317460316, С2 = – 0,09523809522,
С3 = 0,88888888889, С4 = – 3,0476190476.
Находить линейные комбинации для p>5 нет смысла, потому что при решении систем линейных алгебраических уравнений высокого порядка с коэффициентами отмеченного характера, например, методом Гаусса вычислительная погрешность чисел Ci становится больше чем 10–10 и улучшения результата не происходит.
Ниже приведена таблица с результатами численного эксперимента.
n
P = 1
p = 2
p = 3
p = 4
p = 5
16
2,6379284974
2,7160517614
2,7182491138
2,7182815824
2,712818276
32
2,6769901294
2,7176997757
2,7182775238
2,7182818123
2,7182818283
64
2,6973449526
2,7181330868
2,7182812763
2,7182818273
2,7182818284
128
2,7077390197
2,7182442289
2,7182817584
2,7182818283
2,7182818284
256
2,7129916243
2,7182723761
2,7182818196
2,7182818284
512
2,7156320002
2,7182794587
2,7182818273
1024
2,7169557294
2,7182812351
2,7182818285
2048
2,7176184823
2,7182816801
4096
2,7179500813
2,7182817911
8192
2,7181159362
Автор не может не удержаться от комментариев по поводу поведения чисел, приведенных в этой таблице. Во-первых, числа en в первом столбце, начиная с n=32 и оканчивая n=512, уже содержат всю информацию, которая необходима для решения поставленной задачи (нахождение числа e с точностью одиннадцать верных знаков). При этом каждое из них дает приближение к числу e лишь с точностью два верных знака. Во-вторых, используя прием ускорения сходимости с помощью линейной комбинации пяти en при n = 32, 64, 128, 256, 512, мы производим лишь 992 операции возведения в степень, тогда как нахождение приближенного значения числа e без ускорения сходимости дает максимально возможную точность девять верных знаков только на предельных значениях целых чисел типа longint (количество операций возведения в степень – 231).
