ДУ второго порядка часто встречаются в физике: закон Ньютона, уравнение Шредингера, диффузия...

• restart:
• E2:=diff(y(t),t$2)=-g;
• S2:=dsolve({E2,y(0)=y[0],D(y)(0)=v[0]},y(t));

Это очевидный результат.

Заметьте несколько деталей.

Для ввода начальных условий для производной от y применена D(y), что означает dy/dt (D – это общий символ производной для Maple), затем добавили (0), чтобы сказать, что производная вычислена при t=0. Это означает, что diff и D – это производные, но между ними есть тонкое различие, которое надо понимать. Оно обсуждено в разделе о производных (см. D и diff). Поэтому или запоминайте синтаксис команды dsolve (выше), или запомните, где его легко найти. Кроме того, учтите, что применяемая в пакете Physics команда diff может давать другой результат – см. описание команд пакета Physics.

Это второе знаменитое ДУ второго порядка.

• E3:=diff(y(t),t$2)=-omega^2*y(t);

Maple решает его в общем виде с произвольными константами:

• dsolve(E3,y(t));

Предупреждение: эта форма решения опасна, поскольку нет гарантии, что при каждом запуске рабочего листа с командой dsolve в таком виде получим на выходе _C1 с синусом и _C2 с косинусом. В Maple есть собственная логика такого выбора, которая даже для более простых ДУ кажется случайной. Поэтому не давайте Maple выбирать, делайте это сами. Для этого скопируйте решение мышкой в новую строку и выберите неизвестные коэффициенты, например избавьтесь от неудобного y(t)=, заменив его естественным присвоением:

• ysol:= A*sin(omega*t)+B*cos(omega*t);

Теперь каждый раз расчет величин A и B будет срабатывать одинаково, потому что «насильственно» установлено, что A – с синусом, а B – с косинусом.

Если константы определяются начальными условиями, есть форма команды dsolve, которая будет определять их автоматически. Например, нужно решить ДУ E3 с начальными условиями y(0) = 1 и dy/dt(0) = 2. Тогда пишем:

• s:=dsolve({E3,y(0)=1,D(y)(0)=2},y(t));

Применим assign, чтобы Maple дал выражение для y(t):

• assign(s);

Теперь построим график решения:

• plot(y(t),t=0..20);

Не работает. Рамка есть, но нет функции. И нет догадок, что не так. Для отладки надо бы подставить несколько значений t в y(t). Но вспомним, что y(t) выглядит как выражение, но не совсем им является, и что подстановка значений для t запустит его. Рассмотрим y(t) с этой точки зрения и прикинем, можно ли увидеть, что не так:

• y(t);

Команда plot даст численные значения t, но ведь просто включено в выражение и не имеет численного значения. Это и есть ошибка: мы не задали численные значения !

OК, теперь сделаем так:

• plot(subs(omega=1,y(t)),t=0..20);

Решение этой задачи было достаточно легким. Но так же легко Maple решает и более сложные задачи.

Усложним условия и рассмотрим

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

• restart:with(DEtools):
• Eq:=diff(y(t),t$2) + diff(y(t),t)/tau + omega^2*y(t)=0;

и мы получаем общее решение:

• dsolve(Eq,y(t));

Выглядит впечатляюще.

Где же тогда появятся синусы и косинусы? Проблема в том, что Maple не знает, насколько велики и τ. Подумаем с точки зрения физического смысла. Пусть естественная частота порядка 2π, тогда период N = 2π/ = 1 с. Теперь пусть осциллятор (подразумеваем – маятник) находится в моторном масле при 50 градусах ниже нуля. Тогда характеристическое время затухания 0.05 с (трение столь велико, что движение остановится за 0.05 с). Интуиция должна подсказать, что если толкнуть маятник назад и освободить его, то ожидается некоторое покачивание, т. е. маятник медленно проплывет к своему вертикальному положению и останется в нем. Посмотрим, что скажет Maple:

• omega:=2*Pi;tau:=.05;
• Y:=dsolve({Eq,y(0)=1,D(y)(0)=0},y(t));
• plot(rhs(Y),t=0..20);

OK, Maple согласен с интуицией.

Теперь изменим затухание и представим, что произойдет, если нагреть масло или использовать вместо него WD-40 (жидкая смазка) или, может быть, воздух. Маятник начнет колебаться с малым затуханием. Например, если τ = 20, то Maple дает:

• omega:=2*Pi;tau:=20;
• Y:=dsolve({Eq,y(0)=1,D(y)(0)=0},y(t));
• plot(rhs(Y),t=0..20);

Теперь он колеблется дольше и медленно затухает. Но какое чудо в решении превратило экспоненту в синусы и косинусы? Вспомним формулу Эйлера:
т. е. в экспоненте константы стали комплексными, поэтому решение перешло от распада к колебаниям.

Это поднимает вопрос о величине τ, при которой происходит переход от чистого затухания к затухающим колебаниям. Посмотрите на общий вид решения и увидите, что в экспонентах есть квадратный корень Если аргумент под корнем положителен, то оба фундаментальных решения в сумме есть затухающие экспоненты, и получаем затухание. Если аргумент под корнем отрицателен, то результат комплексный и экспонента комплексная, т. е. это синусы и косинусы, тогда получаем колебания. Переход между двумя режимами происходит, если квадратный корень = 0, т. е. когда τ = 1/(2). Этот особый случай называется критическим затуханием.

Пусть Maple решит задачу о гармоническом осцилляторе с затуханием для случая критического затухания при = 1 и нарисует график функции y(t) для t = 0..20 при начальных условиях y(0) = 1 и dy/dt(0) = 0. График будет выглядеть в точности как для осциллятора с затуханием при полном отсутствии намеков на то, что это точно граница между затуханием и колебаниями. Для проверки этого увеличьте τ на 2 % и снова постройте график. Растяните окно вниз, чтобы можно было видеть, что y(t) становится немного отрицательным, показывая начало колебаний (нужны очень серьезные размеры окна).

Примените Maple, чтобы получить общие решения ДУ второго порядка:

(a) ,

(b) ,

(c) .