Шаг 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 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 проигрывать анимацию движения маятника.