Методы и алгоритмы компьютерного решения
дифференциальных уравнений
Введение
Для того, чтобы описать динамику различных процессов, протекающих в природных и в технических системах, составляют, опираясь на физические законы, дифференциальные уравнения. Так, в частности, приходится поступать при исследовании функционирования автоматических систем; работы судовых энергетических комплексов, электрических агрегатов, судовых вспомогательных механизмов, систем навигации и т.д. В ряде случаев эти уравнения допускают линеаризацию и могут быть записаны в виде:
 ,
где y
(t
) – неизвестная функция, a
0
, a
1
,...a
n
– постоянные коэффициенты, а j(x
) – некоторая известная функция независимого аргумента t
, которая обычно выражает внешнее воздействие, оказываемое на систему.
1. Цель контрольной работы
Приобретение навыков алгоритмизации и программирования задач численного интегрирования обыкновенных дифференциальных уравнений и систем с последующим моделированием результатов на персональном компьютере и представлением их в виде таблиц и графиков.
В результате выполнения контрольной работы студент обязан:
1. Научиться решать линейные дифференциальные уравнения численными и символьными методами в рамках пакета компьютерной математики MathCAD.
2. Ознакомиться с основными алгоритмами существующих компьютерных методов.
3. Определить точность этих методов путем сравнения результатов, получаемых путем приближенного и аналитического решений.
2
. Аналитические методы
Общее решение дифференциального уравнения n-го порядка – неизвестная функция y
(t
) – содержит n
произвольных постоянных. Их можно определить, зная начальные условия, накладываемые на неизвестную функцию и на ее производные вплоть до (n-1)-порядка включительно. Аналитически (в символьном виде) такие уравнения решают классическим и операционным методами.
2
.1 Классический метод
В ограниченном числе случаев вида левой части (1) допускает такое преобразование, которое позволяет найти решение путем непосредственного интегрирования, однако в общем случае порядок решения – иной.
Решение неоднородного дифференциального уравнения (с ненулевой правой частью) является суммой общего решения соответствующего однородного дифференциального уравнения y
1
(t
) и частного решения y
2
(t
) неоднородного дифференциального уравнения (1).
Решение однородного уравнения ищем в виде:  . Подстановка его в дифференциальное уравнение приводит к характеристическому алгебраическому уравнению n
-ного порядка:
 ,
которое имеет n корней –  . В частном случае отсутствия кратных корней общее решение может быть записано в виде:
 ,
где Сi
– произвольные постоянные, которые находятся из начальных условий.
Имеются правила, позволяющие определить вид y2
(t
) частного решения в зависимости от вида правой части – функции j(t
).
Последующая подстановка общего решения в исходное дифференциальное уравнение позволяет найти неопределенные константы Ci
в выражении для y
1
(t
).
«Классический» метод анализа процессов в настоящее время используется только в случае простейших систем, поскольку необходимость нахождения частного решения часто приводит к сложным преобразованиям, а также, кроме решения характеристического уравнения дополнительно необходимо составить и решить n
уравнений для определения постоянных интегрирования.
2
.2 Метод операционного исчисления
Суть метода состоит в проведении интегрального преобразования Лапласа функции, входящей в состав дифференциального уравнения, по правилу:
 ,
где s =
a
+ j
×
b
– комплексная переменная величина.
Это преобразование сопоставляет функции действительного переменного функцию комплексного переменного . При этом для линейных дифференциальных уравнений существует изоморфизм
(взаимно-однозначное соответствие) между функциями-оригиналами, входящими в уравнение, и их изображениями (образами Лапласа).
Преобразование Лапласа можно выполнить, используя блок символьных вычислений MathCAD. Этот же блок позволяет выполнить и обратное преобразование Лапласа, в соответствии с соотношением:
 ,
где  , т. е. интегрирование проводится по прямой, лежащей в плоскости комплексного переменного s
и проходящей параллельно мнимой оси jw на расстоянии s
от нее, при этом Лаплас образ Y
(s
) должен иметь особенности слева от этой линии.
Преобразование Лапласа сводит дифференцирование функции оригинала к умножению ее образа на комплексную переменную s
, поэтому решение дифференциального уравнения в пространстве оригиналов сводится к решению алгебраического уравнения в пространстве изображений.
Порядок решения дифференциального уравнения с помощью операционного исчисления представляется следующим:
- выполняя преобразование Лапласа левой и правой части дифференциального уравнения, учитываем начальные условия и переходим от дифференциального уравнения для функции оригинала y
(t
) к алгебраическому уравнению для Лаплас образа – Y
(s
) ;
- решая алгебраическое уравнение, находим в пространстве изображений в явном виде выражение для Y
(s
);
- выполняя обратное преобразование мы находим неизвестную функцию y
(t
).
Все этапы этой процедуры могут быть автоматизированы и выполнены в рамках пакета MathCAD (пример 1).
Следует заметить, что пакет MathCAD далеко не всегда способен выполнить в символьной форме результат обратного Лаплас преобразования. Дело в том, что в блок символьных преобразований пакета заложены правила выполнения данной процедуры для выражений записанных в виде элементарных дробей. Поэтому Лаплас образ предварительно разлагается на элементарные дроби. Однако, если корни полинома в знаменателе представляются в виде комбинации сложных радикалов, то MathCAD «отказывается» работать. В этом случае ему необходимо «помочь» врукопашную выполнив разложения полинома в знаменателе в соответствии с соотношением:
 ,
где s
1
, s
2
,…sn
– корни уравнения  . В примере 1 рассмотрено выполнение обратного преобразования Лапласа и для такого случая.
Рассмотренная методика нахождения аналитического решения дифференциальных уравнений может быть распространена на задачу решения системы дифференциальных уравнений. В этом случае необходимо решить не одно алгебраическое уравнение для Лаплас-образов, а систему алгебраических уравнений с помощью той же процедуры блока решений Given – Find. Отметим, что в отличие от систем компьютерной математики Mathematica 2.2.2 и Maple V R3/R4, которые легко позволяют аналитически решить линейное дифференциальное уравнение с помощью встроенных средств. Система MathCAD предполагает «ручные процедуры» запуска прямого преобразования Лапласа, составления по его результатам алгебраического уравнения и, после его решения, запуска процедуры обратного преобразования Лапласа.
3
. Численный метод решения дифференциальных уравнений
Часто при анализе изучаемых в технических или в природных системах процессов приходится учитывать наличие нелинейного поведения функций, описывающих характеристики их элементов. Это в свою очередь определяет появление нелинейностей в дифференциальных уравнениях, которые теперь уже не могут быть записаны в форме (1). Наличие нелинейностей в дифференциальных уравнениях обуславливает невозможность их точного аналитического решения, а приближенные аналитические методы часто приводят к громоздким выкладкам.
Кроме того, коэффициенты в левой части дифференциального уравнения могут быть определены из эксперимента с ошибкой, что, в значительной степени, обесценивает получаемый точный аналитический результат.
И наконец, точные методы не пригодны для случая, если правая часть дифференциального уравнения представлена не в аналитической форме, а в виде таблицы или графика.
Во всех этих случаях прекрасно продолжают «работать» методы численного решения. В отличие от аналитических, они позволяют получать искомые зависимости для любой из описанных выше ситуаций. Алгоритмы существующих методов численного решения были разработаны сравнительно давно, однако толчок к их применению был обусловлен развитием вычислительной техники. Каждый из существующих численных методов предполагает замену производной на конечное приращение и преобразование дифференциального уравнения в уравнение в конечных разностях.
С этой целью интервал поиска решения разбивается на множество отрезков и решение ищется на каждом из этих кусочков. Ясно, что чем мельче шаг разбиения, тем точнее получается результат. Поэтому, эффективное применение численных методов (при решении реальных, а не учебных задач!) предполагает использование компьютеров с достаточным быстродействием.
Использование для численного решения дифференциальных уравнений компьютерного пакета MathCAD предполагает знание алгоритма работы численных методов для разумного их применения (знание границ применимости, оценки точности, затрат компьютерных ресурсов и др.). Дело в том, что к результатам компьютерных вычислений всегда нужно относиться критически; анализировать их на правдоподобность, и для того, чтобы избежать "подводных камней" при использовании любого стандартного пакета, реализующего численные методы, нужно иметь хотя бы минимальное представление о том, какой именно численный метод реализован для решения той или иной задачи и как он «работает».
Поскольку для оценки точности решения необходим материал для сравнения предлагается рассматривать приближенные (численные) методы решения применительно к тем задачам, для которых ранее было получено аналитическое решение (т.е. линейных уравнений вида (1)).
Обозначим   . Тогда данное уравнение можно преобразовать в следующую систему уравнений первого порядка, разрешенных относительно первой производной (форма Коши):
Поэтому, идея применения численных методов для решения уравнения старших порядков, в принципе, ничем не отличается от идеи численного решения уравнения первого порядка, которые рассматривались в лабораторной работе №10.
Нетрудно расширить применение описанной выше методики на случай системы линейных уравнений. В примерах 2 и 3 представлены реализации метода Эйлера в рамках векторной процедуры и с помощью программы-функции. В обоих случаях получены одинаковые результаты, которые поверяются по полученному ранее (см. пример 1) аналитическому решению. Решение приближенное и точное сильно отличаются и величина относительной ошибки (для выбранного шага) составляет ~13%. Пример 3 нетрудно оптимизировать и придать ему более компактный вид, считая начальные условия и правые части системы уравнений компонентами некоторых векторов. Также допустима доработка программы-функции на случай любого числа уравнений путем изменения числа аргументов программы-функции, а также числа строк в программе: в задании начальных условий и в цикле вычисления массива решений.
Аналогично методу Эйлера векторная и программная реализации вычислений по методу Рунге-Кутта могут быть распространены на случай решения системы дифференциальных уравнений. Прямое решение задачи в этом случае представляется достаточно громоздким и желательно предусмотреть расчета констант метода в рамках внутренних циклов. В пакете MathCAD имеются встроенные функции, решающие подобные задачи. Так, процедуру расчета приближенного решения по методу Рунге-Кутта решается с помощью функций rkfixed и Rkadapt. Пример решения той же системы дифференциальных уравнений дан в примере 4. Отличие в точности полученных решений для данного дифференциального уравнения невелико, однако, если решением дифференциального уравнения является сильно осциллирующая функция, то метод с переменным шагом обеспечивает большую точность.
Примеры решения различных уравнений с помощьюкомпьютерного пакета MathCAD.
Задание к контрольной роботе:
Дифференцируем левые части уравнения методом Лапласа, а затем суммируем их.
Заменяя функцию y(x) на Y, получим выражение, которое потом упрощаем с помощью функции collect.
Находим правую часть уравнения таким же способом
Получаем исходное алгебраическое уравнение которое решаем в символьном виде:
Решение уравнения
Теперь выполняем обратное преобразование и находим решение исходной дифференциальной задачи
Ответ
 . |
Проводим проверку правильности начальных условий, используя функцию subtitute
Графическое представление результатов решения:
|
