Рассмотрим двумерную задачу о движении пушечного ядра в гравитационном поле с учетом сопротивления воздуха. Уравнение Ньютона запишем в виде ДУ первого порядка для координат-импульса (уравнение Гамильтона), где m – масса ядра, g – ускорение свободного падения, b – коэффициент сопротивления воздуха, пропорциональный скорости ядра.

Сопротивление воздуха действует только в направлении x, а гравитация – только в направлении y.

• restart;
• Momx:=m*diff(vx(t),t)=-b*vx(t);

• Momy:=m*diff(vy(t),t)=-m*g-b*vy(t);

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

• IC:=vx(0)=v0*cos(theta),vy(0)=v0*sin(theta);

• sol:=dsolve({Momx,Momy,IC},{vx(t),vy(t)});

• assign(sol);

Можно проверить, что найденные решения на самом деле удовлетворяют дифференциальным уравнениям:

• simplify(Momx);

• evalb(%);

• simplify(Momy);

• evalb(%);

Разложим решение в ряд Тейлора вблизи t = 0. Это покажет, как скорость начинает отличаться от ее начального значения:

• taylor(vx(t),t=0,3);

• taylor(vy(t),t=0,3);

Изменение зависит линейно от времени, причем в нем есть вклад от гравитации. В части, описывающей сопротивление, есть зависимость от компонент начальной скорости. Квадратичные члены – сложнее.

Интересно изучить это решение. Наиболее общая проблема – зависимость траектории от наклона ядра при фиксированном сопротивлении b. В реальности нужно еще учитывать скорость и направление ветра.

Размерности параметров – в СИ (или в СГС). Имеем:

• g:=9.81;

Масса = 1 кг:

• m:=1;

Коэффициент сопротивления (в кг/с):

• b:=0.5;

• vx(t);

• simplify(%);

• simplify(vy(t));

Пусть угол наклона можно изменять, а начальная скорость фиксирована. Ее значения будем измерять в м/с (система СИ):

• v0:=50;

Выберем параметр наклона:

• simplify(vy(t));

Чтобы получить функцию, а не выражение, применим команду unapply:

• Vy:=unapply(simplify(vy(t)),theta,t);

Теперь можно проинтегрировать уравнение движения для компонент вектора. Переменную для времени обозначим s вместо t. Применим команду unapply, чтобы отобразить результат:

• X:=unapply(int(Vx(theta,s),s=0..t),theta,t);

• X(0.5,1);

• Y:=unapply(int(Vy(theta,s),s=0..t),theta,t);

(Результат вычислений посмотрите сами.)

• Y(0.5,1);

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

• plot([X(0.5,t),Y(0.5,t),t=0..10],x=0..X(0.5,10),y=0..20);

Рис. 4.5.1

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

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

Пересечение с поверхностью получается из условия

• t0:=solve(Y(0.5,t)=0,t);

Конечно, в начале движения есть пересечение с y = 0. Порядок, в котором ищутся два решения, можно изменять. Для параметрического графика неважно, в каком порядке отслеживается траектория. Не надо указывать пределы осей, которые автоматически подгоняются под решение. Пример показан на рис. 4.5.2.

• plot([X(0.5,t),Y(0.5,t),t=t0[1]..t0[2]]);

Рис. 4.5.2

Нарисуем семейство решений для разного наклона ядра. Цель – определить диапазон и максимально высоту, достигаемую осколками.

• with(plots):
• imax:=5;

Чтобы различать решения, задаем цветовую шкалу (рис. 4.5.3):

• coltab:=[maroon,red,magenta,yellow,green,blue,violet,
  brown,gray,black];
• for i from 1 to imax do:
• theta:=evalf((i-0.5)/imax*Pi/2);
• t0:=solve(Y(theta,t),t);
• P[i]:=plot([X(theta,t),Y(theta,t),t=t0[1]..t0[2]],
  color=coltab[i]): od:
• display(seq(P[i],i=1..imax));

Рис. 4.5.3

Можно построить последовательность графиков как анимацию.

• display(seq(P[i],i=1..imax),insequence=true);

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