Шаг 1: Выведите уравнение движения

Маятник является простой механической системой, которая следует дифференциальному уравнению. Маятник первоначально покоится в вертикальном положении. Когда маятник перемещается на угол $$\theta$$ и освобождается, сила тяжести тянет его назад к своему положению покоя. Его импульс заставляет его перерегулировать и прийти к углу $$-\theta$$ (если нет сил трения) и так далее. Восстанавливающая сила вдоль движения маятника от силы тяжести $$-mg\sin \theta$$. Таким образом, согласно второму закону Ньютона, масса, умноженная на ускорение, должна равняться $$-mg\sin \theta$$.

syms m a g theta(t)
eqn = m*a == -m*g*sin(theta)

eqn(t) = $$a\hspace{0.17em}m=-g\hspace{0.17em}m\hspace{0.17em}\sin \left(\theta \left(t\right)\right)$$

Для маятника с длиной $$r$$, ускорение маятника боба равно угловому ускорению $$r$$.

$$a=r\frac{d^2\theta }{d{t}^2}$$.

Замена на $$a$$ при помощи subs.

syms r
eqn = subs(eqn,a,r*diff(theta,2))

eqn(t) = 
$$m\hspace{0.17em}r\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }\theta \left(t\right)=-g\hspace{0.17em}m\hspace{0.17em}\sin \left(\theta \left(t\right)\right)$$

Изолируйте угловое ускорение в eqn при помощи isolate.

eqn = isolate(eqn,diff(theta,2))

eqn = 
$$\frac{{\partial }^2}{\partial t^2}\mathrm{ }\theta \left(t\right)=-\frac{g\hspace{0.17em}\sin \left(\theta \left(t\right)\right)}{r}$$

Соберите константы $$g$$ и $$r$$ в один параметр, который также известен как естественная частота.

${\omega }_0=\sqrt{\frac{\mathit{g}}{\mathit{r}}}$.

syms omega_0
eqn = subs(eqn,g/r,omega_0^2)

eqn = 
$$\frac{{\partial }^2}{\partial t^2}\mathrm{ }\theta \left(t\right)=-{{\omega }_0}^2\hspace{0.17em}\sin \left(\theta \left(t\right)\right)$$

Шаг 2: Линеаризируйте уравнение движения

Уравнение движения нелинейно, поэтому трудно решить аналитически. Предположим, что углы малы, и линеаризируйте уравнение с помощью разложения Тейлора $$\sin \theta$$.

syms x
approx = taylor(sin(x),x,'Order',2);
approx = subs(approx,x,theta(t))

approx = $$\theta \left(t\right)$$

Уравнение движения становится линейным уравнением.

eqnLinear = subs(eqn,sin(theta(t)),approx)

eqnLinear = 
$$\frac{{\partial }^2}{\partial t^2}\mathrm{ }\theta \left(t\right)=-{{\omega }_0}^2\hspace{0.17em}\theta \left(t\right)$$

Шаг 3: Решение уравнения движения аналитически

Решить уравнение eqnLinear при помощи dsolve. Задайте начальные условия как второй аргумент. Упростите решение, приняв $${\omega }_0$$ реально использовать assume.

syms theta_0 theta_t0
theta_t = diff(theta);
cond = [theta(0) == theta_0, theta_t(0) == theta_t0];
assume(omega_0,'real')
thetaSol(t) = dsolve(eqnLinear,cond)

thetaSol(t) = 
$${\theta }_0\hspace{0.17em}\cos \left({\omega }_0\hspace{0.17em}t\right)+\frac{{\theta }_{\mathrm{t0}}\hspace{0.17em}\sin \left({\omega }_0\hspace{0.17em}t\right)}{{\omega }_0}$$

Шаг 4: Физическая значимость $${\omega }_0$$

Термин $${\omega }_{0}t$$ называется фазой. Функции косинуса и синуса повторяются каждый $$2\pi$$. Время, необходимое для изменения $${\omega }_{0}t$$ около $$2\pi$$ называется временным периодом.

$$T=\frac{2\pi }{{\omega }_0}=2\pi \sqrt{\frac{r}{g}}$$.

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

Шаг 5: Постройте движение маятника

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

Определите физические параметры:

gValue = 9.81;
rValue = 1;
omega_0Value = sqrt(gValue/rValue);
T = 2*pi/omega_0Value;

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

theta_0Value  = 0.1*pi; % Solution only valid for small angles.
theta_t0Value = 0;      % Initially at rest.

Замените физические параметры и начальные условия в общее решение.

vars   = [omega_0      theta_0      theta_t0];
values = [omega_0Value theta_0Value theta_t0Value];
thetaSolPlot = subs(thetaSol,vars,values);

Постройте график движения гармонического маятника.

fplot(thetaSolPlot(t*T)/pi, [0 5]);
grid on;
title('Harmonic Pendulum Motion');
xlabel('t/T');
ylabel('\theta/\pi');

После нахождения решения для $$\theta (t)$$, визуализируйте движение маятника.

x_pos = sin(thetaSolPlot);
y_pos = -cos(thetaSolPlot);
fanimator(@fplot,x_pos,y_pos,'ko','MarkerFaceColor','k','AnimationRange',[0 5*T]);
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'),'AnimationRange',[0 5*T]);
fanimator(@(t) text(-0.3,0.3,"Timer: "+num2str(t,2)+" s"),'AnimationRange',[0 5*T]);

Введите команду playAnimation воспроизведение анимации движения маятника.

Шаг 6: Определите нелинейное движение маятника, используя постоянные энергетические пути

Чтобы понять нелинейное движение маятника, визуализируйте путь маятника, используя уравнение для общей энергии. Общая энергия сохранена.

Используйте тригонометрический тождества $$1-\cos \theta =2{\sin }^2(\theta /2)$$ и отношение $${\omega }_0=\sqrt{g/r}$$ чтобы переписать масштабированную энергию.

Поскольку энергия сохранена, движение маятника может быть описано постоянными энергетическими путями в фазу пространстве. Фазовое пространство является абстрактным пространством с координатами $$\theta$$ и $$d\theta /dt$$. Визуализируйте эти пути с помощью fcontour.

syms theta theta_t omega_0
E(theta, theta_t, omega_0) = (1/2)*(theta_t^2+(2*omega_0*sin(theta/2))^2);
Eplot(theta, theta_t) = subs(E,omega_0,omega_0Value);

figure;
fc = fcontour(Eplot(pi*theta, 2*omega_0Value*theta_t), 2*[-1 1 -1 1], ...
    'LineWidth', 2, 'LevelList', 0:5:50, 'MeshDensity', 1+2^8);
grid on;
title('Constant Energy Contours in Phase Space ( \theta vs. \theta_t )');
xlabel('\theta/\pi');
ylabel('\theta_t/2\omega_0');

Контуры постоянной энергии симметричны относительно $$\theta$$ ось и $$d\theta /dt$$ ось, и являются периодическими вдоль $$\theta$$ ось. Рисунок показывает две области различного поведения.

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

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

Шаг 7: Решение нелинейных уравнений движения

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

Переписать ОДУ второго порядка как систему ОДУ первого порядка.

syms theta(t) theta_t(t) omega_0
eqs = [diff(theta)   == theta_t;
       diff(theta_t) == -omega_0^2*sin(theta)]

eqs(t) = 
$$\left(\begin{array}{c}\frac{\partial }{\partial t}\mathrm{ }\theta \left(t\right)={\theta }_t\left(t\right)\\ \frac{\partial }{\partial t}\mathrm{ }{\theta }_t\left(t\right)=-{{\omega }_0}^2\hspace{0.17em}\sin \left(\theta \left(t\right)\right)\end{array}\right)$$

eqs  = subs(eqs,omega_0,omega_0Value);
vars = [theta, theta_t];

Найдите большую матрицу M системы и правых сторон уравнений F.

[M,F] = massMatrixForm(eqs,vars)

M = 
$$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

F = 
$$\left(\begin{array}{c}{\theta }_t\left(t\right)\\ -\frac{981\hspace{0.17em}\sin \left(\theta \left(t\right)\right)}{100}\end{array}\right)$$

M и F см. эту форму.

Чтобы упростить дальнейшие расчеты, перепишите систему в форме $$dx/dt=f(t,x(t))$$.

f = M\F

f = 
$$\left(\begin{array}{c}{\theta }_t\left(t\right)\\ -\frac{981\hspace{0.17em}\sin \left(\theta \left(t\right)\right)}{100}\end{array}\right)$$

Преобразование f в указатель на функцию MATLAB при помощи odeFunction. Получившийся указатель на функцию является входом к решателю MATLAB ODE ode45.

f = odeFunction(f, vars)

f = function_handle with value:
    @(t,in2)[in2(2,:);sin(in2(1,:)).*(-9.81e+2./1.0e+2)]

Шаг 8: Решение уравнения движения для замкнутых энергетических контуров

Решить ОДУ для замкнутых энергетических контуров можно используя ode45.

На графике энергетического контура замкнутые контуры удовлетворяют условию $${\theta }_0=0$$, $${\theta }_{t0}/2{\omega }_0\le 1$$. Сохраните начальные условия $$\theta$$ и $$d\theta /dt$$ в переменной x0.

x0 = [0; 1.99*omega_0Value];

Задайте временной интервал от 0 с до 10 с для поиска решения. Этот интервал позволяет маятнику пройти через два полных периода.

tInit  = 0;
tFinal = 10;

Решить ОДУ.

sols = ode45(f,[tInit tFinal],x0)

sols = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [1x45 double]
          y: [2x45 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

sols.y(1,:) представляет угловое перемещение $$\theta$$ и sols.y(2,:) представляет скорость вращения $$d\theta /dt$$.

Постройте график решения с замкнутой траекторией.

figure;

yyaxis left;
plot(sols.x, sols.y(1,:), '-o');
ylabel('\theta (rad)');

yyaxis right;
plot(sols.x, sols.y(2,:), '-o');
ylabel('\theta_t (rad/s)');

grid on;
title('Closed Path in Phase Space');
xlabel('t (s)');

Визуализируйте движение маятника.

x_pos = @(t) sin(deval(sols,t,1));
y_pos = @(t) -cos(deval(sols,t,1));
figure;
fanimator(@(t) plot(x_pos(t),y_pos(t),'ko','MarkerFaceColor','k'));
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'));
fanimator(@(t) text(-0.3,1.5,"Timer: "+num2str(t,2)+" s"));

Введите команду playAnimation воспроизведение анимации движения маятника.

Шаг 9: Решения по открытым энергетическим контурам

Решить ОДУ для открытых энергетических контуров можно используя ode45. На графике энергетического контура открытые контуры удовлетворяют условию $${\theta }_0=0$$, $${\theta }_{t0}/2{\omega }_0>1$$.

x0 = [0; 2.01*omega_0Value];
sols = ode45(f, [tInit, tFinal], x0);

Постройте график решения для открытого энергетического контура.

figure;

yyaxis left;
plot(sols.x, sols.y(1,:), '-o');
ylabel('\theta (rad)');

yyaxis right;
plot(sols.x, sols.y(2,:), '-o');
ylabel('\theta_t (rad/s)');

grid on;
title('Open Path in Phase Space');
xlabel('t (s)');

Визуализируйте движение маятника.

x_pos = @(t) sin(deval(sols,t,1));
y_pos = @(t) -cos(deval(sols,t,1));
figure;
fanimator(@(t) plot(x_pos(t),y_pos(t),'ko','MarkerFaceColor','k'));
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'));
fanimator(@(t) text(-0.3,1.5,"Timer: "+num2str(t,2)+" s"));

Введите команду playAnimation воспроизведение анимации движения маятника.