Как решать дифференциальные уравнения в маткаде

MathCAD — это просто! Часть 21. Продолжаем борьбу с дифференциальными уравнениями

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

Мы с вами уже говорили когда-то о символьном процессоре MathCAD, и я рассказывал вам, что он, мягко говоря, далеко не всемогущ, и уж точно не настолько всемогущ, чтобы решать абсолютно все, что ему подсунет пользователь. В общем-то, причина этого вполне понятна: решение многих математических задач требует творческого подхода — откровенно говоря, лишь в очень малом числе задач оно лежит достаточно близко к поверхности, чтобы его сходу мог найти математический процессор MathCAD’а. Хотя алгоритмы, лежащие в его основе, чрезвычайно сложны (в отличие, кстати, от тех, на которых построены численные способы решения), они вовсе не безупречны и даже близко не приближают компьютер по мыслительным и аналитическим способностям к человеку. Когда речь идет о дифференциальных уравнениях, все сказанное не просто верно, а верно в квадрате. Дифференциальные уравнения сами по себе довольно сложны и коварны, и потому не зря, как я говорил в прошлый раз, вызывают трепет у студентов, имеющих «удовольствие» их изучать. Конечно, существуют выработанные многими поколениями математиков способы решения основных типов дифференциальных уравнений, однако в реальных задачах (особенно физических и технических) иногда возникают уравнения, которые сложно решить с помощью классических формул. И символьный процессор MathCAD’а тоже опускает руки. Как быть в этом случае? Сдаться и признать такое уравнение нерешаемым? Думаю, вряд ли это такая уж хорошая идея. Гораздо эффективнее вооружиться тем инструментом, который придуман специально для решения того, что невозможно решить аналитически — численными методами. Изначально MathCAD проектировался и создавался именно как средство численного решения математических задач, а потому это удается ему особенно хорошо. Решать дифференциальные уравнения в этой мощной математической среде численно ничуть не сложнее, чем аналитически, и сейчас вы сможете в этом воочию убедиться.

Небольшое отступление

Чтобы лучше понимать суть численного решения дифференциальных уравнений, стоит сначала напомнить, что именно представляют собой решения дифференциальных уравнений. Для простоты остановимся на простейшем (простите за тавтологию, но здесь она неизбежна) случае — обыкновенном дифференциальном уравнении с одной переменной первого порядка (т.е. только с первыми производными). Давайте вспомним, что именно мы ищем, решая дифференциальное уравнение? Ищем мы некоторую функцию, которая будет удовлетворять данному уравнению — то есть при подстановке в него превращать его в тождество (верное числовое равенство). Такая функция может существовать, а может и нет. Будем считать, что она существует — ведь в противном случае получается, что дифференциальное уравнение, которое мы пытаемся решить, на самом-то деле нерешаемо. Как мы помним, любое дифференциальное уравнение должно включать в себя хоть одну производную искомой величины — иначе оно не будет дифференциальным, и решать его нужно будет с помощью методов, применяемых к обычным уравнениям. Но, поскольку уравнение содержит производные — значит, его решение верно с точностью до константы. То есть, вообще говоря, решений бесконечное множество, и каждое отличается от другого каким-то постоянным членом. Ведь при дифференцировании любая постоянная обращается в ноль, и, значит, решение с любой константой будет верным.

Если интерпретировать геометрически решения дифференциального уравнения, мы получим семейство кривых, проходящих через все точки плоскости. Для того, чтобы выделить из этих всех кривых какую-то одну, на уравнение накладывают дополнительные ограничения. Формулироваться эти ограничения могут либо в виде задачи Коши, либо в виде граничной задачи. Задача Коши определяет заданную кривую с помощью задания значений самой искомой функции и ее производных (для уравнений высших порядков) в какой-то определенной точке. Граничная же задача формулируется как задание значений функции на границах искомого отрезка. Физическая интерпретация этих двух ограничений такова: задача Коши определяет решение временной задачи, задавая значения переменных в начальный момент времени. Граничная же задача (она также часто называется краевой) нужна для решения пространственных задач: она определяет начальные значения переменных в рамках какой-то пространственной области. Решение этих двух типов задач численно в MathCAD’е отличается, и вы сами сможете убедиться в том, что отличия эти порой довольно-таки существенны. Поэтому мы пойдем, как водится, от простого к сложному и начнем решение дифференциальных уравнений с рассмотрения задачи Коши.

Задача Коши для обыкновенных дифференциальных уравнений

Давайте посмотрим, как можно решать задачу Коши для обыкновенного дифференциального уравнения. В качестве уравнения возьмем простое до ужаса dX(t)/dt + X(t) = 0. Думаю, вы сможете сразу сказать, что за выражение будет решением данного уравнения. Конечно же, это экспонента с показателем минус t. Тем не менее, для того, чтобы удостовериться, что MathCAD правильно решает дифференциальные уравнения численными методами, мы подсунем это простейшее уравнение ему и посмотрим, как он с ним справится. Для того, чтобы показать MathCAD’у, где начинается условие нашей задачи Коши, мы должны воспользоваться ключевым словом MathCAD’а — Given (обязательно писать его с большой буквы, поскольку, как вы помните, среда MathCAD чувствительна к регистру букв в именах переменных, функций и ключевых слов). После того, как записано это слово, мы должны записать начальные условие. Поскольку у нас уравнение первого порядка, то оно у нас будет всего одно — например, можно написать X(0) = 1. Далее (то есть под начальным условием) мы должны записать само дифференциальное уравнение, решением которого мы нагрузим MathCAD. Ну, а для поиска решений мы применим функцию Odesolve. У функции этой предусмотрено два (в нашем случае) параметра. Первый — это имя переменной, по которой будет проводиться интегрирование уравнения (в нашем случае, сами понимаете, выбирать, собственно говоря, не из чего — это будет t), а второй параметр — конечное значение интегрирования, задающее верхнюю границу того интервала, для которого MathCAD будет искать решение нашего уравнения. Собственно, после функции Odesolve уже можно ничего и не писать — если уравнение имеет решения, то оно обязательно будет решено.

Не писать-то, конечно, можно, только это не лучший вариант. Как проверить решение уравнения? Исходя из того, что у нас сейчас простейший одномерный случай, проще всего построить график этого решения. По нему отлично видно, что MathCAD самым что ни на есть замечательным образом справился с поставленной перед ним задачей и решил уравнение именно так, как оно и должно решаться: перед нами — не что иное, как график экспоненциальной функции с отрицательным показателем.

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

Краевая задача для обыкновенных дифференциальных уравнений

Давайте теперь еще коротко поговорим о решении граничной задачи для обыкновенных дифференциальных уравнений. С помощью функции Odesolve, в общем-то, решать ее нужно будет точно так же, как и задачу Коши — с той только, конечно, разницей, что вместо начального условия нужно записать граничные. В качестве простой иллюстрации рассмотрим дифференциальное уравнение d2X(t)/dt2 + dX(t)/dt + X(t) = 0. Граничные условия будут такими: X(0) = 10, X(10) = 0. Конечно, теперь уже такую простую аналитическую формулу решения, как для первого уравнения, нам подобрать не удастся, однако в этом нет ничего страшного, потому что MathCAD сможет ничуть не хуже построить график решения нашего обыкновенного дифференциального уравнения.

Конечно, использование функции Odesolve — не единственный метод решения обыкновенных дифференциальных уравнений в среде MathCAD, однако он достаточно удобен в силу своей простоты. В следующий раз мы с вами, конечно, познакомимся и с другими способами их решения, которые предлагает пользователю эту мощная математическая среда, но пока что достаточно — нехорошо перегружать читателя информацией, особенно по такому сложному предмету, как дифференциальные уравнения. Пусть даже MathCAD и облегчает их решение, все равно лучше не торопиться с тем, чтобы в них въехать.

Подведем итоги

Что ж, давайте подытожим все то, о чем мы с вами успели сейчас поговорить. Как обычно, итоги будут простыми и жизнеутверждающими. Численное решение дифференциальных уравнений — гораздо более универсальный способ их решения, чем аналитическое решение, поскольку численные методы как раз и создавались для тех случаев, когда искать аналитическое решение не слишком удобно или попросту абсурдно. Правда, конечно, в результате численного решения мы не получим простую и красивую формулу (хотя, когда речь идет о дифференциальных уравнениях, весьма наивно ожидать ее получение и в результате аналитического решения). Зато мы получим хоть какой-то результат — а это, согласитесь, гораздо лучше, чем отсутствие вообще любого результата. То есть, как говорит мудрая пословица, лучше синица в руках, чем журавль в небе. Это как раз таки именно тот самый случай. Решать численно дифференциальные уравнения в MathCAD, как вы могли убедиться, совсем не сложно, однако, естественно, везде есть свои тонкости, и численные методы совсем не панацея. Немного об этом рассказано в справке самого MathCAD’а, и именно ее чтением вы можете заняться, чтобы не терять время зря в ожидании следующей статьи из серии «MathCAD — это просто».

SF, spaceflyer@tut.by

Компьютерная газета. Статья была опубликована в номере 35 за 2008 год в рубрике soft

Как решать дифференциальные уравнения в маткаде

БлогNot. MathCAD: решаем основные типы дифференциальных уравнений встроенными функциями

MathCAD: решаем основные типы дифференциальных уравнений встроенными функциями

Решать дифференциальные уравнения (далее ДУ) в MathCAD, составляя собственные подпрограммы-функции не всегда удобно и экономично по времени, хотя и полезно на этапе обучения. Опишем в этой заметке способы решения основных типов ДУ с помощью стандартных средств пакета, ограничимся простыми примерами.

1. ДУ с разделяющимися переменными. Общая постановка задачи: y’=f(x,y)=g(x)*h(y) , y(x0)=y0 . То есть, f(x,y) допускает представление в виде произведения функций от x и от y .

Для решения уравнения достаточно задать его правую часть как пользовательскую функцию MathCAD, определить интервал поиска решения [x0,x1] , начальное условие y0 и применить стандартную функцию Odesolve . Покажем этот процесс на примере уравнения y’=2x-y+x 2 , x∈[0,2] , y(0)=0 с известным решением y(x)=x 2 :

Решение дифференциальных уравнений с разделяющимися переменными

Знак «равно» в записи уравнений, конечно же, жирный (панель Boolean или сочетание клавиш Ctrl+=).

Функция Odesolve вернула именно функцию y , её нужно смотреть от аргумента, например, y(1)= .

И ещё 2 особенности:

Зная точное решение, графически сравним с ним найденное решение. Как видно на графике, MathCAD справился с задачей отлично.

Графики полученного и точного решения совпадают

2. Неоднородное ДУ первого порядка. В общем виде такое уравнение можно записать как y’=a(x)*y+b(x) . Оно решается аналитически по формуле, которую можно найти в любой книге по решению обыкновенных ДУ:

Формула для решения неоднородного ДУ первого порядка

Здесь С – константа интегрирования. Остаётся применить формулу к конкретному уравнению (возьмём для примера задачу y’+2xy=x*e -x 2 sin(x) , y(0)=1 ) и оценить её символьно:

Аналитическое решение неоднородного ДУ первого порядка в MathCAD

Здесь мы получаем решение в общем виде. Нижний оператор оценён символьно (см. панель «Символика»), а аргумент t используется, так как x в документе «уже занят» (для корректной работы символьной оценки переменные не должны быть определены заранее).

После подстановки начального условия получим частное решение y(t,Y₀) , а для проверки решения будет достаточно подставить полученную функцию в исходное уравнение и упростить его символьной функцией simplify . Полученный результат в нашем случае совпал с заданной в условии правой частью. Также для y(x,Y₀) , как и для любой функции, можно построить график на нужном интервале изменения x.

Проверка решения неоднородного ДУ первого порядка и построение графика

3. Неоднородное ДУ второго порядка. В общем виде имеем уравнение y» + p(x)*y’ + g(x)*y = f(x) плюс набор краевых условий, количество которых соответствует порядку задачи, например y(0)=. y'(0)=. или y(0)=. y(1)=.

Возьмём уравнение, которое мы мучили вот здесь, решим его стандартными средствами, сравним с известным точным решением и построим график:

Решение неоднородного ДУ второго порядка в MathCAD

Здесь при вызове Odesolve второй параметр, равный единице — это правая граница интервала, третий параметр, равный 10, задаёт количество интервалов. Точное решение u(t) взяли по ссылке. Как видим, даже на 10 интервалах всё очень хорошо совпадает.

4. Система ДУ. Подход к решению системы ДУ покажем на примере. Пусть задана система дифференциальных уравнений

x’ = a*x — y — (x 2 + y 2 )*x,
y’ = a*y + x — (x 2 + y 2 )*y,
x(0)=0, y(0)=1, a=-0.2

Чтобы решить эту систему стандартной функцией rkfixed , нужно задать для неё вектор начальных значений x = (x₀, y₀) и вектор правых частей D(t,x) .

После этого задача решится вызовом rkfixed , второй и третий параметры ( 0, 20 ) задают интервал по времени t , на котором ищется решение, четвёртый параметр 100 означает количество точек на интервале.

Функция вернёт матрицу решений системы, в которой количество строк соответствует количеству точек на интервале, а количество столбцов — количеству уравнений в системе.

Для построения графика достаточно отобразить зависимость столбцов Z_i,1 , Z_i,2 от Z_i,0 , i=0..99 :

Решение системы ДУ в MathCAD функцией rkfixed

Решение дифференциальных уравнений в MathCad 15

В свободном поле mathcad введите оператор Given. Этот оператор запускает процесс ввода исходных данных для корректной работы функции odesolve. После этого найдите панель под названием Calculus. В этой панели нам понадобятся кнопки Derivative и Nth Derivative. Эти кнопки вводят заготовки для дифференциального уравнения. С помощью клавиатуры введите уравнение, как показано на рисунке 1. Знак равенства необходимо использовать из панели Boolean

Рис. 1. Ввод исходных данных для решения дифференциального уравнения

Далее введите начальные приближения. Количество начальных приближений зависит от порядка дифференциального уравнения. В нашем случае их будет 2. Введите приближения, как показано на рисунке 2. Обратите внимание, что для ввода значания первой производной вам нужно использовать символ «верхний апостроф«. Если вы его не можете ввести с клавиатуры вручную воспользуйтесь приложением windows Capter Map или используйте комбинацию клавиш Alt + 96 или Alt + 39

Рис. 2. Ввод первого приближения для решения дифференциального уравнения

Теперь, после начального приближения введите любую переменную (например y) и присвойте ей функцию Odesolve, как показано на рисунке 3. В качестве параметров функции Odesolve используется переменная t и интервал интегрирования. В нашем случае интервал равен 15

Рис. 3. Ввод функции odesolve для решения дифференциального уравнения

Можно отобразить функцию y на графике, где в качестве аргумента будет переменная t. Этот график и будет являться решением дифференциального уравнения. Обратите внимание, что график строится в пределах интервала интегрирования. Особенности оформления и отображения графиков в mathcad 15 смотрите в соответствующем разделе

Рис. 4. Вывод результата решения дифференциального уравнения на график

После корректного решения дифференциального уравнения функцию y(t) можно использовать далее в расчетах

Рис. 5. Результат решения дифференциального уравнения в mathcad 15 и старше

Donec eget ex magna. Interdum et malesuada fames ac ante ipsum primis in faucibus. Pellentesque venenatis dolor imperdiet dolor mattis sagittis. Praesent rutrum sem diam, vitae egestas enim auctor sit amet. Pellentesque leo mauris, consectetur id ipsum sit amet, fergiat. Pellentesque in mi eu massa lacinia malesuada et a elit. Donec urna ex, lacinia in purus ac, pretium pulvinar mauris. Curabitur sapien risus, commodo eget turpis at, elementum convallis elit. Pellentesque enim turpis, hendrerit tristique.

Lorem ipsum dolor sit amet, consectetur adipiscing elit. Duis dapibus rutrum facilisis. Class aptent taciti sociosqu ad litora torquent per conubia nostra, per inceptos himenaeos. Etiam tristique libero eu nibh porttitor fermentum. Nullam venenatis erat id vehicula viverra. Nunc ultrices eros ut ultricies condimentum. Mauris risus lacus, blandit sit amet venenatis non, bibendum vitae dolor. Nunc lorem mauris, fringilla in aliquam at, euismod in lectus. Pellentesque habitant morbi tristique senectus et netus et malesuada fames ac turpis egestas. In non lorem sit amet elit placerat maximus. Pellentesque aliquam maximus risus, vel venenatis mauris vehicula hendrerit.

Interdum et malesuada fames ac ante ipsum primis in faucibus. Pellentesque venenatis dolor imperdiet dolor mattis sagittis. Praesent rutrum sem diam, vitae egestas enim auctor sit amet. Pellentesque leo mauris, consectetur id ipsum sit amet, fersapien risus, commodo eget turpis at, elementum convallis elit. Pellentesque enim turpis, hendrerit tristique lorem ipsum dolor.

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

Mathcad предлагает новый способ для решения обыкновенных дифференциальных уравнений, разрешенных относительно старшей производной. Для этих целей служит уже известный нам блок given совместно с функцией odesolve. Дифференциальное уравнение совместно с начальными или граничными условиями записывается в блоке given. Производные можно обозначать как штрихами (Ctrl+F7), так и с помощью знака производной . Пример использования функции для решения задачи Коши приведен ниже.

У искомой функции явно указан аргумент, знак производной стоит перед скобкой.

Функция odesolve имеет три аргумента. Первый аргумент — независимая переменная, вторая — граница интервала, на котором ищется решение, последний аргумент — шаг, с которым ищется решение. Последний аргумент может быть опущен.

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

Решение уравнений в частных производных

Одним из методов решения дифференциальных уравнений в частных производных является метод сеток. Идея метода заключается в следующем. Для простоты, ограничимся случаем только функции двух переменных, и будем полагать, что решение уравнения ищется на квадратной области единичного размера. Разобьем область сеткой. Шаг сетки по оси x и по оси y, вообще говоря, может быть разный. По определению частная производная равна

Если рассматривать функцию только в узлах сетки, то частную производную можно записать в форме

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

Такое выражение называется левой конечной разностью. Можно получить центральную конечную разность, найдя среднее этих выражений.

Теперь получим выражения для вторых производных.

В данном случае для нахождения производной мы использовали симметричные точки. Однако, очевидно, можно было бы использовать точки с несимметричным расположением.

1. Уравнения гиперболического типа

В качестве примера рассмотрим решение волнового уравнения (уравнения гиперболического типа).

Уравнение будем решать методом сеток. Запишем уравнение в конечных разностях

Полученное уравнение позволяет выразить значение функции u в момент времени через значения функции в предыдущие моменты времени.

Такая разностная схема называется явной, так как искомая величина получается в явном виде. Она устойчива, если.

Зададим начальные условия: смещение струны U в начальный и последующий моменты времени описывается синусоидальной функцией.

(Совпадение смещений при j=0 и j=1 соответствует нулевой начальной скорости).

Зададим граничные условия: на концах струны смещение равно 0 в любой момент времени

Будем полагать коэффициент

Записываем уравнение в конечных разностях, разрешенное относительно

Представляем результат на графике

2. Уравнения параболического типа

Еще один пример использования конечных разностей — уравнение диффузии.

Это уравнение параболического типа. Явная разностная схема для этого уравнения имеет вид

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

и диапазон изменения пространственной и временной координат:

Задаем начальные и граничные условия

Уравнение в конечных разностях имеет вид

Представляем результаты на графике. (Для большей наглядности изображена только центральная часть)

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

Рассмотрим применение неявной схемы на примере уравнения теплопроводности

Запишем неявную разностную схему для этого уравнения

Здесь первый индекс соответствует пространственной, а второй — временной координате. В отличие от явной схемы, для вычисления в правой части уравнения используются значения функции на том же самом временном шаге. Вводя обозначение , уравнение (5) можно переписать в виде

или в матричной форме

Задаем количество узлов сетки (в данном случае оно одинаково для обеих переменных)

Задаем значения параметров

и начальное распределение температуры в области

Формируем матрицы уравнения (7)

Находим решение системы

3. Решение уравнений Лапласа и Пуассона

Для решения уравнений Пуассона и Лапласа (частный случай, когда ) — уравнений эллиптического типа — предназначена функция relax (a, b, c, d, e, f, u, rjac), реализующая метод релаксации. Фактически, эту функцию можно использовать для решения эллиптического уравнения общего вида

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

В частности, для уравнения Пуассона коэффициенты .

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

При наличии источников разностная схема имеет вид

Метод релаксации сходится достаточно медленно, так как фактически он использует разностную схему (3) с максимально возможным для двумерного случая шагом .

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

Функция relax возвращает квадратную матрицу, в которой:

расположение элемента в матрице соответствует его положению внутри квадратной области,

это значение приближает решение в этой точке.

Эта функция использует метод релаксации для приближения к решению.

Вы должны использовать функцию relax, если Вы знаете значения искомой функции u (x, y) на всех четырех сторонах квадратной области.

a, b, c, d, e — квадратные матрицы одного и того же размера, содержащие коэффициенты дифференциального уравнения.

f — квадратная матрица, содержащая значения правой части уравнения в каждой точке внутри квадрата

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

rjac — Параметр, управляющий сходимостью процесса релаксации. Он может быть в диапазоне от 0 до 1, но оптимальное значение зависит от деталей задачи.

Задаем правую часть уравнения Пуассона — два точечных источника

Задаем значения параметров функции relax

Задаем граничные условия и начальное приближение — нули во всех внутренних точках области