Документация

Оглавление

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)$$

Ускорение массы в конце маятника

Замените 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$$ в один параметр, названный собственной частотой,

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. Задайте начальные условия в качестве второго аргумента.

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

thetaSol(t) = 
$$\frac{0.5000\hspace{0.17em}{e}^{-1\hspace{0.17em}{\omega }_0\hspace{0.17em}t\hspace{0.17em}i}\hspace{0.17em}\left({\omega }_0\hspace{0.17em}{\theta }_0+1\hspace{0.17em}{\theta }_{\mathrm{t0}}\hspace{0.17em}i\right)}{{\omega }_0}-\frac{0.5000\hspace{0.17em}{e}^{1\hspace{0.17em}{\omega }_0\hspace{0.17em}t\hspace{0.17em}i}\hspace{0.17em}\left(-{\omega }_0\hspace{0.17em}{\theta }_0+1\hspace{0.17em}{\theta }_{\mathrm{t0}}\hspace{0.17em}i\right)}{{\omega }_0}$$

Упростите решение путем принятия $${\omega }_0$$ действительное использование assume.

assume(omega_0,'real')
thetaSol(t) = dsolve(eqnLinear,cond)

thetaSol(t) = 
$${\theta }_0\hspace{0.17em}\mathrm{потому\; что}\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$$ называется фазой. Косинусная функция $\mathrm{потому\; что}{\omega }_0\mathit{t}$ повторяет каждый $$2\pi$$. Время должно было измениться $${\omega }_{0}t$$ $$2\pi$$ называется периодом времени

От уравнения период времени прямо пропорционален к длине маятника. Однако $$t_0$$ не зависит от массы, потому что ее момент инерции и ее вес оба прямо пропорциональны к его массе.

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

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

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

Задайте физические параметры.

gValue = 9.81; % m / s^2
rValue = 1;    % m
omega_0Value = sqrt(gValue/rValue);
t_0 = 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);

figure;
fplot(thetaSolPlot(t*t_0)/pi, [0 1]);

fplot(thetaSolPlot(t*t_0));
grid on;
title('Harmonic Pendulum Motion');
xlabel('t/t_0');
ylabel('\theta/\pi');

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

6. Понимание нелинейного движения маятника Используя постоянные энергетические пути

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

Полная энергия

Мы можем использовать тригонометрическую идентичность

Используйте отношение $${\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)\\ -9.8100\hspace{0.17em}\sin \left(\theta \left(t\right)\right)\end{array}\right)$$

M и F относятся к форме

Чтобы упростить дальнейшие вычисления, перепишите систему в форме

f = M\F

f = 
$$\left(\begin{array}{c}{\theta }_t\left(t\right)\\ -9.8100\hspace{0.17em}\sin \left(\theta \left(t\right)\right)\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; 2*omega_0Value];

Выберите временной интервал, который позволяет маятнику проходить полный период. Это может быть найдено методом проб и ошибок.

tInit  = 0;
tFinal = 3.365*t_0;

Решите ОДУ.

[t, x] = ode45(f, [tInit tFinal], x0);

$$\theta$$ возвращен в первом столбце x, и $$d\theta /dt$$ возвращен во втором столбце x.

Масштабируйте измерение времени периодом $$t_0$$,

t  = t/t_0;

Сделайте значения $$\theta$$ повторитесь на интервале $$(-\pi ,\pi )$$ с помощью wrapToPi затем масштабируйтесь $$\theta$$ $$\pi$$.

x(:,1) = wrapToPi(x(:,1))/pi;

Шкала $$d\theta /dt$$ $$2{\omega }_0$$.

x(:,2) = x(:,2)/(2*omega_0Value);

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

figure;

yyaxis left;
plot(t, x(:,1), '-o');
ylabel('\theta/\pi');

yyaxis right;
plot(t, x(:,2), '-o');
ylabel('\theta_t/2\omega_0');

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

Это - то, как движение маятника появилось бы:

9. Решения на открытых энергетических контурах

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

От энергетического контурного графика открытые контуры удовлетворяют условие $${\theta }_0=0$$, $${\theta }_{t0}/2{\omega }_0>1$$.

x0 = [0; 1.001*2*omega_0Value];
tFinal = 1.465*t_0;

[t, x] = ode45(f, [tInit, tFinal], x0);

t = t/t_0;
x(:,1) = wrapToPi(x(:,1))/pi;
x(:,2) = x(:,2)/(2*omega_0Value);

figure;

yyaxis left;
plot(t, x(:,1), '-o');
ylabel('\theta/\pi');

yyaxis right;
plot(t, x(:,2), '-o');
ylabel('\theta_t/2\omega_0');

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

Это - то, как движение маятника появилось бы:

10. Заключение

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

Функции помощника

function lambda = wrapToPi(lambda)
%wrapToPi Wrap angle in radians to [-pi pi]
%
%   lambdaWrapped = wrapToPi(LAMBDA) wraps angles in LAMBDA, in radians,
%   to the interval [-pi pi] such that pi maps to pi and -pi maps to
%   -pi.  (In general, odd, positive multiples of pi map to pi and odd,
%   negative multiples of pi map to -pi.)

q = (lambda < -pi) | (pi < lambda);
lambda(q) = wrapTo2Pi(lambda(q) + pi) - pi;
end

function lambda = wrapTo2Pi(lambda)
%wrapTo2Pi Wrap angle in radians to [0 2*pi]
%
%   lambdaWrapped = wrapTo2Pi(LAMBDA) wraps angles in LAMBDA, in radians,
%   to the interval [0 2*pi] such that zero maps to zero and 2*pi maps
%   to 2*pi. (In general, positive multiples of 2*pi map to 2*pi and
%   negative multiples of 2*pi map to zero.)

positiveInput = (lambda > 0);
lambda = mod(lambda, 2*pi);
lambda((lambda == 0) & positiveInput) = 2*pi;
end