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

Суть метода состоит в проведении интегрального преобразования Лапласа функции, входящей в состав дифференциального уравнения, по правилу:

,

где s = + j - комплексная переменная величина.

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

Преобразование Лапласа можно выполнить, используя блок символьных вычислений MathCAD. Этот же блок позволяет выполнить и обратное преобразование Лапласа, в соответствии с соотношением:

,

где , т. е. интегрирование проводится по прямой, лежащей в плоскости комплексного переменного s и проходящей параллельно мнимой оси j на расстоянии от нее, при этом Лаплас образ Y(s) должен иметь особенности слева от этой линии.

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

Порядок решения дифференциального уравнения с помощью операционного исчисления представляется следующим:

- выполняя преобразование Лапласа левой и правой части дифференциального уравнения, учитываем начальные условия и переходим от дифференциального уравнения для функции оригинала y(t) к алгебраическому уравнению для Лаплас образа - Y(s) ;

- решая алгебраическое уравнение, находим в пространстве изображений в явном виде выражение для Y(s);

- выполняя обратное преобразование мы находим неизвестную функцию y(t).

Все этапы этой процедуры могут быть автоматизированы и выполнены в рамках пакета MathCAD (пример 1).

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

,

где s₁, s₂,…s_n - корни уравнения . В примере 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.

Пример №1

Пример №2

Пример №3

Пример №4

Задание к контрольной роботе:

Дифференцируем левые части уравнения методом Лапласа, а затем суммируем их.

Заменяя функцию y(x) на Y, получим выражение, которое потом упрощаем с помощью функции collect.

Находим правую часть уравнения таким же способом

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

Решение уравнения

Теперь выполняем обратное преобразование и находим решение исходной дифференциальной задачи

Ответ

.

Проводим проверку правильности начальных условий, используя функцию subtitute

Графическое представление результатов решения: