Дифференциальные уравнения

Решение обыкновенных дифференциальных уравнений Приближенные методы решения дифференциальных уравнений Уравнения в частных производных Заключительные замечания Контрольные вопросы

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

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()

На заметку

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

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

На заметку

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

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

На заметку

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

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

На заметку

В выражении для 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().

На заметку

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

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

На заметку

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

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

На заметку

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

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

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

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

присваивается значение, на самом деле

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

Как сказано выше, второй параметр,

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

с использованием операторного ряда при

Удобство подобного подхода с использованием операторного ряда при работе с начальными условиями состоит в то, что нет необходимости отдельно выделять начальную точку для аргумента. Однако это палка о двух концах — процедура сможет работать только с такими начальными условиями, которые записаны относительно самой функции, но не ее производных. При работе с производными возникает проблема, связанная с попыткой действия оператора производной 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, а в области вывода — ц). Ниже приведены соответствующие команды.

в области вывода содержит знак

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

При определении коэффициентов разложения используются

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

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

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

Решение задачи

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

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

и начальных условий. Уравнение решается

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

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

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

и не набирать его при

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

к которой определяется данная функциональная

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

в таком случае дается следующей

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

Полученное выражение можно упростить, причем

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

к которому стремится скорость при

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

в пределе бесконечно больших времен

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

Данное уравнение можно решить

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

с математической точки зрения, Произвольными

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

в начальный момент отклонение маятника

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

и _С2 присвоить соответствующие значения

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

в выражении экспонент смущать не

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

ю равной 1, переходим фактически

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

если точнее, то поверхность) зависимости

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

с увеличением силы трения амплитуда

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

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

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

в виде интеграла, который вычислить

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

В результате получено два решения,

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

Однако не всегда удается получить

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

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

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

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

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

Решение задачи

Преобразуем это уравнение, перенеся все

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

Теперь можно собрать слагаемые при

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

выбираются так, чтобы коэффициенты при

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

к нулю коэффициент при первой

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

прежде чем решать уравнение Eq_l,

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

в какой они вводились. Причина

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

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

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

В этом месте функция

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

в приведенном выше вызове изменить

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

В данном случае, перед присваиванием,

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

Последовательность действий такая

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

в самом начале задачи разложению

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

Задача несколько усложняется, если малый

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

в этом случае определяются через

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

Решать уравнение будем таким образом,

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

Формировать ряд по малому параметру

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

Слагаемые добавлять будем

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

В рамках этого цикла на

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

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

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

Преобразуем это уравнение, перенеся все

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

Далее левую часть уравнения Eq

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

к нулю, получаем последовательность уравнений

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

й степени epsilon определяется процедурой

Коэффициент при 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). Если в дальнейшем при вызове процедуры хотя бы один из ее параметров будет иметь иной тип, появится сообщение об ошибке. Если напрямую тип не указывать, будет предпринята попытка выполнить процедуру и, скорее всего, ошибка возникнет непосредственно в процессе выполнения — такие ошибки отслеживать сложнее.

Первой определяется переменная х, которой

Первой определяется переменная х, которой в качестве значения присваивается операнд указанной аргументом функции 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)). Общий результат формируется как сумма от действия каждого из слагаемых.

с помощью созданной выше процедуры

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

а малый параметр epsilon при

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

Еще один распространенный метод поиска

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

Начальные условия для этого уравнения

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

Если решать данную задачу

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

в виде ряда, используем все

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

в ряде определяется значением переменной

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

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

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

что последнее слагаемое имеет степень

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

Результат преобразования возьмем за основу,

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

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

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

в окрестности точки 0, где

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

Для решения этого уравнения воспользуемся

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

В данном случае функции _F1

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

Далее воспользуемся тем, что производная

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

Решение задачи

с точностью до знака аргумента

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

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

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

искать решение уравнения нужно

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

В последнем выражении присутствует уже

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

С другой стороны, это функция

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

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

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

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

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

Значение функции

Значение функции 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. После этого начинается процесс отражения от границы. На десятом кадре видно, как деформируется набегающая волна.
Двенадцатый кадр позволяет убедиться, что деформация струны у границы практически отсутствует. На четырнадцатом кадре видно отраженную волну, которая уходит от стенки. На шестнадцатом кадре представлена полностью сформированная отраженная волна. После этого обе волны движутся в одном "направлении (с одинаковой скоростью).
Рассмотренный выше подход называется методом распространяющихся волн или методом Даламбера. Данный метод приемлем, в основном, при ешении задач для бесконечных и полубесконечных областей. Ситуация несколько усложняется, если рассматривать струну конечной длины. Об этом следующая задача.

Решение будем

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

в результате выполнения команды выражении

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

Решение задачи

Параметр Я. должен быть таким,

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

В этом выражении переменная среды

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

После этого определим собственные функции

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

Решение искалось

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

находятся из начального условия. Согласно

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

в частности, следует, что искомые

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

Решение задачи

Проанализируем полученное решение, отобразив его

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

Кроме того, следует учесть, что

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

с помощью процедуры animate

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

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

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

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

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

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

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

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

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

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

Составить дифференциальное уравнение, описывающее падение

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

и решить дифференциальное уравнение для

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

В данном случае укажем уравнение,

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

Найти общее решение уравнения

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

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

Найти решение задачи Коши
Задаем непосредственно уравнение.

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

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

в виде многочлена второго порядка

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

в виде многочлена второго порядка

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

в виде ряда для уравнения

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

Рассматриваемое уравнение является уравнением гиперболического

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

Уравнение не изменилось, однако здесь

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

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

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

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

---

---