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

МГТУ ГА

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

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

на тему:

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

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

МОСКВА 2008

Содержание

Введение

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

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

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

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

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

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

Заключение

Литература

Приложения

Приложение 1

Приложение 2

Приложение 3

Введение

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

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

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 неизвестных параметров l1, l2, ..., ln, которые называются собственными значениями. Для единственности решения на интервале [a,b] необходимо задать n + m граничных условий.

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

(2.3)

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

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

1) Выбираем в области интегрирования упорядоченную систему точек a=x1<x2<…<xn<b, называемую сеткой. Точки xi называют узлами разностной сетки, разность между соседними узлами h=xi-xi-1 – шаг сетки. Формула для вычисления шага равномерной сетки, заданной на интервале [a,b]:

, (2.4)

где nx – количество узлов заданной сетки.

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

Для получения конечно-разностной схемы удобно использовать интегроинтерполяционный метод, согласно которому необходимо проинтегрировать уравнение (2.3) на каждом интервале [xk, xk+1] и разделить полученное выражение на длину этого интервала:

(2.5)

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

ε(h)=|| || (2.6)

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

, (2.7)

где С0 – константа устойчивости, p – порядок апроксимации.

Поэтому для увеличения точности решения необходимо уменьшить шаг сетки h.

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

Одношаговые схемы

– Метод Эйлера

Заменяем интеграл в правой части уравнения (2.5) по формуле левых прямоугольников:

(2.8)

Получим:

, (2.9)

где k=0,1,2,…,n.

Схема явная устойчивая. В силу того, что формула для левых прямоугольников имеет погрешность второго порядка, точность ε(h) первого порядка.

– Неявная схема 1-го порядка

Используя формулу правых прямоугольников, получим:

(2.10)

Эта схема неразрешима в явном виде относительно , поэтому проводится итерационная процедура:

, (2.11)

где s=1,2,… - номер итерации. Обычно схема сходится очень быстро – 2-3 итерации. Неявная схема первого порядка эффективнее явной, так как константа устойчивости С0 у неё значительно меньше.

– Метод Эйлера-Коши

Вычисления проводятся в два этапа : этап прогноза и этап коррекции.

На этапе прогноза определяется приближенное решение на правом конце интервала по методу Эйлера:

(2.12)

На этапе коррекции, используя формулу трапеций, уточняем значение решения на правом конце:

(2.13)

Так как формула трапеций имеет третий порядок точности, то порядок погрешности апроксимации – равен двум.

– Неявная схема 2-го порядка (метод Эйлера-Коши)

Используя в (2.5) формулу трапеций, получим:

(2.14)

Схема не разрешена в явном виде, поэтому требуется итерационная процедура:

, (2.15)

где s=1,2,… – номер итерации. Обычно схема сходится за 3-4 итерации.

Так как формула трапеций имеет третий порядок точности, то погрешность апроксимации – второй.

Схемы с дробным шагом

– Схема предиктор-корректор (Рунге-Кутта) 2-го порядка

Используя в (2.5) формулу средних, получим:

,(2.16)

где – решение системы на середине интервала [xk, xk+1] . Уравнение явно разрешено относительно , однако в правой части присутствует неизвестное значение . Поэтому сначала расчитывают (предиктор):

. (2.17)

Затем расчитывают (корректор) по формуле (2.16). Схема имеет первый порядок погрешности.

– Схема Рунге-Кутта 4-го порядка

Используя в (2.5) формулу Симпсона, получим:

(2.18)

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

Сначала рассчитывают предиктор вида:

(2.19)

затем корректор по формуле:

(2.20)

Поскольку формула Симпсона имеет пятый порядок погрешности, то точность ε(h) – четвёртого порядка.

Многошаговые схемы

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

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

Если заменим в (2.5) подинтегральное выражение интерполяционным многочленом Ньютона, построенного по узлам , то после интегрирования на интервале получим явную экстраполяционную схему Адамса. Если заменим в (2.5) подинтегральное выражение на многочлен Ньютона, построенного по узлам , то получим неявную интерполяционную схему Адамса.

– Явная экстраполяционная схема Адамса 2-го порядка

(2.21)

Схема двухшаговая, поэтому необходимо для расчётов найти по схеме Рунге-Кутта 2-го порядка , после чего , , … вычисляют по формуле (2.21)

– Явная экстраполяционная схема Адамса 3-го порядка

(2.22)

Схема двухшаговая, поэтому необходимо сперва найти и по схеме предиктор-корректор 4-го порядка, после чего , , … вычисляют по формуле (2.22).

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

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

(3.1)

интерполяционным многочленом Ньютона 2-го порядка, вида:

(3.2)

После интегрирования полученного выражения на интервале , приходим к уравнению неявной схемы Адамса 3-го порядка:

. (3.3)

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

, (3.4)

где s=1,2,3,… – номер итерации. Условие выхода из цикла итерационной процедуры:

, (3.5)

где ε – заданная погрешность.

Начальное приближение задаётся формулой для явной экстраполяционной схемы Адамса 2-го порядка:

. (3.6)

Схема устойчива, сходится быстро. Чаще всего достаточно одной итерации. Порядок погрешности ε(h) неявной схемы Адамса третьего порядка равен четырём.

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

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

Блок-схема алгоритма условно разбита на 11 блоков.

Главная функция программы (блоки 1,2,5) отвечает за обработку события создания формы, взаимодействие со стандартным компонентом TСhart, а также за реализацию решения системы дифференциальных уравнений неявной схемой Адамса 3-го порядка. Блок-схема алгоритма решения задачи Коши разбита условно на 35 блоков:

1-й блок отвечает за ручной ввод интервала [a,b], на котором ищется решение системы; количества шагов сетки nx; шаг вывода результатов на экран np; строк u1 и v1, соответствующих уравнениям системы; значения искомых функций в начале заданного интервала; допустимая погрешность e.

Во втором блоке происходит вычисление шага h и установка текущего узла на x=a. Блок 3 – функция преобразования исходных записей уравнений системы в равносильные им строки с постфиксной формой записью математических операций (см. далее «алгоритм обратной польской записи»). В качестве аргументов функции выступают введённые ранее строки u1 и v1.

Блоки 4-15 – расчет первых 2-х точек заданной сетки методом Рунге-Кутта 4-го порядка. В данных блоках и далее используется пользовательская функция FPR, рассчитывающая значения вводимых пользователем уравнений в узлах заданной сетки. В качестве аргументов функции выступают: уже преобразованные в обратную польскую запись строки, задающие уравнения системы; текущее значение x; значения искомых функций на предыдущем шаге (условно обозначаем ).

В блоках 16-34 в цикле (16) рассчитываются значения искомых решений в узлах 2-nx заданной сетки неявной схемой Адамса 3-го порядка. Цикл 21-29 осуществляет итерационную процедуру неявной схемы. Условие выхода из этого цикла – выполнение неравенства de<e, где de – наибольший из модулей , e – заданная точность. Поскольку на экран выводятся значения искомых функций не во всех узлах, а только в узлах с номером, кратным шагу вывода nx, вводимым с клавиатуры, то блоки 33-34 осуществляют выбор этих узлов.

Преобразование в обратную польскую запись происходит по следующим правилам:

Рассматриваем поочередно каждый символ:

1. Если этот символ - число (или переменная), то просто помещаем его в выходную строку.

2. Если символ - знак операции (+, -, *, / ,^), то проверяем приоритет данной операции. Операция возведения в степень имеет наивысший приоритет (равен 4). Операции сложения и вычитания имеют меньший приоритет (равен 2). Наименьший приоритет (равен 1) имеет открывающая скобка.

Получив один из этих символов, мы должны проверить стек:

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

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

3. Если текущий символ - открывающая скобка, то помещаем ее в стек.

4. Если текущий символ - закрывающая скобка, то извлекаем символы из стека в выходную строку до тех пор, пока не встретим в стеке открывающую скобку (т.е. символ с приоритетом, равным 1), которую следует просто уничтожить. Закрывающая скобка также уничтожается.

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

Согласно этим правилам создан модуль ”Unit3.cpp”, содержащий функцию преобразования строки в обратную польскую запись OPZ (блок 3 в блок-схеме алгоритма), алгоритм которой приведён в приложении. В данном модуле использованы также вспомогательные функции PUSH, PRIOR, DEL. Функция PUSH записывает в стек, на вершину которого указывает HEAD, символ a. Возвращает указатель на новую вершину стека. Функция PRIOR вычисляет приоритет текущего символа, естественно, лишь в том случае, если текущий символ – математическая операция. Функция DEL удаляет символ с веpшины стека. Возвpащает удаляемый символ. Изменяет указатель на веpшину стека.

Для работы с полученной обратной польской записью создана функция(блок 4), организованная в виде подключаемого модуля