Фельдман Л.П. Сходимость и оценка погрешности параллельных одношаговых блочных методов моделирования динамических систем с сосредоточенными параметрами

Автор: Л.П. Фельдман, д-р техн. наук

Донецкий национальный технический университет

Украина, 83001, Донецк, ул. Артема, 66.

тел. (062)3010856, e-mail: feldman@r5.dgtu.donetsk.ua

Сходимость и оценка погрешности параллельных одношаговых блочных методов моделирования динамических систем с сосредоточенными параметрами

Введение

Доклад содержит обобщение результатов исследований [1,2,3,5], посвященных параллельным методам численного решения задачи Коши для системы обыкновенных дифференциальных уравнений и является продолжением ранее опубликованных работ [6,7,8,9,10]. В нем приводится доказательство сходимости приближенного решения для m- шаговых k точечных блочных методов, что представляет обобщение ранее опубликованных результатов. В [7] приведено доказательство сходимости для одношаговых k точечных, а в [9] –для многошаговых многоточечных методов при k=m Рассмотрены также и методы оценки локальной погрешности решения блочными разностными схемами и даны практические рекомендации их использования для более широкого набора параллельных разностных схем.

1. Параллельные вычислительные схемы многошаговых блочных методов

Для упрощения изложения рассмотрим вначале решение задачи Коши для одного обыкновенного дифференциального уравнения первого порядка.

x" = f(t,x), x(t₀) = x₀ (1)

Рис. 1. Схема разбиения на блоки для одношагового k-точечного метода

блочными методами. Множество M точек равномерной сетки {t_m}, и t_M = T с шагом τ разобьем на N блоков, содержащих k точек каждый, при этом kN£M. В каждом блоке введем номер точки и обозначим через t_n,i точку n блока с номером i. Точку назовём началом блока n, а – концом блока. Очевидно, что имеет место t_n,k = t_n+1,0. Условимся начальную точку t₀ в блок не включать. При численном решении задачи Коши для каждого следующего блока новые k значений функции вычисляются одновременно. Поэтому блочные методы особенно удачно реализуются на параллельных вычислительных системах [1].

Различают два типа блочных методов: одношаговый блочный метод m=1 и многошаговый блочный метод 1 < m. В первом методе только значение в последней точке предшествующего блока используется для вычислений приближенных значений в следующем блоке, тогда как во втором методе используются все или несколько точек предшествующего блока.

Пусть u_n,0 – приближенное значение решения задачи Коши (1) в точке t_n,0 – начальной точке обрабатываемого блока. Предполагается, что в пределах одного блока точки сетки находятся на равных расстояниях: t_n,i= t_n + iτ, i = 1,…,k В общем случае уравнения многошаговых разностных методов для блока из k точек при использовании вычисленных значений приближенного решения в m предшествующих блоку узлах, с учетом введенных выше обозначений можно записать в виде:

, , , (2)

где .

Формулы (2) определяют m-шаговый k-точечный разностный метод.

Для одношаговых разностных методов разностные уравнения имеют вид:

, , (3)

Формулы (3) определяют одношаговый k-точечный разностный метод. Определить коэффициенты a_i,j и b_i,j формул (2) и (3) можно интегро-интерполяционным методом. Построим интерполяционный многочлен L_m+k–1(t) с узлами интерполяции i = –m – 1, m – 2,…,0, 1, k и соответствующим им значениям правой части уравнения (1) F_n,i = f(t_n,i, u_n,i). Проинтегрировав его в пределах (t_n, t_n + iτ), ,

получим уравнения (2) для выбранных m и k Например, формулы для двухточечного двухшагового метода имеют вид

u_n,1= u_n– τ( F_-1+n,1– 13F_n,0 – 13F_n,1 + F_n,2)

u_n,2= u_n– τ (F_n,0 + 4F_n,1 + F_n,2)

2. Погрешность аппроксимации блочных методов.

Выражения для невязoк m-шагового k-точечного разностного метода (2) на решении x(t) исходного дифференциального уравнения имеют вид

(4)

где .

Разлагая и в ряды Тейлора в окрестности точки t_n, получим

Подставляя эти разложения в выражение (4) для невязки, будем иметь

(5)

Сгруппируем члены с одинаковыми производными и изменим порядок суммирования в последнем выражении, тогда получим

(6)

Отсюда следует, что погрешность аппроксимации имеет порядок p, если выполнены условия

(7)

Система уравнений (7) для каждого фиксированного i содержит p уравнений и k+m неизвестных a_i,j, и b_i,j, i, . Положим, чтобы p = k + m, тогда из системы (7) при фиксированном i можно будет определить неизвестные коэффициенты b_i,j, и a_i,j, . Поступая аналогично для каждого , определим все неизвестные коэффициенты b_i,j, и a_i,j, i, , которые в дальнейшем будем считать элементами матриц А и В соответственно. Отсюда следует, что наивысший порядок аппроксимации m-шагового k-точечного блочного метода равен m + k. Его погрешность в соответствии с (6) определяется формулой

(8)

Элементы матриц Bи A можно найти, решая систему (7) для любых заданных k и m.

3. Сходимость и оценка погрешности многошаговых блочных методов.

Обозначим через: u_n – вектор значений приближенного решения в точках n-го блока с компонентами u_n,i, i = 1, 2, …, k; u"_n – вектор, компонента u"_n,i = f(t_n,i, u_n,i), которогоравна значению правой части уравнения (1) в точке n-го блока, соответствующей значению приближенного решения; v"_n-1 – вектор, компонента v"_n-1,j = f(t_{n-1, j}, u_n-1,j), которого равна значению правой части уравнения (1) в точке (n-1)-го блока (если k ≥ m), соответствующей значению приближенного решения.

Запишем систему (2) в векторной форме

D^-1(u_n– u_n,0e) / τ = Bv"_n–1+ Au"_n, (9)

где D = (d_ii) – диагональная матрица с элементами d_ii = i, i = 1, 2, …, k, B – матрица с элементами b_i,j, , , A – матрица с элементами a_i,j, i, , e – единичный вектор-столбец

Обозначим через x_n – вектор значений точного решения задачи (1) в точках t_n блока n. Получим уравнение, которому удовлетворяет вектор z_n= u_n– x_n погрешностей в блоке. Подставим в левую часть уравнения (9) u_n= z_n + x_n добавим к правой части и вычтем из нее выражение By"_n–1 + Ax"_n, где y"_n–1,i = f(t_n–1, x_n–1, i) . Тогда уравнение для погрешности примет вид

D^-1(z_n– z_n,0e) / τ = – D^-1(x_n– x_n,0e) / τ +By"_n–1+ Ax"_n +B(v"_n-1– y"_n–1) +A(u"_n– x"_n).

Входящее в правую часть выражение

r_n=–D^-1(x_n– x_ne) / τ + By"_n–1+ Ax"_n(10)

представляет собой вектор невязок разностных уравнений (9) на точном решении уравнения (1). Поскольку разностные уравнения (9) аппроксимируют исходное уравнение (1) в точках блока с порядком O(τ^k⁺^m), то имеет место оценка

||r_n|| = O(τ^k⁺^m). (11)

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

ω_n= B(v"_n–1– y"_n–1) +A(u"_n– x"_n). (12)

Тогда уравнение для погрешности запишется короче

z_n= z_n,0e + τD(r_n – ω_n). (13)

Вектор-функция ω_n зависит нелинейно от погрешности z_n. Вид этой зависимости определяется функцией f(t, x). В дальнейшем будем предполагать, что f(t, x) удовлетворяет условию Липшица по второму аргументу, т.е.

|f(t, x₁) – f(t, x₂)| ≤ L|x₁ – x₂| (14)

для всех t, x₁, x₂ из рассматриваемой области. Согласно формуле конечных приращений Лагранжа для соответствующих компонент имеем

u"_n,i– x"_n,i= f(t_n,i, u_n,i) – f(t_n,i, x_n,i) = l_n,iz_n,i, i =

v"_n–1,i– y"_n–1,i= f(t_n–1,i, u_n–1,i) – f(t_n–1,i, x_n–1,i) = l_n–1,iz_n–1,i, i =

где

l_n,i= f_x(t_n,i, x_n,i+ θz_n,i), 0 < θ <1, i = ,

l_n–1,i= f_x(t_n–1,i, x_n–1,i+ ηz_n–1,i), 0 < η <1, i = .

Подставляя последние в (12), получим

ω_n= BD_m _n–1 + AD_k z_n,

где = v"_n-1– y"_n–1, D_m – диагональная матрица с элементами d^m_i,i= l_n–1,i, i = и D_k – диагональная матрица с элементами d^k_i,i= l_n,i.

Для норм матриц в силу (14) справедливы следующие оценки

||D_m|| ≤ L, ||D_k|| ≤ L (15)

Заменим ω_n в уравнении (13) полученным для него выражением, запишем его в виде

z_n= z_n,0e + τDr_n + τDBD_m _n–1+ τDAD_kz_n

Введем норму ||z_n|| = |z_n,i|. Далее, учитывая (11) и (15) и то, что для

= |z_{n–1,k–m+i}| || z_n–1||, получим неравенство

||z_n|| ≤ Сτ^k⁺^m⁺¹ + (1 + kτ||B||)||z_n–₁|| + τkL||A||||z_n||,

которое преобразуем к виду

(1 – τkL||A||)||z_n|| ≤ Сτ^k⁺^m⁺¹+ (1 + kτ||B||)||z_n–₁||.

Если на τ наложить ограничение

τ < τ₀ = , (16)

то оценка, связывающая нормы погрешностей соседних блоков, примет вид

||z_n|| ≤ {Сτ^k⁺^m⁺¹+ (1 + kτ||B||)||z_n–₁||} /(1 – τkL||A||) , (17)

т.к. в силу условия (16) имеет место (1 – τkL||A||) > 0. Подставляя последовательно в (17) значения погрешностей для блоков n–1, n–2, …, 1, получим

≤ /

/ (1 – τkL||A||),

где ||z₀|| = |z_i| – норма погрешности приближенного решения, полученного каким-либо способом на начальном участке интегрирования для первых m узлов. Упростим последнее выражение

||z_n|| ≤

Последнее неравенство справедливо для любого 0 < n ≤ N; из этого следует:

||z_n|| ≤ . (18)

Таким образом, если правая часть уравнения (1) f(tx) удовлетворяет условию Липшица по второму аргументу с константой L и r_n – невязка m-шагового k-точечного блочного метода (2), определенная согласно (18) с оценкой (11), то при

τ < τ₀ = и τkn ≤ T

для погрешности метода имеет место оценка (18).

Следствие. Если разностное уравнение (2) аппроксимирует исходное уравнение (1) и , то решение разностной задачи (2) сходится при τ → 0 к решению исходной задачи (1), причем порядок точности совпадает с порядком аппроксимации (11).

4. Алгоритм решения разностных уравнений многошаговым блочным методом.

Итерационные формулы параллельного решения системы разностных уравнений многошагового многоточечного метода (2) для произвольного блока n получим, использовав (9):

u_n^s+1= u_ne + τDBv"_n–1+ τDA(u"_n)^s, s = , (19)

где s – номер итерации. По формуле (19) пересчитываются компоненты вектора приближенного решения u_n^s⁺¹, принадлежащие блоку n. При этом компоненты вектора , определенные в узлах предшествующего блока, сохраняют свои значения при этих вычислениях, а компоненты вектора вычисляются вновь на каждой итерации. Чтобы начать решение m-шаговым разностным методом, необходимо определить каким-либо одношаговым разностным методом значения приближенного решения u_0,i в первых m-1 точках, примыкающих к начальной точке t₀ отрезка интегрирования. Таким образом будут определены значения компонент вектора ν"₀. Теперь следует задать исходное приближение для компонент вектора u⁰₁, значения которых могут быть, например, найдены по m-шаговым формулам Адамса-Башфорта. В общем случае m-шаговые формулы Адамса-Башфорта для значений приближенного решения в узлах k-точечного блока на первой итерации в векторной форме могут быть записаны в виде

u¹_n= u_ne + τB₁v"_n–1. (20)

Можно показать, что выполнение k шагов вычислений по формуле (19) обеспечит получение приближенного решения с локальной ошибкой порядка .

Приведем результаты решения следующей задачи Коши

x"= –10(t–1)x , x(0 )=1, (21)

трехшаговым четырехточечным методом с шагом t = 0.017.

a b

Рис 2. Графики приближенного решения уравнения (21) трехшаговым четырехточечным методом a и глобальной погрешности b

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

5. Методы с контролем на шаге.

Для оценки локальной погрешности при решении одношаговым многоточечным методом используем следующий подход. Решаются две задачи на одной сетке c одним и тем же шагом t:

первая – одношаговым k-точечным методом;

вторая – одношаговым (k+1)-точечным методом.

Второе решение необходимо для оценки локальной погрешности, поэтому основным является k-точечный метод. Используются значения приближенных решений в совпадающих k узлах основного блока. Лишняя k+1 точка является начальным приближением в расчетах решения в следующем блоке. Локальная погрешность приближенного решения одношаговым k-точечным методом в i-ом узле блока определяется формулой

u_n,i^(k)– x(t_n,i) ≈ φ^(k+2)(t_n,0, x_n,0)τ^k+2, i =

и для (k+1)-точечного метода локальная погрешность в том же узле определяется формулой

u_n,i^(k+1)– x(t_n,i) ≈ φ^(k+3)(t_n,0, x_n,0)τ^k+3, i = .

Вычитая из верхнего соотношения нижнее, получим представление главного члена погрешности k-точечного метода на шаге в виде

(22)

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

Для иллюстрации приведем оценки локальной ошибки, полученные по результатам решения задачи (21) с шагом t = 0.017, одношаговыми трехточечным (основное решение) и четырехшаговым методами (корректирующее решение). Оценка модуля пошаговой погрешности, полученная как разность решений одношаговыми трех- и четырехточечным блочными методами, и действительная локальная погрешности отличаются на величину порядка τ⁵, что обеспечивает надежность используемой оценки.

a b

Рис.3. Графики a- оценки локальной погрешности Etи b-глобальной погрешности Ep при решении задачи (21) одношаговым трехточечным методом

Из сравнения графиков, представленных на рис. 3, следует, что оценка локальной погрешности приближенного решения уравнения (21) трехшаговым четырехточечным методом меньше (в 15 раз) соответствующей накопленной погрешности. Это полностью согласуется с теоретическими оценками.

Для оценки локальной погрешности при решении многошаговым многоточечным методом используется аналогичный подход. Решаются две задачи на одной сетке c одним и тем же шагом t: первая – m-шаговым k-точечным методом; вторая – (m+1)-шаговым k-точечным методом. Второе решение используется для оценки локальной погрешности, поэтому основным является m-шаговый k-точечный метод. Оцениваются значения приближенных решений в совпадающих k узлах основного блока. Локальная погрешность приближенного решения m-шаговым k-точечным методом в i-ом узле блока определяется формулой

(23)

а для (m+1)-шагового k-точечного метода локальная погрешность в том же узле определяется формулой

(24)

Вычитая из (23) соотношение (24), получим представление главного члена погрешности одношагового k-точечного метода на шаге в виде

Таким образом, для главного члена погрешности получаем оценку

(25)

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

Рис. 4. График оценки локальной погрешности при решении уравнения (21) двухшаговым четырехточечным методом.

В рассматриваемом случае оценка локальной погрешности приближенного решения двухшаговым четырехточечным методом практически совпадает со значениями накопленной погрешности. Таким образом, получается завышенная оценка локальной погрешности. Это можно объяснить тем, что решение трехшаговым четырехточечным методом с машинной точностью совпадает с точным решением задачи (21).

6. Заключение. Оценка эффективности блочных алгоритмов.

При численном решении задачи Коши для сравнительной характеристики методов можно рассматривать различные показатели. В случае произвольной правой части уравнения о трудоемкости метода естественно судить по числу обращений для вычисления значений правой части уравнения на каждый узел сетки. Для оценки эффективности одношаговых блочных методов найдем отношение времени выполнения алгоритма Рунге-Кутта на однопроцессорной ЭВМ ко времени выполнения одношагового блочного алгоритма соответствующего порядка на параллельной ВС. Определим время выполнения алгоритма Рунге-Кутта k+1 порядка точности на одном процессоре. Обозначим через t_f– время вычисления значения функции f(tx), t_ad, t_mul – время выполнения операции сложения и умножения соответственно. Время последовательного вычисления приближенных значений решения с точностью O( ) во всех k узлах блока составит

T_s=k (k + 1)² t_f+ k² (t_ad+ t_mu).

Для параллельного выполнения вычислений по формулам (3) закрепим за каждым узлом блока процессор. При его реализации на k процессорах можно одновременно вычислять значения F_n,i,s,, азатем также одновременно получить по формулам (3) значения u_n,i,sдля каждого фиксированного s Объединим процессоры в кольцо, чтобы иметь возможность одновременной передачи данных соседним процессорам. Обозначим через t_ta– время передачи числа соседнему процессору. Время параллельного вычисления приближенных значений решения с той же точностью для всех узлов блока составит

T_p=( k+1)t_f+ k(t_ad+t_mu) + k(k – 1)t_ta.

Ускорение параллельного одношагового k-точечного алгоритма можно будет теперь вычислить по формуле

W(k)= T_s/ T_p.

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

W(k)_pot= k

Для оценки ускорения m-шагового k-точечного блочного метода сравним время его выполнения на мультипроцессорной системе со временем выполнения алгоритмаm-шагового метода Адамса-Башфорта на однопроцессорной ЭВМ. Метод Адамса-Башфорта можно рассматривать как многошаговый одноточечный блочный метод. Последовательное k-кратное применение формул Адамса-Башфорта позволяет вычислить приближенное решение в тех же k узлах блока, в которых параллельно за k итераций может быть вычислено решение m-шаговым k-точечным блочным методом. В этом случае время вычисления будет приблизительно одинаково. Точность приближенного решения, полученного m-шаговым k-точечным блочным методом, имеет порядок , а точность приближенного решения, полученного по m шаговой формуле Адамса-Башфорта, имеет порядок . Поэтому для получения решения с одинаковой точностью для метода Адамса-Башфорта надо выбрать шаг сетки мельче в раз, чем шаг для m-шагового k-точечного метода. Здесь M – число узлов сетки на отрезке решения задачи методом Адамса-Башфорта.

Таким образом, ускорение параллельного m-шагового k-точечного алгоритма равно

W(mk) .

Аналогично могут быть получены оценки эффективности решения задачи Коши для системы обыкновенных дифференциальных уравнений параллельными блочными методами. Так, например, для наиболее употребительной формулы Адамса-Башфорта с m = 4 получим

W(4,4) .

Если на отрезке интегрирования взято сто узлов M = 100, то W(4,4)=5.62.

Литература

Worland P.B. Parallel method for the numerical solution of ordinary differential equations.IEEE Trans.Comp.C-2,10 (1976) 1045-8.
Дж. Холл,Дж. Уатт. Современные численные методы решения обыкновенных дифференциальных уравнений. М.: Мир,1979.-312с.
Системы параллельной обработки. Под ред. Ивенса Д. – М.:Мир,1985. – 416с
Молчанов И.Н. Введение в алгоритмы параллельных вычислений. АН УССР, Инст. Кибернетики им.В.М. Глушкова. Киев: Наукова думка. 1990,-128с
Фельдман Л.П. Параллельные интерполяционные алгоритмы численного решения задачи Коши для обыкновенных дифференциальных уравнений на SIMD компьютере. Наук. пр. ДонДТУ. Серія: Проблеми моделювання i автоматизації проектування динамічних систем, випуск 10: – Донецьк:, 2000, с. 15-22.
L.P. Feldmann Implementierung und Effizienzanalyse von parallelen blockartigen Simulationsalgorithmen für dynamische Systeme mit konzentrierten Parametern. In: Möller, D.P.F. (Hrsg.): Tagungsband 14. ASIM-Symposium Simulationstechnik in Hamburg, September 2000, SCS-Europe BVBA, Ghent/Belgium 2000, S. 241-246.
L.P. Feldman, O.A. Dmitrieva, S. Gerber. Abbildung der blockartigen Algorithmen auf Parallelrechnerarchitekture. In: Tavangarian,D., Grützner,R. (Hrsg.): Tagungs-band 15. ASIM-Symposium Simulationstechnik in Rostock, September 2002, SCS-Europe BVBA, Ghent/Belgium 2002, S.359-364.