Методы и алгоритмы компьютерного решениядифференциальных уравнений
Введение
Для того, чтобы описатьдинамику различных процессов, протекающих в природных и в технических системах,составляют, опираясь на физические законы, дифференциальные уравнения. Так, вчастности, приходится поступать при исследовании функционированияавтоматических систем; работы судовых энергетических комплексов, электрическихагрегатов, судовых вспомогательных механизмов, систем навигации и т.д. В рядеслучаев эти уравнения допускают линеаризацию и могут быть записаны в виде:
/>,
где y(t) –неизвестная функция, a0, a1,...an– постоянные коэффициенты, а j(x) – некоторая известная функция независимого аргумента t,которая обычно выражает внешнее воздействие, оказываемое на систему.
1.Цель контрольной работы
Приобретение навыковалгоритмизации и программирования задач численного интегрирования обыкновенныхдифференциальных уравнений и систем с последующим моделированием результатов наперсональном компьютере и представлением их в виде таблиц и графиков.
В результате выполненияконтрольной работы студент обязан:
1. Научиться решать линейныедифференциальные уравнения численными и символьными методами в рамках пакетакомпьютерной математики MathCAD.
2. Ознакомиться с основными алгоритмамисуществующих компьютерных методов.
3. Определить точность этих методовпутем сравнения результатов, получаемых путем приближенного и аналитическогорешений.
2. Аналитические методы
Общее решениедифференциального уравнения n-го порядка – неизвестная функция y(t)– содержит n произвольных постоянных. Их можно определить, знаяначальные условия, накладываемые на неизвестную функцию и на ее производныевплоть до (n-1)-порядка включительно. Аналитически (в символьном виде) такиеуравнения решают классическим и операционным методами.
2.1 Классический метод
В ограниченном числеслучаев вида левой части (1) допускает такое преобразование, которое позволяетнайти решение путем непосредственного интегрирования, однако в общем случаепорядок решения – иной.
Решение неоднородногодифференциального уравнения (с ненулевой правой частью) является суммой общегорешения соответствующего однородного дифференциального уравнения y1(t)и частного решения y2(t) неоднородногодифференциального уравнения (1).
Решение однородногоуравнения ищем в виде: />. Подстановка егов дифференциальное уравнение приводит к характеристическому алгебраическомууравнению n-ного порядка:
/>,
которое имеет n корней – />. В частном случаеотсутствия кратных корней общее решение может быть записано в виде:
/>,
где Сi– произвольные постоянные, которые находятся из начальных условий.
Имеются правила,позволяющие определить вид y2(t) частного решения взависимости от вида правой части – функции j(t). Последующая подстановка общего решения в исходноедифференциальное уравнение позволяет найти неопределенные константы Ciв выражении для y1(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 «отказывается» работать. В этомслучае ему необходимо «помочь» врукопашную выполнив разложения полинома в знаменателев соответствии с соотношением:
/>,
где s1,s2,…sn – корни уравнения />. В примере 1 рассмотреновыполнение обратного преобразования Лапласа и для такого случая.
Рассмотренная методиканахождения аналитического решения дифференциальных уравнений может бытьраспространена на задачу решения системы дифференциальных уравнений. В этомслучае необходимо решить не одно алгебраическое уравнение для Лаплас-образов, асистему алгебраических уравнений с помощью той же процедуры блока решений Given– Find. Отметим, что в отличие от систем компьютерной математики Mathematica2.2.2 и Maple V R3/R4, которые легко позволяют аналитически решить линейноедифференциальное уравнение с помощью встроенных средств. Система MathCADпредполагает «ручные процедуры» запуска прямого преобразования Лапласа,составления по его результатам алгебраического уравнения и, после его решения,запуска процедуры обратного преобразования Лапласа.
3. Численный метод решения дифференциальных уравнений
Часто при анализе изучаемых втехнических или в природных системах процессов приходится учитывать наличиенелинейного поведения функций, описывающих характеристики их элементов. Это всвою очередь определяет появление нелинейностей в дифференциальных уравнениях,которые теперь уже не могут быть записаны в форме (1). Наличие нелинейностей вдифференциальных уравнениях обуславливает невозможность их точногоаналитического решения, а приближенные аналитические методы часто приводят кгромоздким выкладкам.
Кроме того, коэффициенты в левойчасти дифференциального уравнения могут быть определены из эксперимента сошибкой, что, в значительной степени, обесценивает получаемый точный аналитическийрезультат.
И наконец, точные методы не пригодныдля случая, если правая часть дифференциального уравнения представлена не ваналитической форме, а в виде таблицы или графика.
Во всех этих случаях прекраснопродолжают «работать» методы численного решения. В отличие от аналитических,они позволяют получать искомые зависимости для любой из описанных вышеситуаций. Алгоритмы существующих методов численного решения были разработанысравнительно давно, однако толчок к их применению был обусловлен развитиемвычислительной техники. Каждый из существующих численных методов предполагаетзамену производной на конечное приращение и преобразование дифференциальногоуравнения в уравнение в конечных разностях.
С этой целью интервал поиска решенияразбивается на множество отрезков и решение ищется на каждом из этих кусочков.Ясно, что чем мельче шаг разбиения, тем точнее получается результат. Поэтому,эффективное применение численных методов (при решении реальных, а не учебныхзадач!) предполагает использование компьютеров с достаточным быстродействием.
Использование для численного решениядифференциальных уравнений компьютерного пакета MathCAD предполагает знаниеалгоритма работы численных методов для разумного их применения (знание границприменимости, оценки точности, затрат компьютерных ресурсов и др.). Дело втом, что к результатам компьютерных вычислений всегда нужно относиться критически;анализировать их на правдоподобность, и для того, чтобы избежать«подводных камней» при использовании любого стандартного пакета,реализующего численные методы, нужно иметь хотя бы минимальное представление отом, какой именно численный метод реализован для решения той или иной задачи икак он «работает».
Поскольку для оценки точности решениянеобходим материал для сравнения предлагается рассматривать приближенные(численные) методы решения применительно к тем задачам, для которых ранее былополучено аналитическое решение (т.е. линейных уравнений вида (1)).
Обозначим />/>. Тогда данное уравнениеможно преобразовать в следующую систему уравнений первого порядка, разрешенныхотносительно первой производной (форма Коши):
/>
Поэтому, идея применения численныхметодов для решения уравнения старших порядков, в принципе, ничем не отличаетсяот идеи численного решения уравнения первого порядка, которые рассматривались влабораторной работе №10.
Нетрудно расширитьприменение описанной выше методики на случай системы линейных уравнений. Впримерах 2 и 3 представлены реализации метода Эйлера в рамках векторнойпроцедуры и с помощью программы-функции. В обоих случаях получены одинаковыерезультаты, которые поверяются по полученному ранее (см. пример 1)аналитическому решению. Решение приближенное и точное сильно отличаются и величинаотносительной ошибки (для выбранного шага) составляет ~13%. Пример 3 нетруднооптимизировать и придать ему более компактный вид, считая начальные условия иправые части системы уравнений компонентами некоторых векторов. Также допустимадоработка программы-функции на случай любого числа уравнений путем изменениячисла аргументов программы-функции, а также числа строк в программе: в заданииначальных условий и в цикле вычисления массива решений.
Аналогично методу Эйлеравекторная и программная реализации вычислений по методу Рунге-Кутта могут бытьраспространены на случай решения системы дифференциальных уравнений. Прямоерешение задачи в этом случае представляется достаточно громоздким и желательнопредусмотреть расчета констант метода в рамках внутренних циклов. В пакетеMathCAD имеются встроенные функции, решающие подобные задачи. Так, процедурурасчета приближенного решения по методу Рунге-Кутта решается с помощью функцийrkfixed и Rkadapt. Пример решения той же системы дифференциальных уравнений данв примере 4. Отличие в точности полученных решений для данногодифференциального уравнения невелико, однако, если решением дифференциальногоуравнения является сильно осциллирующая функция, то метод с переменным шагомобеспечивает большую точность.
Примеры решения различныхуравнений с помощью компьютерного пакета MathCAD.Пример №1
/>
/>
/>
/> Пример №2
/>
/>
/>
/> Пример №3
/>
/>
/>
/>
/>
/> Пример №4
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
Задание к контрольнойроботе:
/>
/>
/>
Дифференцируем левыечасти уравнения методом Лапласа, а затем суммируем их.
/>
/>
/>
Заменяя функцию y(x) на Y, получимвыражение, которое потом упрощаем с помощью функции collect.
/>
Находим правую часть уравнения такимже способом
/>
Получаем исходное алгебраическоеуравнение которое решаем в символьном виде:
/>
/>
Решение уравнения
/>
Теперь выполняем обратноепреобразование и находим решение исходной дифференциальной задачи
/>
Ответ
/>
/> .
Проводим проверку правильностиначальных условий, используя функцию subtitute
/>
/>
Графическоепредставление результатов решения:
/>