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

Лабораторная работа по выч. математике №4 «Решение обыкновенных дифференциальных уравнений. Метод Адамса»

CАНКТ-ПЕТЕРБУРГСКИЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ, МЕХАНИКИ И ОПТИКИ

Лабораторная работа по выч. математике №4

«Решение обыкновенных дифференциальных уравнений.

Выполнил: Припадчев Артём

Задание: составить подпрограмму для решения ОДУ первого порядка используя многошаговый метод Адамса. Разгонные точки вычислить методом Рунге-Кутта 4-го порядка. Вычисление правых частей реализовать отдельной подпрограммой. Найти решение заданного уравнения с точностью e, контролируя точность на каждом шаге вычислений, построить график решения.

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

В формуле Симпсона для приближенного вычисления определенного интеграла используются значения подинтегрального выражения в трех точках. В интеграле их всего две, поэтому введем дополнительную точку в середине отрезка [xi+1 xi]

тогда можно переписать так:

Полученное выражение является неявным, так как в правой части содержатся еще не определенные значения функции yi+h/2 и yi+1. Чтобы воспользоваться этой формулой, надо использовать некоторое приближение для вычисления этих значений

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

Алгоритм Рунге-Кутты третьего порядка — РК3 (погрешность порядка h3):

(6.8)

Алгоритм Рунге-Кутты четвертого порядка — РК4 (погрешность порядка h4):

(6.9)

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

Рассмотренный ранее метод Рунге-Кутты использует значение функции на одном предшествующем шаге, поэтому они относятся к так называемым одношаговым методам. Точность вычислений можно увеличить, если использовать при нахождении решения в некотором узле xi информацию о значениях функции, полученных в нескольких (k) предыдущих узлах сетки интегрирования (xi-1, xi-2 … xi-k).

Если используются значения в k предыдущих узлах, то говорят о k-шаговом методе интегрирования уравнения. Одним из способов построения многошаговых методов заключается в следующем. По значениям функции, вычисленным в k предшествующих узлах, строится интерполяционный полином степени (k-1) — , который используется при интегрировании дифференциального уравнения по выражению (6.3). Интеграл при этом выражается через квадратурную формулу:

,

где лl – квадратурные коэффициенты.

Очевидно, что при k=1 в качестве частного случая получается формула Эйлера. Значения квадратурных коэффициентов для k от 2 до 4 приведены в таблице.

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

Автор: Дмитриева O.А.
Донецкий национальный технический университет
83000, Донецк, ул. Артема, 58
e-mail: dmitriv@r5.donntu.org
http://www.imamod.ru/magazin/pdf/12/1205-081r.pdf

В течение последних двух-трех десятилетий пиковая производительность параллельных вычислительных систем возросла на несколько порядков. Радикально изменились технологическая база и архитектура. Но главным препятствием к внедрению практически всех параллельных архитектур является отсутствие параллельных алгоритмов. Речь, в первую очередь, идет о больших системах дифференциальных уравнений, как обыкновенных, так и в частных производных, являющихся основой больших современных научных и инженерных задач. Многие численные методы решения приходится пересматривать, а от некоторых полностью отказываться. В то же время, целый ряд вопросов, которые были несущественны или вообще не возникали при проведении последовательных вычислений, приобрели исключительную важность для эффективного и правильного использования вычислительных систем. Опыт эксплуатации параллельных систем показал [1]. что для эффективного их применения нужно радикально менять структуру численных методов. Эти обстоятельства послужили основой для адаптации известных методов решения систем обыкновенных дифференциальных уравнений (ОДУ) на параллельных системах с SIMD архитектурой.

Численное решение систем ОДУ методами Адамса-Башфорта и Адамса-Моултона на решетке процессоров

Численное решение задачи Коши для системы обыкновенных дифференциальных урав¬нений с постоянными коэффициентами

можно получить последовательно по шагам с помощью, например, следующих формул Адамса-Башфорта и Адамса-Моултона [2]

где h — выбранный шаг интегрирования.

Приведенный метод в целом является явным. Сначала по формуле Адамса-Башфорта вычисляется значение х_k+1 (р) являющееся прогнозом для х_k+1. Затем х_k+1 (р) используется для расчета скорректированного значения х_k+1 (k) , вычисляемого по формуле Адамса-Моултона. Представляемый метод имеет четвертый порядок точности, хотя можно получить эти методы сколь угодно высокого порядка, используя все большее число предыдущих точек, а, следовательно, интерполяционный полином более высокой степени. При этом с ростом точности формулы становятся все более громоздкими, но принцип остается тем же. Многошаговые методы порождают проблему, которая не возникала при использовании одношаговых методов, поскольку для начала счета им необходимо несколько разгоночных значений. При этом, однако, существенно, чтобы эти стартовые значения были вычислены с той же степенью точности, с которой будет работать многошаговый метод. Для этой цели предполагается использование метода Рунге-Кутты четвертого порядка точности, подробное описание реализации ко¬торого на параллельных архитектурах типа SIMD представлено в [3].

Модель, на которую ориентируется решение, имеет следующие особенности: используется вычислительная система SIMD структуры с квадратной сеткой, содержащей m*m процессорных элементов (ПЭ); каждый процессорный элемент может выполнить любую арифме¬тическую операцию за один такт; временные затраты, связанные с обращением к запоминаю¬щему устройству, отсутствуют. Для простоты изложения рассматривается случай, когда количество процессорных элементов в каждой строке совпадает с размерностью задачи. Для эффективной работы методов Адаме а в каждый процессорный элемент, имеющий индексы i,j, пересыпается элемент исходной матрицы матрицы А, приведенной предварительно к виду (3).

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

Из вектора неизвестных х формируется новая матрица U (см. (5)), элементы которой также пересылаются в соответствующие процессорные элементы. Таким образом, в каждом процессорном элементе с индексами i,j оказывается соответствующее значение g_i,j и элемент матрицы u_i,j. Первоначально выполняется умножение g_i,j*U_i,j во всех процессорных элементах. Затем осуществляется одиночный циклический сдвиг полученных значений и сложение содержимого процессорных элементов. Количество позиций, на которое производится очередной сдвиг, определяется элементами следующего ряда: 2 0 ,2 1 ,2 2 . 2 k-1 , где k — ближайшее целое сверху [log₂m]. Учитывая, что в процессорных элементах уже находятся рассчитанные ранее значения для предыдущих точек, произведем необходимые операции умножения и сложения.

Легко видеть, что на такой сетке в каждом ПЭ i-й строки получается новое значение x_i (p) на (n+1)-м шаге. Вычисление всех значений x_n+1 (p) определяется временем, затрачиваемым на одно умножение t_умн, суммой времен, необходимых для осуществления сдвигов и сложения сдваиванием , а также временами 4t_умн — умножение хранящихся значений на новые коэффициенты, 3t_сл — суммирование полученных результатов, 4t_умн+3t_сл — время для умножения на h/24 и суммирования со значением, полученным на предыдущем шаге. Всего получаем

Но необходимо учесть, что это время характеризует только нахождение значений век¬тора неизвестных прогнозируемых значений. При этом, получающиеся в сетке процессорных элементов значения расположены не в том порядке, как при изначальной засылке (5), а имен¬но в виде (6)

Для того, чтобы выполнить расчет скорректированных значений, нужно потратить какое-то время на переупорядочивание элементов, для приведения их к виду (5). На это потребуется m-1 циклических сдвигов по столбцам (в первом — 0, во втором — 1, . в m столбце необходимо выполнить 2 k-1 сдвигов). Итого дополнительно потребуется времени (m-1)*t_сдв. После этого для получения скорректированных значений нужно повторить ту же последовательность операций, но с другими значениями. Общее время, затрачиваемое на нахождение вектора неизвестных х_n+1 на решетке процессоров m*m, оценивается как

Для оценок качества рассмотренного алгоритма используются наиболее распространенные критерии [4]: коэффициент ускорения и коэффициент эффективности , где T_i(N) — время вычисления вектора x_n+1 на однопроцессорной ЭВМ с быстродействием арифметического процессора и объемом ОЗУ, равным суммарному объему всех ЗУ арифметических процессоров, и с необходимым числом внешних устройств, имеющих скорости обмена такие же, как в многопроцессорной вычислительной системе типа SIMD. При расчете значений вектора x_n+1 на однопроцессорной ЭВМ потребуется время, равное

Чтобы оценить параллелизм алгоритма, можно принять [4], что t_оп=t_умн=t_сл=10t_сдв. Тогда

Эффективность методов Адамса при решении на SIMD структуре

Полученные результаты значительно отличаются от характеристик потенциального параллелизма методов Адамса (ускорение ускорение примерно равно m 2 , эффективность — 1), поскольку при расчете этих характеристик не учитывались времена, затрачиваемые на сдвиги, хотя, как оказалось, эти времена существенно влияют на характеристики параллелизма (рис.1).

Реализация методов Адамса на линейке процессоров

Еще одним часто встречающимся способом коммутации параллельных вычислительных систем является линейка процессоров. Рассмотрим решение задачи (1) на модели SIMD-структуры, построенной из последовательно соединенных процессорных элементов (последний связан с первым) при условиях, оговоренных выше. Из всех рассмотренных способов первоначальной засылки значений в процессорные элементы оптимальным оказался способ, когда в 1-й процессорный элемент записываются значения i-й строки матрицы А: a_i1,a_i2. a_im и значения векторов неизвестных, полученных на предыдущих шагах, первоначально: x₍₀₎, x₍₁₎, x₍₂₎, x₍₃₎ , а на (n+1)-м шаге x_(n), x_(n-1), x_(n-2), x_(n-3). При таком подходе одновременно все процессоры могут начать проводить вычисления прогнозируемых значений по формулам (2). Через время m*t_умн+(m-1)t_сл+4t_умн+3t_сл+t_умн+t_сл=2m*t_оп+8t_оп в каждом i-м процессорном элементе получаем прогнозируемое значение x_n+1(i) (p) . Эти действия необходимо выполнить только один раз, перед началом основных вычислений и хранить в каждом i-м процессорном элементе получившиеся значения:

Теперь необходимо передать полученные прогнозируемые значения по типу «все-всем», так как в каждом процессорном элементе нам необходимо иметь значения векторов x_n+1 (p) . Время, затрачиваемое на обмены, определится, как (m-1)t_сдв, после чего в каждом ПЭ будут содержаться векторы прогнозируемых значений в виде (11) и потребуется еще время на переупорядочивание элементов и приведение их к виду (12). При этом значения, хранящиеся в первом ПЭ останутся без изменений, во втором на их переупорядочивание потребуется время (m-1)t_сдв, в третьем (m-2)t_сдв в последнем t_сдв. Учитывая, что все сдвиги могут выполняться одновременно, нужное расположение значений получается за время (m-1)t_сдв, тогда общее время, затрачиваемое на распространение «все-всем» осуществится за время 2(m-1)t_сдв.

Вторая группа сдвигов, потребовавшаяся на переупорядочивание значений, может быть исключена, если провести первоначальную перестановку элементов в строках, т.е. записать в каждый (i-Й процессорный элемент i-ю строку матрицы G(3). Тогда после получения в каждом ПЭ прогнозируемого значения, передача «все-всем» может быть осуществлена за время (m-1)t_сдв. Перед расчетом скорректированного значения х_n+1 (k) из ПЭ удаляются значения х_n-3 как уже не использующиеся. Для расчета скорректированного значения потребуется время m*t_умн+(m-1)t_сл+4t_умн+3t_сл+t_умн+t_сл=2m*t_оп+8t_оп , рассылка «все-всем» (m-1)t_сдв Приняв t_сдв=0,1*t_оп, можно сказать, что полный цикл расчета значений для одной точки закончится через

Ускорение определим как

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

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

1. Воеводин В.В. Информационная структура алгоритмов. -М.: Изд-во МГУ, 1997, 139с.

2. Ортега Дж„ Пул У. Введение в численные методы решения дифференциальных уравнений / Пер. с англ.; Под ред. Абрамова А.А.-М.: Наука, Гл. ред. физ.-мат. лит., 1986, 288 с.

3. Фельдман Л.П. Параллельные алгоритмы моделирования динамических систем, описываемых обык¬ новенными дифференциальными уравнениями. Научн. тр. ДонГТУ, Вып. 6, Донецк: ДонГТУ, 1999, 330с.

Курсовая работа: Решение систем дифференциальных уравнений при помощи неявной схемы Адамса 3-го порядка

Кафедра вычислительной математики и программирования

Пояснительная записка к курсовому проекту

«Решение систем дифференциальных уравнений при

помощи неявной схемы Адамса 3-го порядка»

1. Постановка задачи

2. Описание математических методов решения

3. Описание используемого метода

4. Описание блок-схемы

5. Описание программы

6. Анализ результатов

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

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

1. Постановка задачи

Необходимо решить с заданной степенью точности задачу Коши для системы дифференциальных уравнений на заданном интервале [a,b]. Добиться погрешности на втором конце не более 0,0001. Результат получить в виде таблицы значений приближенного и точного решений в точках заданного интервала. Построить графики полученных решений и сравнить их с точным решением.

– система дифференциальных уравнений вида:

– интервал, на котором ищется решение: [a,b]

– погрешность, с которой ищется решение: е

– формулировка задачи Коши в начальной точке заданного интервала: u(a)=u, v(a)=v

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

– шаг вывода на экран значений искомых функций в узлах заданной сетки: np

– таблица значений приближенного и точного решений в узлах заданной сетки;

– графики полученных и точных решений.

2. Описание математических методов решения задачи

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

(2.1)

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

Первый тип – это задачи Коши, или задачи с начальными условиями. Для таких задач кроме исходного уравнения в некоторой точке a должны быть заданы начальные условия, т.е. значения функции u1(a),…, um(a):

u1(a)=,…, um(a)= (2.2)

Ко второму типу задач относятся так называемые граничные, или краевые задачи, в которых дополнительные условия задаются в виде функциональных соотношений между искомыми решениями. Количество условий должно совпадать с порядком n уравнения или системы. Если решение задачи определяется в интервале xÎ[a,b], то такие условия могут быть заданы как на границах, так и внутри интервала.

Третий тип задач для систем дифференциальных уравнений – это задачи на собственные значения. Такие задачи отличаются тем, что кроме искомых функций u1(x),…, um(x) в уравнения входят дополнительно n неизвестных параметров l₁ , l₂ , . l_n , которые называются собственными значениями. Для единственности решения на интервале [a,b] необходимо задать n + m граничных условий.

Рассмотрим подробнее задачу Коши. Воспользуемся компактной записью задачи (2.1), (2.2) в векторной форме:

(2.3)

Требуется найти на интервале [a,b].

Задачу Коши удобнее всего решать методом сеток. Метод сеток состоит в следующем :

1) Выбираем в области интегрирования упорядоченную систему точек a=x1

#pragma resource «*.dfm»

char *opz(char *); // ф-ия преобразования в обратную польскую запись;

double fpr(char *str,double u, double v,double x); // обратныйходпольской

int p=1,s=1,j=1,o=0; // записи;

__fastcall TForm1::TForm1(TComponent* Owner)

void __fastcall TForm1::N5Click(TObject *Sender)

void __fastcall TForm1::Button3Click(TObject *Sender)

void __fastcall TForm1::N7Click(TObject *Sender)

void __fastcall TForm1::N2Click(TObject *Sender) // очисткаформы