Компьютерные лабораторные занятия по теоретической механике. Ч. 2 (110

Ответ:

Задача 9.24

Составить уравнение Гамильтона-Якоби для частицы, движущейся в однородном поле тяжести и найти полный интеграл этого уравнения (в декартовой системе координат).

Ответ:

Задача 9.25

Найти полный интеграл уравнения Гамильтона-Якоби и закон движения (в квадратуре) математического маятника массы m длины l.

Решение.

Из функции Гамильтона составим уравнение Гамильтона-Якоби:

> H:p^2/(2*m*l^2)-m*g*l*cos(phi);

Так как уравнение не содержит переменную t в явном виде, то решение ищем в виде

S('; t) = Et + S0('):

> eqS:S=-E*t+S0(phi);

> subst(eqS,eqHJ)$ eq:ev(%,diff);

Решим дифференциальное уравнение:

> eqv:solve(eq,diff(S0(phi),phi));

> ode2(eqv[2],S0(phi),phi);

> subst(%,eqS);

Здесь постоянная E появилась при решении дифференциального уравнения. Согласно методу Гамильтона-Якоби, производная @S=@E òî-

же равна постояной, которую удобно обозначить t0:

> diff(part(%,2),E)=-t0;

31

Ответ:

Задача 9.26

Найти полный интеграл уравнения Гамильтона-Якоби и закон движения для линейного гармонического осциллятора массы m с частотой

!.

Ответ:

r

2E

x(t) = m!2 sin !(t t0):

10. Гидродинамика

10.1. Идеальная жидкость

Задача 10.1

Цилиндрический сосуд вместе с находящейся в нем несжимаемой жидкостью вращается в однородном поле тяжести ( g) вокруг своей оси

с постоянной угловой скоростью !. Плотность жидкости , внешнее дав-

ление на жидкость равно нулю. Проверить выполнение уравнения непрерывности.

Решение.

Считаем, что вектор угловой скорости ! направлен вдоль оси Oz.

> load(vect);

> r:[x,y,z]; Omega:[0,0,omega];

Поле скоростей внутри сосуда найдем по формуле: v = [! r] :

> v:express(Omega~r);

Итак, вектор скорости зависит от координат следующим образом:

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

32

которой подчиняется движение жидкости в данной задаче, применим, например, при = 1:

> plotdf([v[1],v[2]],[x,y],[parameters,"omega=1"]);

Искомое поле скоростей отображено на рисунке. Уравнение непрерывности

при постоянной плотности жидкости соответствует div j = 0, что проверяется введением следующих строк:

> div(v);

> ev(express(%),diff);

Задача 10.2

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

Решение.

Векторная конструкция в левой части уравнения Эйлера

@@tv + (vr) v = 1 grad p + g

33

запрограммируем ее отдельно (в сумме i = 1; 2; 3):

> vnabla(t):=sum(v[i]*diff(t,r[i]),i,1,3);

Далее предполагается, что в переменную v уже введен вектор скорости (10.1). Подставляем уравнение Эйлера в переменную eq, считая, что

ускорение свободного падения направлено против оси Oz:

> G:[0,0,-g];

> eq:diff(v,t)+vnabla(v)=-express(grad(p))/rho+G;

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

(обозначив dp1, dp2, dp3), это удобно сделать в цикле:

> for n from 1 thru 3 do

dp[n]:solve(part(eq,1)[n]=part(eq,2)[n],'diff(p,r[n]))[1];

Итоговое решение собирается из интегралов и общей постоянной P0:

> P:integrate(part(dp[1],2),x)

+integrate(part(dp[2],2),y)+integrate(part(dp[3],2),z)+p0;

Итак, поле давления равно

Уравнение свободной поверхности жидкости найдем из условия p = 0, например, при

! = 1; p0 = 10; = 1; g = 10 :

> Z:subst([omega=1,P0=10,rho=1,g=10],solve(P=0,z)[1]);

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

команд plot3d èëè wxplot3d :

> plot3d(part(Z,2), [x,-5,5], [y,-5,5])$

Для наглядного отображения границ цилиндра удобно выполнить построение в цилиндрической системе координат:

> load(draw);

> draw3d(cylindrical(sqrt(20*z-20),z,0,2.5,az,0,2*%pi));

Задача 10.3

Сосуд заполнен водой до высоты H. В дне сосуда открыли отверстие малой площади . Найти, через какое время вся вода вытечет из сосуда, если известна его форма:

34

а) цилиндр радиуса R;

б) конус с углом 2 при вершине; в) перевернутый конус с углом 2 при вершине;

г) сфера радиуса R, причем H = 2R.

Решение для пункта б).

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

равна 2gh:

> v:sqrt(2*g*h);

Изменение объема воды в сосуде dV за бесконечно-малый промежуток времени dt можно записать двумя способами:

dV = S dh; dV = v dt;

ãäå S площадь свободной поверхности воды. Выразим производную времени по изменяющейся высоте:

dhdt = vS ;

здесь знак минус означает уменьшение высоты воды в сосуде со временем. Площадь поверхности воды в цилиндре равна

S = r2; r = (H h) tg :

> r:(H-h)*tan(alpha);

> S:%pi*r^2;

Подставим площадь в дифференциальное уравнение для t:

> eq:'diff(t,h)=-S/sigma/v;

и решим его с учетом начальных условий h = H ïðè t = 0:

> ic1(ode2(eq, t, h),h=H,t=0);

Упростив радикалы, найдем, при каком времени t высота воды станет равно нулю:

> radcan(%);

> subst(h=0,%);

10.2. Вязкая жидкость

Задача 10.4

Стационарный поток несжимаемой вязкой жидкости движется между неподвижными параллельными плоскостями, находящимися на рас-

стоянии l друг от друга. Вдоль направления течения поддерживается

35

постоянный перепад давления, равный p= L ( L длина отрезка, на котором давление изменяется на p). Пренебрегая силой тяжести, найти

поле скоростей. Является ли такое течение потенциальным? Решение.

Обозначив p0 давление на плоскости x = 0, запишем зависимость давления от координаты x:

p

p = p0 L x:

> p:p0-dp/dL*x;

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

только от координаты y:

> load(vect);

> r:[x,y,z]; v:[vx,0,0];

> depends(vx,y);

В уравнении Навье-Стокса для несжимаемой жидкости

@@tv + (vr) v = grad p + g + v

конструкцию с оператором grad снова вычисляем по формуле (10.2). В

отсутствие поля тяготения имеем

> rho*(diff(v,t)+sum(v[i]*diff(v,r[i]),i,1,3))=

-grad(p)+eta*laplacian(v);

В последнем равенстве выполним действие дифференциальных операторов и приравняем между собой x-компоненты левой и правой частей:

> ev(express(%),diff);

> eq:part(%,1)[1]=part(%,2)[1];

Решим полученное уравнение (eq) относительно vx(y):

> sol:ode2(eq,vx,y);

Из условия неподвижности стенок подставляем граничные условия:

> bc2(sol,y=0,vx=0,y=l,vx=0);

factor(%);

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

> v[1]:part(%,2);

и возьмем rot v:

> ev(express(curl(v)),diff);

36

factor(%);

Для наглядного отображения поля скоростей, например, при = 0:1

èl = 10, снова применим пакет plotdf:

> load(plotdf);

plotdf([v[1],v[2]],[x,y],

[parameters,dp=1,eta=0.1,l=10,dL=1"],[y,0,10]);"

1 p

Ответ: v = 2 Ly(l y) ex.

Поле скоростей непотенциально, так как rot v = (2y l) p ez 6= 0.

2 L

Задача 10.5

Решить предыдущую задачу при условии, что p= L = 0, а нижняя плоскость движется со скоростью u.

11.Функции стандартных библиотек Maxima

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

> example(integrate);

выдаст несколько типовых примеров применения команды integrate, а после ввода строки

37

> describe(integrate);

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

11.1.Аналитические преобразования алгебраических выражений

Наиболее часто встречающимися на практике аналитическими преобразованиями являются разложения на множители и раскрытие скобок. Эти преобразования осуществляются соответственно функциями factor

èexpand

> expand((y+4)^5);

y5 + 20 y4 + 160 y3 + 640 y2 + 1280 y + 1024

> factor(x^4-256);

(x 4) (x + 4) x2 + 16

> factor(2016);

25 32 7

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

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

> z1:1/(x+1)^2-2/(x+1)+2/(x+2);

> z2:ratsimp(z1);

> z3:partfrac(z2,x);

38

Для упрощения тригонометрических выражений используются следующие функции: trigsimp тригонометрическое упрощение, trigrat (синтаксис вызова trigrat(expr)) приводит заданное тригонометрическое выражение expr к канонической упрощ¼нной квазилинейной форме. Это выражение рассматривается как рациональное, содержащее sin, cos, tan, аргументы которых есть линейные формы некоторых переменных. Всегда, когда возможно, заданное выражение линеаризуется.

Для раскрытия скобок в тригонометрическом выражении используется функция trigexpand, которая использует формулы преобразования сумм двух углов для представления введенного выражения в как можно более простом виде, где в качестве аргумента только одна переменная. Также отметим функцию тригонометрического приведения trigreduce. Она преобразует тригонометрическое выражение к сумме элементов, каждый из которых содержит только sin или cos.

> ev(sin(3*x)/sin(x)+x,trigexpand,expand);

> trigreduce(%);

> trigsimp(%);

> trigexpand(sin(y+10*x));

2 cos (2 x) + x + 1

cos (10 x) sin (y) + sin (10 x) cos (y)

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

> k1:(log(x+x^2)-log(x))^a/log(1+x)^(a/2);

> radcan(k1);

log x2 + x log (x) a

a

log (x + 1)2

a

log (x + 1)2

> k2:log(1+2*a^x+a^(2*x))/log(1+a^x);

39

> radcan(k2);

log a2 x + 2 ax + 1

log (ax + 1)

2

> k3:(%e^x-1)/(%e^(x/2)+1);

> radcan(k3);

ex 1

x

e2 + 1

x

e2 1

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

> fif:(sin(x)+x**2)/(x+5)/(x-5);

x2 + sin (x) (x 5) (x + 5)

> num(fif);

sin (x) + x2

> denom(fif);

(x 5) (x + 5)

Функция combine приводит сумму дробей к общему знаменателю

> z:1*f/2*b+2*c/3*a+3*f/4*b+c/5*b*a;

> combine(z);

4 (3 a b c + 10 a c) + 75 b f

60

Для преобразования выражений, содержащих вложенные квадратные (или более высоких степеней) корни используется функция rootscontract

> rootsconmode:true;

> rootscontract(x^(1/2)*y^(1/4));

> rootscontract(x^(1/2)*y^(1/3));

true

40