Многошаговые методы решения обыкновенных дифференциальных уравнений

Если в один и тот же момент времени выполняются одновременно несколько операций обработки данных, то такие вычисления называются параллельными[1]. Потребность решения сложных прикладных задач с большим объемом вычислений и принципиальная ограниченность максимального быстродействия «классических» — по схеме фон Неймана — ЭВМ привели к появлению многопроцессорных вычислительных систем. Использование таких средств вычислительной техники позволяет существенно увеличивать производительность ЭВМ при любом существующем уровне развития компьютерного оборудования. При этом, однако, необходимо «параллельное» обобщение традиционной — последовательной — технологии решения задач на ЭВМ. Так, численные методы в случае многопроцессорных вычислительных систем должны проектироваться как системы параллельных и взаимодействующих между собой процессов, допускающих исполнение на независимых процессорах. Применяемые алгоритмические языки и системное программное обеспечение должны обеспечивать создание параллельных программ, организовывать синхронизацию и взаимоисключение асинхронных процессов.

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

1. Типовые структуры линий связи параллельных ВС

Топология ВС – это способ соединения компьютеров в единую ВС.

К числу типовых топологий относят[2]:

— Линейка ( linear array or farm ) – представляет собой линейный массив процессоров (рис.1).

Рисунок 1 – Сеть процессоров с топологией линейка

Система, в которой каждый процессор имеет линии связи только с двумя соседними (с предыдущим и последующим) процессорами; такая схема является, с одной стороны, просто реализуемой, а с другой стороны, соответствует структуре передачи данных при решении многих вычислительных задач (например, при организации конвейерных вычислений);

— Кольцо ( ring ) – данная топология аналогична топологии линейка с тем учетом, что первый и последний компьютер соединены (рисунок 2). Связи в могут быть однонаправленными или двунаправленными. [6] Кольцевая структура сохраняет достоинства, а также сокращает максимальное расстояние между процессорами – n /2 [5].

Рисунок 2 – Сеть с топологией кольцо

— Решетка (2D-решетка, матрица, сетка – mesh) – система, в которой граф линий связи образует прямоугольную сетку, то есть процессоры расположены в виде правильной двумерной решетки и каждый процессор (кроме крайних) соединен с четырьмя соседями (рисунок 3)

Рисунок 3 – Сеть с топологией решетка

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

— Гиперкуб ( hypercube ) – данная топология представляет частный случай решетки, когда по каждой размерности сетки имеется только два процессора (то есть гиперкуб содержит n = 2 N процессоров при размерности N ).

Рисунок 4 – Сеть с топологией 3-х мерный гиперкуб

Также может быть схема соединения процессоров в четырехмерном кубе. Эта схема приведена на рисунке 5.

Рисунок 5 – Сеть с топологией 4-х мерный гиперкуб

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

— Клика или полный граф ( completely — connected graph or clique ) – система, в которой связи процессоров образуют полный граф.

Рисунок 6 – Сеть с топологией клика

— Звезда ( star ) – система, в которой все процессоры имеют линии связи с некоторыми управляющим процессором (рисунок 7). Данная топология является эффективной при организации централизованных схем параллельных вычислений.

Рисунок 7 – Сеть с топологией звезда

— «Толстое дерево» («fat-tree», hypertree) – архитектура процессоры в которой локализованы в листьях дерева, в то время как внутренние узлы дерева скомпонованы во внутреннюю сеть. Поддеревья могут общаться между собой, не затрагивая более высоких уровней сети. (рисунок 8, 9)

Рисунок 8 – Сеть с архитектурой «Fat-tree»

Рисунок 8 – Сеть с архитектурой «Fat-tree» (вид сверху)

При этом данная топология является наиболее эффективной.

Реальные высокопроизводительные параллельные вычислительные системы обычно используют несколько различных схем соединения процессоров. Это позволяет сочетать лучшие качества известных топологий и минимизировать недостатки. А также, поскольку способ соединения процессоров друг с другом больше влияет на производительность кластера , чем тип используемых в ней процессоров, то может оказаться более целесообразным создать систему из большего числа дешевых компьютеров, чем из меньшего числа дорогих[8].

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

2. Численные методы решения задач Коши для ОДУ

Обыкновенным дифференциальным уравнением (ОДУ) называется уравнение вида [3] (в него входят: независимая переменная t, неизвестная x и ее производные по t )

где, как правило, обозначается значения неизвестной функции буквой x, независимой переменной — t (и интерпретировать ее как время), производных от x по t — x’, x», . x(m). Также используется сокращенное обозначение J(m)x = (x, x’, . x(m)) — этот вектор называют струей, или джетом m-го порядка функции x в точке t. В дифференциальные уравнения может входить также набор C = (C1, C2, . Cp) произвольных постоянных (параметров) [9].

В свою очередь, порядком ОДУ называется порядок старшей производной искомой функции.

Общим интегралом уравнения называют функцию , связывающую независимую переменную t, искомую функцию x(t) и n констант интегрирования с помощью уравнения

т.е. решение x(t) входитнеявным образом, причем количество констант интегрирования равно порядку ОДУ

Общим решением ОДУ называется функция

связывающая независимую переменную t и n констант интегрирования C i , т.е. решение x(t) определяется явным образом.

Для определения констант интегрирования С i задаются дополнительные условия в количестве, равном количеству констант интегрирования или порядку ОДУ.

Если в дополнительные условия подставить исходное уравнение и решить полученную систему относительно С i , а затем подставить в общий интеграл, то получим частное решение или частный интеграл

Аналогичные процедуры с общим решением (3) дают частный интеграл

Если все дополнительные условия задаются в одной точке x 0 , то совокупность ОДУ и дополнительных условий называют задачей Коши для рассматриваемого ОДУ. В этом случае дополнительные условия называют начальными условиями.

Если дополнительные условия задаются более чем в одной точке, то совокупность ОДУ и дополнительных условий называют краевой задачей Коши для рассматриваемого ОДУ, а дополнительные условия называют граничными или краевыми условиями.

3. Многошаговые методы решения задач Коши для ОДУ

Методы решения ОДУ бывают одношаговые и многошаговые. К одношаговым относятся: метод Эйлера, метод Рунге – Кутты и др., а к многошаговым: линейные многошаговые разностные методы, в том числе методы Адамса-Башфорта и методы Адамса-Мултона.

В многошаговых методах повышение точность вычисления происходит за счет использования информации о поведении решения на предыдущих шагах.

Общая схема построения m-шаговых разностных методов, используемых для приближенного решения задачи Коши [4]

(1)

выглядит следующим образом

(2)

где — числовые коэффициенты, не зависящие от n, причем .Уравнение (2) следует рассматривать как рекуррентное соотношение, выражающее новое значение через найденные ранее значения .

Наиболее простые и наиболее часто встречающиеся вычислительные правила для многошаговых методов имеют вид:

(3)

Среди правил (5) особенно широко известны явные (экстраполяционные) при и неявные (интерполяционные) при методы Адамса.

3.1 Явные многошаговые методы Адамса-Башфорта

Все m –шаговые методы Адамса, полученные из формулы (3), положив в ней получили название явные методы Адамса-Башфорта или э кстраполяционные методы Адамса

(4)

Определить коэффициенты формулы (4) можно следующим образом. Пусть известно приближенное решение в m узлах сетки . Следовательно, в этих точках отрезка известно и значение правой части дифференциального уравнения (1) при i=n-m+1, . n-1,n, причем будет уже функцией только одной переменной . Заменим в (6) функцию интерполяционным многочленом , например, многочленом Ньютона и вычислим значение , проинтегрировав (1) на отрезке . Находим

(5)

Проводя интегрирование, получим разностную схему для решения дифференциального уравнения. Порядок схемы определяется величиной остаточного члена интерполяционного полинома. По существу, интерполяционный многочлен в формуле (5) используется вне области, т.е. в данном случае это экстраполяционный многочлен. Поэтому полученный таким образом метод часто называют экстраполяционным методом Адамса — Башфорта. Получим формулы Адамса — Башфорта для различного числа точек m. Зададим таблицу , заданных значений приближенного решения, по которым можно вычислить значения и составить таблицу .

Рассмотрим случай при m =2. В этом случае соответствующий многочлен Ньютона будет иметь вид:

Подставив его в формулу (7) мы получим:

А после интегрирования формула для вычисления приближенного решения будет иметь вид:

3. 2 Неявные многошаговые методы Адамса-Мултона

Неявные m-шаговые методы Адамса определяются формулой (6) при . Их также называют интерполяционные методы Адамса-Мултона.

(6)

Построение интерполяционных m- шаговых методов Адамса-Мултона может быть осуществлено совершенно аналогично случаю экстраполяционных формул. Для этого следует использовать интеграл (5) , в котором интерполяционный многочлен построен по следующим узлам интерполяции .

Приведем интерполяционные формулы Адамса различных порядков точности:

двухшаговый неявный разностный метод Адамса- Мултона

(7)

трехшаговый неявный разностный метод Адамса – Мултона

(8)

Определение порядка аппроксимации неявных методов Адамса-Мултона выполняется аналогично определению невязки методов Адамса-Башфорта

Погрешность аппроксимации трех шагового метода Адамса – Мултона равна

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

(9)

действительно является уравнением относительно неизвестного значения . Поэтому интерполяционный метод Адамса называют неявным.

На практике обычно не решают уравнение (9), а используют совместно явную и неявную формулы (метод Адамса-Бошфорта-Мултона), что приводит к методу прогноза и коррекции. Одним из широко используемых методов прогноза и коррекции является объединение методов Адамса четвертого порядка точности

(10)

При такой последовательности вычислений этот метод является явным

4 Параллельные методы численного решения задачи Коши для ОДУ

В данной главе будут рассмотрены основы распараллеливания одношаговых и многошаговых алгоритмов методов решения задачи Коши для ОДУ [5].

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

Множество M точек равномерной сетки < t m >, m = и t m =T с шагом ? разобьем на N блоков, содержащих к точек каждый, при этом k NхM. В каждом блоке введем номер точки i = и обозначим через t n , i точку n блока с номером i. Точку t n ,0 назовем началом блока n, а t n , k — концом блока. Очевидно, что имеет место t n , k = t n +1,0 . Условимся начальную точку t 0 = t 1,0 в блок не включать. При численном решении задачи Коши для каждого следующего блока новые k значений функции вычисляются одновременно. Поэтому блочные методы особенно удачно реализуются на параллельных вычислительных системах.

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

, , n =1,2,…

5 Оценка параллелизма многошаговых блочных методов

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

Здесь шаг сетки решения задачи методом Адамса-Башфорта, а шаг сетки решения этой же задачи. Если , например, шаг для решения задачи 4 -шаговым 4 -точечным блочным методом выбран равным =0.01, то для обеспечения той же точности решения методом Адамса-Башфорта необходимо взять шаг =0.0063.В этом случае время решение задачи m -шаговым k -точечным блочным методом будет меньше в 16 раз времени решения методом Адамса-Башфорта. В приведенных здесь оценках не учитывалось время обменов между процессорами при решении на линейке.

ЭВМ и сравним времена решения задачи m -шаговым k -точечным блочным методом на однопроцессорной линейке k процессоров . Закрепим за каждым узлом блока процессор. Используя обозначения, введенные в п. 7.3. получим для времени вычисления на одном процессоре с локальной точностью O( во всех k узлах блока по формулам

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

Ускорение параллельных многоточечных алгоритмов можно вычислиь по формуле

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

и эффективность будет равна

Отсюда следует, что при решении на однопроцессорной ЭВМ решение m -шаговым k -точечным блочным методом потребует меньше времени, чем решение этой же задачи соответствующим методом Адамса — Башфорта.

Выводы

В данный момент все чаще и чаще приходится касаться вопросов вычислительной мощности систем. Одним из способов повышения этой мощности является организация параллельных систем.

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

Общая формулировка многошаговых методов

Многошаговые методы

Решения задачи коши

9.1. Метод “предиктор-корректор” 2

9.2. Общая формулировка многошаговых методов 4

9.3. Устойчивость и сходимость разностных методов 7

9.1. Метод “предиктор-корректор”

В одношаговых методах значение y_n₊₁ определяется лишь значениями t_n, y_n и значением шага приращения аргумента h. В многошаговых методах используется также информация в предыдущих точках y_n_-1, y_n_-2,

Многошаговые методы строятся на равномерной сетке

, где h – шаг сетки.

Обозначим: , .

Линейный m-шаговый метод можно описать формулой:

. (9.1)

Многошаговые методы не являются самостартующими. Для использования m-шагового метода необходимо предварительно задать m предыдущих значений y. Однако после того как многошаговый метод стартует, он работает быстрее, чем одношаговый. Например, при использовании метода Рунге-Кутта-Фельберга на каждом шаге приращения аргумента требуется шесть раз вычислять правую часть дифференциального уравнения . В то же время при использовании многошагового метода на каждом шаге приращения вычисляется только одно новое значение правой части.

Еще одним ограничением использования многошаговых методов является наличие равномерной сетки. Однако это ограничение легко преодолевается, например, с помощью интерполяции.

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

Порядок многошагового метода может выбираться автоматически и динамически изменяться.

Если , то многошаговый метод называется явным.

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

Частным случаем многошаговых методов являются методы Адамса:

(9.2)

Методы Адамса также могут быть явными и неявными.

Трудности использования неявной формулы преодолевается в методе Предиктор–Корректор. В этом методе значение неизвестной функции в каждой точке сетки вычисляется дважды. Сначала вычисляется новое значение неизвестной функции с помощью менее точной явной формулы – формулы предиктор. Затем это значение уточняется с помощью неявной формулы – формулы корректор; при этом для вычисления величины используется найденное значение

Простейшим примером метода Предиктор–Корректор является второй модифицированный метод Эйлера, который можно описать с помощью последовательности формул:

– предиктор,

– корректор.

Метод Эйлера является одношаговым методом. Для построения многошаговых методов вновь используем утверждение: решение задачи Коши эквивалентно решению интегрального уравнения

. (9.3)

Для приближенного решения задачи Коши аппроксимируем подынтегральную функцию в формуле (9.3) с помощью интерполяционного многочлена. Интерполяционный многочлен Ньютона степени m для интерполирования назад от точки x_n можно записать следующим образом

. (9.4)

В формуле (9.4) , а обозначает конечную разность порядка k, отсчитываемую от точки x_i:

, . (9.5)

(9.6)

Найдем вначале формулы для двухшагового метода m=2. В качестве подынтегральной функции в формуле (9.6) будем использовать интерполяционный многочлен Ньютона первой степени.

Для вывода формулы предиктор вычислим интеграл (9.6) на отрезке экстраполяции . После замены переменной интегрирования получим

(9.7)

и .

Подставив найденное значение интеграла в формулу (9.3), получим экстраполяционную формулу – формулу предиктор:

(9.8)

При выводе формулы корректор используем многочлен Ньютона для интерполирования от точки x_n₊₁. Соответственно, в формуле (9.7) все индексы нужно увеличить на 1, а пределы интегрирования изменить на :

(9.9)

Подставив найденное значение интеграла в формулу (9.3), получим интерполяционную формулу – формулу корректор:

(9.10)

С помощью двухшагового метода предиктор-корректор найдем решение начальной задачи:

Точное решение представляет собой экспоненту .

Выберем шаг приращения .

Начальное условие: .

Дополнительное условие: .

С помощью формул (9.8) и (9.10) вычислим значение .

По формуле предиктор получаем: .

По формуле корректор уточняем решение: .

Для сравнения приведем расчеты значения , исходя из точки , по формуле второго модифицированного метода Эйлера:

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

Вывод формулы предиктор.

Подставив найденное значение интеграла в формулу (9.3), получим с учетом выражений (9.5) для конечных разностей формулу предиктор – формулу Адамса-Башфорта (Adams-Bashforth):

(9.11)

Вывод формулы корректор.

Аналогично предыдущему с учетом (9.3) и (9.5) получим формулу корректор – формулу Адамса-Мультона (Adams-Moulton):

(9.12)

Общая формулировка многошаговых методов

Для решения задачи Коши

(9.13)

введем сетку с постоянным шагом и сеточные функции:

– точное решение задачи Коши;

– приближенное решение;

Линейным m-шаговым разностным методом называется система разностных уравнений

, (9.14)

где – числовые коэффициенты, независящие от n, причем .

Уравнение (9.14) представляет собой рекуррентное соотношение, выражающее новое значение через найденные ранее значения .

Заметим, что коэффициенты уравнения (9.14) определены с точностью до множителя. Для устранения неоднозначности будем считать, что выполнено условие нормировки

(9.15)

Метод называется явным, если .

Если , то метод называется неявным. В этом случае для нахождения приходится решать нелинейное уравнение

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

Погрешностью аппроксимации на решении или невязкой разностного метода (9.14) называется функция

(9.16)

Функция невязки получается в результате подстановки точного решения дифференциального уравнения в разностное уравнение (9.14).

Подставим эти разложения в формулу для невязки

Изменим порядок суммирования. При этом во второй сумме будем суммировать по индексу i от 1 до p, соответственно уменьшив всюду индекс i на единицу.

Объединим суммы, выделив в первой сумме отдельно слагаемое с j=0.

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

и при (9.17)

Заметим что, если в системе (9.17) отбросить последние s уравнений, то порядок метода понизится на s.

Итак, получили систему уравнений для коэффициентов a_k, b_k. К этим уравнениям нужно добавить еще условие нормировки (9.15): . Получаем, таким образом p+2 уравнений. Чтобы система уравнений не была переопределена, нужно, чтобы количество уравнений не превышало количество неизвестных.

Для неявного метода имеем 2m+2 неизвестных: . Должно быть:

– для неявного метода порядок аппроксимации не превышает 2m.

Для явного метода имеем 2m+1 неизвестных: . Следовательно и порядок аппроксимации не превышает 2m-1.

Преобразуем систему уравнений для коэффициентов a_k, b_k. С учетом условия нормировки уравнение при можно переписать в виде . Учтем также, что коэффициент a₀ входит только в одно уравнение, так что разрешим это уравнение относительно a₀.

Окончательно получаем систему уравнений для коэффициентов линейного m-шагового разностного метода общего вида:

(9.18)

На практике часто используются методы Адамса, которые описываются формулой

(9.19)

Для методов Адамса в системе (9.18) остаются только уравнения, определяющие коэффициенты b_k.

(условие нормировки) и уравнения: (9.20)

Имеем p уравнений. В случае неявного метода Адамса имеется m+1 неизвестных: . Следовательно, наивысший порядок неявного метода равен p=m+1. В случае явного метода имеем m неизвестных, следовательно, максимальный порядок явного метода Адамса равен p=m.

Рассмотрим варианты двухшаговых методов m=2. Максимальный порядок двухшагового метода равен p=4, так что можно написать следующую систему уравнений.

Ограничившись первыми четырьмя уравнениями этой системы, можно построить методы 2-го порядка. Если использовать пять уравнений, получим методы 3-го порядка. Система из шести уравнений определяет единственный двухшаговый метод 4-го порядка.

Варианты методов второго порядка.

1) Положим , тогда . Получаем явный метод Адамса: .

2) Положив , получим метод Милна: .

3) Положив , построим неявный метод: .

Варианты методов третьего порядка.

1) Положив , получим явный метод .

2) Если выбрать , построим неявный метод .

3) Выбрав , получим неявный метод Адамса .

Метод четвертого порядка получим, используя все шесть уравнений:

Построим трехшаговые методы Адамса.

Явный трехшаговый метод имеет порядок p=3. В соответствии с формулами (9.20) для нахождения коэффициентов b_k имеем систему уравнений, включающую условие нормировки и еще два уравнения для двух значений i: .

Решение этой системы в среде Mathcad приведено на рис. 9.1. Здесь M – матрица системы, C – столбец свободных коэффициентов. Функция lsolve находит решение системы линейных алгебраических уравнений. Вектор B включает найденное решение. Чтобы получить точные значения найденных коэффициентов используем знак аналитических вычислений (стрелка). Получаем явный трехшаговый метод Адамса третьего порядка:

Неявный метод имеет четвертый порядок, и для коэффициентов имеем систему из трех уравнений, где i принимает значения 2, 3 и 4. Значение b₀ определяем из условия нормировки. Решение в системе Mathcad приведено на рис. 9.2. Получаем вновь формулу Адамса-Мультона, которая ранее была найдена путем непосредственного вычисления интеграла от интерполяционного многочлена Ньютона (9.12):

Пример 9.4. На рис. 9.3 приведено решение системы уравнений для коэффициентов явного и неявного четырехшаговых методов Адамса. Решение аналогично решению в примере 9.2 для трехшаговых методов. Обозначения M, C и B относятся к явному методу; обозначения MI, CI, BI и bI₀ используются для описания неявного метода. Получаем в итоге:

явный метод четвертого порядка (формула Адамса-Башфорта)

;

неявный метод пятого порядка

источники:

http://poisk-ru.ru/s28953t3.html