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

Сразу следует отметить, что с обыкновенными дифференциальными уравнениями и системами этих уравнений Maple справляется достаточно неплохо. Если уравнение в принципе решается, то Maple, скорее всего, его решит. Полезной в этом случае будет процедура dsolve(), параметрами которой указываются уравнение (система уравнений), начальные условия (если такие имеются), а также функция (или набор функций для системы уравнений), относительно которой это уравнение (систему) следует решать. Рассмотрим пример.

Задача 5.1

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

В первую очередь определяем уравнение, присвоив его в качестве значения переменной Eq (запись D@@2 означает, что оператор D действует дважды).

На заметку
Уравнение есть прямым следствием второго закона Ньютона. Если координатную ось выбрать так, чтобы она была направлена вниз, т.е. в направлении к поверхности Земли, а начало отсчета совместить сточкой начала движения парашютиста, то зависимость координаты парашютиста от времени будет тем самым определять расстояние, которое парашютист пролетел. Ускорение равно, как известно, второй производной от координаты по времени. Кроме того, на парашютиста действуют две силы: сила тяжести, равная произведению массы парашютиста на ускорение свободного падения (д) — она направлена вниз, ее направление совпадает с направлением координатной оси и поэтому проекция данной силы на координатную ось положительна, и сила сопротивления воздуха — направлена вверх и поэтому ее проекция отрицательна. После сокращения формулы на массу парашютиста из второго | закона Ньютона получаем нужное уравнение (а — коэффициент пропорциональности в выражении для силы сопротивления воздуха, нормированный на массу парашютиста).

Предполагаем, что парашютист начинает движение без начальной скорости и в начальный момент координата парашютиста равна нулю. Начальные условия записываем в виде последовательности равенств в переменную In_Con.

После этого уравнение можно решать.

Первым аргументом процедуры dsolve() указано множество, состоящее , из уравнения и начальных условий. Уравнение решается относительно функции x(t).

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

Полученное решение, тем не менее, является достаточно громоздким. Его можно несколько упростить.

Чтобы задать функциональную зависимость x(t) согласно полученному выше равенству и не набирать его при этом непосредственно с клавиатуры, Воспользуемся процедурой unapply(). Эта процедура позволяет выделить функциональную зависимость, так что не придется прибегать даже к копированию выражения в буфер обмена.

Первым параметром процедуры unapply() указывается выражение, из которого "извлекается" функциональная зависимость, а вторым параметром — переменная, по отношению к которой определяется данная функциональная зависимость. Другими словами, результатом выполнения процедуры unapply() есть оператор, действие которого на указанную вторым параметром переменную продуцирует выражение, указанное первым параметром.

Оператор v определяем как результат действия оператора дифференцирования D на функциональный оператор х (в результате зависимость v(t) будет определять скорость парашютиста в момент времени t).

Скорость парашютиста в таком случае дается следующей закономерностью.

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

Далее исследуем динамику парашютиста по прошествии существенных промежутков времени. Вычислительному ядру Maple сообщим, что параметры д и а положительны.

После этого вычисляем предел, к которому стремится скорость при стремлении времени к бесконечности.

Таким образом, в пределе бесконечно больших времен скорость парашютиста выходит на стационарное значение — движение будет равномерным.

Еще один интересный пример — математический маятник с трением. Об этом в следующей задаче.

Задача 5.2

Составить и решить дифференциальное уравнение для математического маятника с трением.

В отличие от уравнения для маятника без трения, в данном случае будет присутствовать слагаемое, пропорциональное скорости маятника, т.е. первой : производной. Для удобства коэффициент пропорциональности обозначим через a Таким образом, рассматриваем следующее уравнение.

Данное уравнение можно решить в общем виде; в этом случае начальные , условия не указываются.

Переменные среды _С1 и _С2 являются, с математической точки зрения, Произвольными константами, которые в общем случае определяются из дополнительных условий — как правило, начальных.

Далее с помощью уже известной процедуры unapply() определяем соответствующую функциональную зависимость.

Если предположить, что в начальный момент отклонение маятника от положения равновесия равнялось А, а начальная скорость была равна нулю, то константы С1 и С2 можно определить как решение системы уравнений (обычных, не дифференциальных!).

Чтобы переменным _С1 и _С2 присвоить соответствующие значения (при решении такие значения находятся, но не присваиваются), можно воспользоваться процедурой assign(). Если в качестве аргумента этой процедуры указано равенство, то, согласно этому равенству, вычислительным ядром Maple будет выполнено присваивание (формально можно полагать, что знак равенства заменяется оператором присваивания). Ниже в качестве аргумента указана переменная среды — ссылка на результат выполнения последней команды. Это множество, элементами которого являются два равенства, определяющие переменные _С1 и _С2. В результате выполнения приведенной ниже команды переменные получат значения.

Наличие в выражении экспонент смущать не должно: если параметр a меньше частоты в (это соответствует затухающим колебаниям), то показатели в экспонентах становятся комплексными и решения выражаются, на самом деле, через синус и косинус, правда с экспоненциальными множителями.

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

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

Далее строим график (если точнее, то поверхность) зависимости координаты маятника от времени и параметра а, определяющего трение в системе.

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

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

Задача 5.3

Найти решение задачи Коши.

В данном случае укажем уравнение, равно как и начальное условие, непосредственно аргументом процедуры dsolve().

Как можно видеть выше, получаем решение.

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

Задача 5.4

Найти общее решение уравнения y(x)+exp(-x)y=sin(x).

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

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

Задача 5.5

Найти решение задачи Коши

Задаем непосредственно уравнение.

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

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

Задача 5.6

Найти общее решение системы уравнений

Однако не всегда удается получить точное решение. Иногда приходится довольствоваться и приближенным. Именно об этом далее пойдет речь.

Приближенные методы решения дифференциальных уравнений

Среди приближенных методов решения дифференциальных уравнений достаточно распространенным является метод разложения по малому параметру. Идея, положенная в основу метода, проста: в уравнении (или системе) выделяется малый параметр, а решение ищется в виде ряда по этому параметру.

На заметку
К сожалению, далеко не в каждом уравнении такой малый параметр можно выделить.

В качестве примера рассмотрим следующую задачу.

Задача 5.7

Найти приближенное решение в виде многочлена второго порядка по малому параметру для задачи Коши.

Задаем исходное уравнение.

Поскольку это уравнение первого порядка, начальное условие только одно.

Строго говоря, данное уравнение вычислительным ядром Maple решается точно. Ниже приведена соответствующая команда, однако, без указания результата. Причина проста -— результат этот весьма нетривиален. Желающие могут убедиться в этом самостоятельно. Последнее, кстати, является свидетельством того, что точный результат искать не всегда полезно; зачастую достаточно ограничиться приближенным решением — оно может оказаться вполне приемлемым по точности и в то же время простым и наглядным.

Решение ищем в виде ряда по малому параметру — в данном случае это е. Поэтому у(х) представляем в следующем виде.

Внимание!
Представленная выше команда, с помощью которой функции у(х) присваивается значение, на самом деле функцию не определяет. Если в командной строке ввести команду у(х), в области вывода появится уО(х) + еу\(х) + егу2(х). Но если команду заменить, скажем, на y(t), результат будет y(t). Другими словами, присваивание выполняется "на уровне названий", а у(х) в данном случае является не результатом действия оператора у() на переменную х, а названием переменной.

Уравнение, таким образом, будет иметь следующий вид.

Преобразуем это уравнение, перенеся все слагаемые, содержащие малый параметр, в левую часть, а слагаемые, не содержащие малый параметр, — в правую. Сделать это можно с помощью процедуры isolate!). При вызове процедуры вычислительным ядром Maple предпринимается попытка в равенстве или алгебраическом выражении, указанном первым параметром, выразить подвыражение, указанное вторым параметром, т.е. фактически указанное первым параметром выражение решается относительно второго параметра-выражения. Таким образом, получаем следующее.

Теперь можно собрать слагаемые при соответствующих степенях малого параметра (epsilon).

Неизвестные функции yO(x), yl(x) и у2(х) выбираются так, чтобы коэффициенты при степенях малого параметра (не всех, а только до второго включительно, до которого раскладывается в ряд искомая функция) равнялись нулю.

Отсюда получаем первое уравнение (для определения уО(х)), приравняв правую часть уравнения Eq (это коэффициент при нулевой степени параметра epsilon) к нулю.

Приравнивая к нулю коэффициент при первой степени малого параметра, получаем еще одно уравнение.

В этом уравнении, помимо yl(x), присутствует и уО(х). Поэтому, прежде чем решать уравнение Eq_l, необходимо сначала решить уравнение Eq 0, найти тем самым уО(х), после чего можно будет решать Eq_l относительно yl(x).

На заметку
Коэффициент при первой степени параметра epsilon определяется процедурой coef f (). Первым ее аргументом указана левая часть уравнения Eq (lhs(Eq)) — именно в левой части содержатся слагаемые с epsilon. Второй параметр (epsilon) является указанием на то, что нужно выделить коэффициент при первой степени epsilon. Приравняв этот коэффициент к нулю, получаем нужное уравнение.

Практически точно так же получаем последнее дифференциальное уравнение — приравниваем к нулю коэффициент при второй степени epsilon. Отличие от предыдущего случая состоит в том, что вторым параметром в процедуре coeff() указано epsilonA2.

Решать уравнения Eq_0, Eq_l и Eq 2 следует строго в той последовательности, в какой они вводились. Причина очевидна — каждое последующее уравнение содержит, помимо неизвестной функции, еще и функции, которые определяются предыдущими уравнениями.

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

Итак, получаем следующее.

Ниже, согласно полученному решению, выполняется присваивание.

Внимание!
В этом месте функция уО(х) "технически" становится переменной!

Ситуация такая же, как и с функцией у(х). Например, если ввести команду уО(х), получим ожидаемый результат.

Однако если в приведенном выше вызове изменить аргумент (yO(t)), нужного значения (1/t) не получим.

Решаем теперь уравнение Eg 1 и находим yl(x).

В данном случае, перед присваиванием, раскладываем выражение для yl(x) на сумму дробей (процедура expand(%)).

Наконец, находим у2(х). Последовательность действий такая же.

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

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

Задача 5.8

Найти приближенное решение в виде многочлена второго порядка по малому параметру для задачи Коши 1 + х + еу, у(0) = sin(s).

Как и прежде, в первую очередь задаем само уравнение.

Начальные условия в этом случае определяются через малый параметр.

Решать уравнение будем таким образом, чтобы, при необходимости, можно было легко изменить порядок разложения (в условии сказано, что искать решение следует в виде многочлена второго порядка по малому параметру). Для этого вводим переменную N, которая и будет определять максимальный порядок при разложении решения в ряд по малому параметру. Инициализируем эту переменную, присвоив ей значение 2 (именно это требуется в условии).

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

Слагаемые добавлять будем с помощью оператора цикла.

В рамках этого цикла на каждом очередном шаге значение у(х) увеличивается на слагаемое, равное произведению параметра epsilon в соответствующей степени на название, формируемое объединением символа "у" и индекса — объединяется такая конбтрукция процедурой cat(). В результате получаем следующее.

Исходное уравнение при этом будет иметь такой вид.

Преобразуем это уравнение, перенеся все слагаемые в левую часть.

Далее левую часть уравнения Eq раскладываем в степенной ряд по параметру epsilon в окрестности нуля (команда series(lhs(Eq),epsilon=0,N+l)). Поскольку нас интересуют степени epsilon вплоть до N, в процедуре series () третьим параметром указано N+1 (на единицу больше) — этот параметр определяет порядок остатка, а сам ряд имеет порядок на единицу меньше. Разложение преобразуем в полиномиальный вид (процедура convert()). Чтобы уравнение осталось уравнением, полученное выражение приравниваем к нулю.

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

Названия этих переменных являются объединением имени "Eq" и индекса i, который изменяется в диапазоне от 0 до N.

Коэффициент при i-й степени epsilon определяется процедурой coeff(). Степень epsilon указывается третьим параметром процедуры.

На следующем этапе выполняем разложение в ряд по малому параметру начальных условий. В данном случае начальное условие одно, и записано оно в переменную inCon. Эта переменная — равенство. Поэтому переносим все слагаемые в одну сторону — это делается командой (lhs(InCon)-rhs(InCon)) — и полученное выражение раскладывается в ряд по epsilon в окрестности нуля.

Теперь выражение трансформируется в полином.

В этом выражении присутствует функция у(0) (ее значение в точке 0), которая ищется в виде ряда. Поэтому следует выполнить соответствующую замену. Для этого введем переменную s, которой присвоим первоначальное значение у0(0).

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

Теперь выполним замену.

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

Решаем первое уравнение (приняв во внимание начальные условия).

Описываем процедуру spDiffSol(), которая будет иметь пять параметров: решаемое дифференциальное уравнение Eq, начальные условия InCon, функция fnc (с указанием аргумента!), относительно которой решается уравнение, малый параметр epsilon и целое неотрицательное число п, определяющее порядок разложения по малому параметру при нахождении решения.

После каждого из перечисленных параметров стоит двойное двоеточие (::), после которого указывается тип параметра. Так, equation (этот тип имеют первый и второй параметры) определяет уравнение. Третий параметр имеет тип функция (function). Малый параметр имеет символьный тип (symbol) и, наконец, последний, пятый, параметр — целое неотрицательное число (тип nonnegint). Если в дальнейшем при вызове процедуры хотя бы один из ее параметров будет иметь иной тип, появится сообщение об ошибке. Если напрямую тип не указывать, будет предпринята попытка выполнить процедуру и, скорее всего, ошибка возникнет непосредственно в процессе выполнения — такие ошибки отслеживать сложнее.

Внимание!
Как сказано выше, второй параметр, определяющий начальные условия, имеет тип раввнство. Это справедливо только в том случае, если такое условие одно. Известно, что [для однозначного решения дифференциального уравнения число начальных условий должно равняться порядку этого уравнения. Таким образом, процедура может использоваться для решения задачи Коши только для уравнений первого порядка. Из сказанного, Однако, не следует, что с помощью создаваемой процедуры нельзя будет решать уравнения более высоких порядков. Просто в этих случаях такие решения придется искать в 1 общем виде, а присутствующие в решениях произвольные константы потом самостоятельно находить из начальных условий.

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

Первой определяется переменная х, которой в качестве значения присваивается операнд указанной аргументом функции fnc (op(fnc)), т.е. ее аргумент. Переменной у присваивается нулевой операнд функции fnc (команда op(0,fnc)); для функциональных зависимостей таким операндом является имя функции. Эти переменные, следовательно, определены так, чтобы тождественно выполнялось равенство y(x)=fnc. Кроме того, для удобства малый параметр запишем в переменную.

На следующем этапе в теле процедуры определяется другая процедура eps() — локальная, т.е. вне процедуры spDiffSol данная процедура будет недоступна! — при действии которой на аргумент t (т.е. процедурой определяется оператор) этот аргумент умножается на первый параметр процедуры г, возведенный в степень, определяемую вторым параметром процедуры i. После определения процедуры eps() формируется сумма F. Сумма эта операторная: слагаемыми суммы являются операторы. Каждый такой оператор есть композиция (оператор композиции) двух операторов: a[i] и eps(z,i), где i — индексная переменная, изменяющаяся в диапазоне от 0 до п. Действие оператора F на аргумент будет, согласно его описанию, следующим: на аргумент сначала действует оператор a[i], а результат этого действия умножается на z(в этом состоит действие оператора eps(z,i)). Общий результат формируется как сумма от действия каждого из слагаемых.

На заметку
В выражении для F использован оператор композиции §. Если G и f — функциональные операторы (т.е. G(x) и f (x) являются функциями), то оператор L, определенный как L:=G§f, действует следующим образом: L(x)=G(f (х)).

Описанные выше функциональные операторы a[i] будут формировать коэффициенты разложения искомой функции в ряд по малому параметру, и именно их нужно определить в процессе решения.

Далее вводится переменная eq, значение которой формируется следующим образом. Сначала в исходном решаемом уравнении Eq все слагаемые переносятся в левую часть (команда lhs(Eq)-rhs(Eq)), а в полученном таким образом выражении функция (fnc) заменяется с помощью процедуры subs() разложением в ряд. Разложение в ряд получается в результате действия оператора F на аргумент функции fnc (переменная х, а результат равен F(x)). После этого выражение eq раскладывается в ряд (процедура series()) по малому параметру в окрестности нуля (это нужно сделать на тот случай, если в выражении присутствуют отличные от степенных функции малого параметра). Остаток ряда должен иметь порядок п+1. Далее выполняется преобразование в полиномиальный вид (процедура convert ()). Из этого выражения будут сформированы уравнения для определения коэффициентов a[i].

Правая часть уравнения inCon переносится влево, и это выражение присваивается переменной S. После этого в выражении S оператор у() заменяется оператором F(). Для того чтобы оператор F возымел действие, вызывается процедура value () и результат ее выполнения присваивается переменной S. Таким образом, переменная S является тем выражением, которое определяет начальное условие и в котором действие оператора у() заменено действием операторного ряда F().

Внимание!
Удобство подобного подхода с использованием операторного ряда при работе с начальными условиями состоит в то, что нет необходимости отдельно выделять начальную точку для аргумента. Однако это палка о двух концах — процедура сможет работать только с такими начальными условиями, которые записаны относительно самой функции, но не ее производных. При работе с производными возникает проблема, связанная с попыткой действия оператора производной D на оператор eps ().

Выражение с начальными условиями далее раскладывается в ряд по малому параметру и трансформируется в полином.

После выполнения описанных выше преобразований с учетом начального условия, приступаем непосредственно к решению задачи. Для этого задействуем оператор цикла. В теле этого оператора первой командой формируется уравнение для определения очередного коэффициента a[i](x). С этой целью приравниваем к нулю полиномиальные коэффициенты разложения выражения eq по малому параметру. Полученные таким образом уравнения присваи-I ваем в качестве значения переменным, названия которых формируются путем | объединения базового имени eqn и индексной переменной i, принимающей значения в диапазоне от 0 до п. Точно так же поступаем и с начальными условиями, только в этом случае к нулю приравниваются коэффициенты в выражении S, а соответствующие уравнения записываются в переменные КО, IC1, .., 1Сп. Названия переменных создаются, как и ранее, посредством оператора конкатенации (| |). Наконец, уравнение (eqn| |i) и соответствующее ему граничное условие (1С 11 i) объединяются во множество, которое присваивается в качестве значения переменной Eq_and_ic| |i.

Далее следует процедура dsolve(), с помощью которой уравнение eqn| |i с граничным условием IC||i (они являются элементами множества Eq_and_IC11i) решается относительно функции a[i)(x). Поскольку результат решения представляется в виде равенства, после нахождения решения следует, согласно равенству-решению, определить функцию a[i](x) (если точнее, то оператор a[ij). Результат действия оператора a[i] на аргумент х определяется правой частью равенства, полученного в результате решения уравнения. Ссылка на этот результат выполнена посредством переменной среды %, а правая часть равенства, соответственно, возвращается командой rhs(%). Чтобы из этого выражения выделить оператор, используем процедуру unapply(). На этом команды оператора цикла заканчиваются.

По завершении описанного выше оператора цикла все уравнения решены, и в переменную Res записывается результат вычисления F(x).

Внимание!
Если не использовать процедуру eval(), а просто указать Res:=F(x), то результат выполнения процедуры будет представляться через коэффициенты a[i], без непосредственного вычисления результата их действия на аргумент х.

Последней командой процедуры является Res, в результате выполнения которой на экран выводится приближенное решение уравнения.

После того как определена процедура, логично проверить ее на предыдущей задаче., решение которой уже известно. Для "чистоты эксперимента" заново определим уравнение и начальные условия, заменив неизвестную функцию на f, аргумент — на t, а малый параметр — на mu (в области ввода — это mu, а в области вывода — ц). Ниже приведены соответствующие команды.

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

Это уравнение второго порядка, а малый параметр epsilon при квадратичном по отклонению слагаемом определяет степень ангармонизма (если параметр epsilon равен нулю, определяемые уравнением Eq2 колебания называются гармоническими).

Так как мы имеем уравнение второго порядка, для нахождения однозначного решения следует задать два начальных условия. Сначала запишем начальное условие для отклонения.

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

Задача 5.9

Найти решение в виде ряда для уравнения ху"(х)+у'(х)+ху(х) = чальными условиями у(0) = 1(0) = 0. с начальными условиями y(0)=1, y(0)=0.

Это уравнение имеет решение — функцию Бесселя j(x). Ниже попытаемся найти приближенное решение и сравнить его с точным. Как обычно, задаем уравнение.

Начальные условия для этого уравнения имеют следующий вид.

Если решать данную задачу с помощью процедуры dsolvef), получим следующий результат.

Чтобы получить приближенное решение в виде ряда, используем все ту же процедуру dsolve() практически с теми же параметрами; в конце добавлена опция series. Решение будет представлено рядом.

Количество слагаемых в ряде определяется значением переменной среды Order; значение этой переменной задает порядок остатка. Например, если нужно, чтобы последнее слагаемое было степени 9 по х, переменной Order присваиваем значение на единицу больше (т.е. 10).

В этом случае приближенное решение будет следующим.

То, что последнее слагаемое имеет степень 8, объясняется просто: в разложении присутствуют только четные степени х, т.е. коэффициенты при нечетных степенях — нули.

Полученное приближенное решение сравним с точным. Для этого преобразуем приближенное решение в полиномиальный вид.

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

Ниже показаны кривые для точного и приближенного решений.

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

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

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

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

Задача 5.10

Найти решение уравнения u(х,0)=0, <x<+,t=>0. (нижний индекс означает производную), удовлетворяющее начальным условиям u(х,0)=0, <x<+,t=>0.

Рассматриваемое уравнение является уравнением гиперболического типа. Подобные задачи возникают при описании процесса распространения волн, например, по бесконечной струне. Функция u(x,t) определяет отклонение точки струны с координатой х в момент времени t. Соответствующее волновое уравнение (а — параметр задачи) имеет следующий вид.

Для решения этого уравнения воспользуемся процедурой pdsolve() и поручим следующее.

В данном случае функции _F1 () и _F2 () являются произвольными дважды дифференцируемыми функциями. Таким образом, общее решение уравнения iBqn представляется в виде суперпозиции двух функций с соответствующими аргументами. Следовательно, чтобы полностью решить задачу, необходимо [определить вид этих функций. Функции определяются из начальных условий. Но прежде задаем и () как функцию двух параметров х и t.

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

На заметку
Производная по второму аргументу функции u(x,t) вычисляется посредством оператора эенцирования с указанием в квадратных скобках индекса переменной, по которой числяется производная: D(2](u).

Полученное таким образом дифференциальное уравнение решаем относительно функции _Fl(x).

Видим, что функции _F1 и _F2 с точностью до знака аргумента и конанты С1 совпадают. Константу можно положить равной нулю (несложно доказать, что общности метода это не ограничит), а функцию _F1 обозначим как F.

Тогда естественно определить функцию F2 следующим образом.

Следовательно, искать решение уравнения нужно в таком виде.

В последнем выражении присутствует уже только одна неизвестная функция F. При этом первое слагаемое F(at+x) описывает волну, распространяющуюся влево, а слагаемое F(-at+x) — волну, которая распространяется вправо. Непосредственно функцию F будем искать из оставшегося неиспользованным начального условия для значения функции u(x,t) в начальный момент времени.

Согласно полученному выражению для функции u(x,t), в начальный момент (при t=0) она равна следующему.

С другой стороны, это функция f (х), т.е. начальное отклонение струны. Для определенности возьмем функцию f (x) в таком виде.

Первоначальный профиль струны, таким образом, имеет форму симметричного треугольника.

На заметку
Функция Heaviside(x) равна 1 при х>0 и 0 — в противном случае.

Функция F тогда равна следующему.

Чтобы отобразить решение графически, присвоим значение символьному параметру а.

Значение функции u(x,t) будет таким.

Динамику системы будем отслеживать с помощью процедуры отображения анимации animate() из пакета plots. Параметры этой процедуры практически такие же, как и у процедуры plot(), но есть некоторые особенности. Второй параметр (в данном случае — это t) является не пространственной, а временной переменной. Другими словами, он определяет интервал, на протяжении которого будет отображаться анимация. Поэтому, чтобы задать диапазон отображения по вертикальной координате, используют опцию view и указывают нужный диапазон (здесь от 0 до 0.5). Назначение прочих опций такое же, как и в процедурах двухмерной графики.

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

Деформация бесконечной струны. Второй кадр

Деформация бесконечной струны. Третий кадр

Деформация бесконечной струны. Четвертый кадр

Деформация бесконечной струны. Пятый кадр

Чтобы "оживить" картинку, следует ее выделить (щелкнуть на ней кнопкой мыши) и выбрать на появляющейся в этом случае контекстной панели анимации нужные кнопки (их описание можно найти в главе 1). То же можно сделать с помощью опций раскрывающегося меню. К сожалению, в книге нет возможности демонстрировать анимацию. Поэтому выше приведены только некоторые кадры. Первым показан профиль струны в начальный момент времени. Второй и третий кадры дают представление о том, как реализуется суперпозиция распространяющихся в разные стороны волн. На четвертом кадре волны разошлись и дальше распространяются независимо друг от друга. Это особенно хорошо видно на пятом кадре.

Практически так же решается задача и для полубесконечной струны — условие задачи такое же, но только в этом случае 0 < х < +оо. Кроме того, следует задать значение функции u(0,t) на левой границе. Если левый конец струны закреплен, то это значит, что u(0,t)=0. Именно такую ситуацию И рассмотрим.

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

Начальное отклонение возьмем в виде выгнутой вверх параболы.

Функцию f(x) получаем из функции fl(x) путем нечетного продолжения Последней, т.е. f (x)=fl(x) при х>0 и f (x)=-fl(-x) — в противном случае.

Это решение отображаем в динамическом режиме.

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

Полубесконечная струна в антисимметричном отображении. Второй кадр

Полубесконечная струна в антисимметричном отображении. Четвертый кадр

Полубесконечная струна в антисимметричном отображении. Шестой кадр .

Полубесконечная струна в антисимметричном отображении. Восьмой кадр

Выше на первом кадре показано начальное распределение отклонения J(x,0). Полезной является только правая (относительно вертикальной координатной оси) часть графика; левая часть вспомогательная. На втором кадре можно видеть причудливую фигуру, образующуюся при расхождении волн. Симметричная относительно начала координат структура является мнимой; это условная волна, при наличии которой колебания бесконечной струны осуществляются так же, как и полубесконечной.

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

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

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

Колебания полубесконечной струны. Первый (начальный) кадр

Колебания полубесконечной струны. Второй кадр

Колебания полубесконечной струны. Третий кадр

Колебания полубесконечной струны. Четвертый кадр

Колебания полубесконечной струны. Пятый кадр

Колебания полубесконечной струны. Восьмой кадр

Колебания полубесконечной струны. Десятый кадр

Колебания полубесконечной струны. Двенадцатый кадр

Колебания полубесконечной струны. Четырнадцатый кадр

Колебания полубесконечной струны. Шестнадцатый кадр

Начальное состояние такое же, но весь процесс выглядит теперь несколь-иначе. Кадры с первого по пятый позволяют в деталях проследить дина-1ку расхождения двух волн.

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

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

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

Задача 5.11

Решить задачу о колебаниях конечной струны

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

Определим уравнение.

Решение будем искать методом разделения переменных. Этот метод подразумевает, что поиск решения осуществляется в виде произведения нескольких (в данном случае двух — по количеству переменных) функций, каждая из которых зависит от одного аргумента. В этом случае при вызове процедуры pdsolve() следует использовать опцию HINT, указав ее значение равным X(x)*T(t); именно в таком виде будем искать функцию u(x,t)=X(x)*T(t) (X(x) и T(t) являются неизвестными функциями одной переменной). Таким образом, имеем следующее.

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

Задаем уравнение для функции Х(х), заменив в нем для удобства переменную среды

На заметку
Стоит обратить внимание на способ, которым с помощью оператора формирования последовательности ($) задана вторая производная. Здесь использована та особенность, что результатом выполнения команды '$' (х,2) является последовательность х,х.

Решаем это уравнение относительно Х(х), приняв во внимание одно из начальных условий, а именно: поскольку u(0,t)=X(0)T(t)=0, то Х(0)=0. Получаем следующее.

Параметр Я. должен быть таким, чтобы выполнялось и условие Х(1)=0. Но прежде чем решать соответствующее уравнение (относительно X), присвоим переменной среды _EnvAllSolutions, отвечающей за поиск всех решений уравнения, значение true.

В этом выражении переменная среды _Z1 "нумерует" собственные числа.

Внимание!
Зыше переменная среды _Z1 в области вывода содержит знак тильды. Это значит, что переменная может принимать далеко не любое значение; на нее наложены ограничения Что это за ограничения, можно узнать, воспользовавшись процедурой about (). В данном [случае в результате выполнения команды about(_Zl) появится сообщение Originally renamed _Z1": is assumed to be: integer, что в переводе значит следующее I исходном варианте _Z1, переименовано в _2Г: предполагается, что имеет целочисленный тип (integer)".

Поэтому можем задать зависимость, определяющую собственные числа краевой задачи (именно так называются найденные выше значения для .lambda).

После этого определим собственные функции — такие функции, которые соответствуют собственным числам задачи.

Решение искалось такое, чтобы оно автоматически удовлетворяло одному из начальных условий (равенство нулю производной в начальный момент). Поскольку возможные значения для lambda определены выше (зависимость nu(n)), заменяем параметр lambda на nu(n) и определяем базовые функции, соответствующие собственным числам.

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

Отсюда, в частности, следует, что искомые коэффициенты совпадают с коэффициентами разложения функции f (х) в ряд Фурье по синусам на интервале от 0 до 1. Ниже эти коэффициенты будут вычислены, однако прежде следует указать, что индекс п коэффициента А является целым положительным числом.

Внимание!
При определении коэффициентов разложения используются индексы. Поскольку все присвоения символьные, массив в этом случае не задается, равно как и функция от ин-зкса. Другими словами, вызовом А[ 1 ] выражение для первого коэффициента ряда почить не удастся. Чтобы это было возможно, следует определять коэффициенты как функции индекса.

Выражение для функции u(x,t) будет иметь следующий вид.

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

Кроме того, следует учесть, что ряд для функции u(x,t) бесконечный, поэтому его следует ограничить — нужно оставить конечное число слагаемых. Соответствующее выражение определим следующим образом (коэффициенты записаны в явном виде, N — число слагаемых в ряде).

Теперь с помощью процедуры animate () воспроизводим процесс колебаний струны.

Колебания струны конечной длины. Первый (начальный) кадр

Колебания струны конечной длины. Четвертый кадр

Колебания струны конечной длины. Седьмой кадр

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

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

Заключительные замечания

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

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

Контрольные вопросы

1. Какая процедура используется при решении обыкновенных дифференциальных уравнений? a)diff() 6) solve(] в) dsolve(); г) int(); 2. Какие из приведенных ниже команд корректны? Каков результат их выполнения? а) dsolve(D(y)(t)=0,y(t)); б) solve(D(y)(t)=0,t); в) dsolve(y(t)=5,y(t)); г) solve(y(t)=5,y(t)); 3. Действие процедуры isolate() состоит в том, что а) указанное первым параметром выражение заменяется вторым параметром-выражением; б) указанное вторым параметром выражение заменяется первым параметром-выражением; в) указанное первым параметром выражение решается относительно второго параметра-выражения; г) правая и левая части уравнения меняются местами. 4. Переменная S описана как S:=x*sin(x). Каков результат выполнения команды F:=unapply(S,x)? а) F:=S; б) F:=x->x*sin{x); в)S:=x; г) команда некорректна 5. Если имеет место Fl:=x->xA2, F2:=x->sin(x) и F:=F2§F1, что получится в результате выполнения команды F(x)? а) xA2*sin(x); б) sin(х)А2; в)sin(xA2); г) команда некорректна 6. Какая процедура используется при решении дифференциальных уравнений в частных производных? а) dsolve (); в) pdsolve (); б) solve (); г) pdesolve (); 7. Какая из перечисленных ниже процедур используется для отображения двухмерной анимации? a) plot (); 6)plot3d(); в) animate (); г) animate3d()