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

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

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=\frac{{}^{\mathrm{}}}{{}^{\mathrm{}}}$$rd2startdt2.

Заменить $$\mathrm{на\; с}$$ помощью 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$$ в один параметр, который также известен как собственная частота.

${}_{\mathrm{\lambda \; 0}}\sqrt{\frac{\text{=}\mathit{}}{\mathit{}}}$гр.

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: Линеаризация уравнения движения

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

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. Укажите начальные условия в качестве второго аргумента. Упростите решение, предположив, что $${}_{\mathrm{group0}}$$ является реальным, используя 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: Физическая значимость$${}_{\mathrm{}}$$

Термин $${}_{\mathrm{}}\mathrm{start0t}$$ называется фазой. Косинусная и синусоидальная функции повторяются каждые $$\mathrm{}\mathrm{2\delta }$$. Время, необходимое для изменения $${}_{\mathrm{}}\mathrm{start0t}$$ на $$\mathrm{}\mathrm{2\delta }$$, называется периодом времени.

$$Т\frac{\mathrm{=}}{{}_{\mathrm{}}}\mathrm{2íà\; 0}\sqrt{\frac{=}{}}$$2.drg.

Период времени $$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');

Найдя решение для ($$\mathrm{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: Определение нелинейного движения маятника с использованием путей постоянной энергии

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

$$\mathrm{}\mathrm{}\mathrm{}{\mathrm{}}^{\mathrm{}}\mathrm{}$$$${}_{\mathrm{}}\sqrt{}$$Чтобы переписать масштабированную энергию, используйте тригонометрическое тождество 1-cosstart= 2sin2 (start/2) и уравнение λ 0 = g/r.

Поскольку энергия сохраняется, движение маятника может быть описано путями постоянной энергии в фазовом пространстве. Фазовое пространство представляет собой абстрактное пространство с координатами, $$и$$ $$d$$Визуализация этих путей с помощью 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$$ оси и $$\mathrm{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 см. эту форму.

Для упрощения дальнейших вычислений переписать систему в виде $$\mathrm{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.

Из графика энергетических контуров замкнутые контуры удовлетворяют условию $${}_{\mathrm{start0}}\mathrm{=}$$0$${,}_{\mathrm{}}\mathrm{}{}_{\mathrm{}}\mathrm{}$$θt0/2ω0≤1. Запишите в $$$$$$$$ переменную начальные условия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,:) представляет угловое смещение$$$$ sols.y(2,:) представляет собой угловую скорость $$\mathrm{d/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. Из графика энергетических контуров, открытые контуры удовлетворяют условию $${}_{\mathrm{start0}}\mathrm{=}$$0$${,}_{\mathrm{}}\mathrm{}{}_{\mathrm{}}\mathrm{startt0/2ü\; 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 для воспроизведения анимации движения маятника.