Симулируйте движение периодического Swing маятника

В этом примере показано, как симулировать движение математического маятника с помощью Symbolic Math Toolbox™. Выведите уравнение движения маятника, затем решите уравнение аналитически для маленьких углов и численно для любого угла.

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

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

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

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

a=rd2θdt2.

Займите место a при помощи subs.

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

mr2t2 θ(t)=-gmsin(θ(t))m*r*diff (тета (t), t, 2) ==-g*m*sin (тета (t))

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

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

2t2 θ(t)=-gsin(θ(t))rdiff (тета (t), t, 2) == - (g*sin (тета (t)))/r

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

ω0=gr.

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

2t2 θ(t)=-ω02sin(θ(t))diff (тета (t), t, 2) ==-omega_0^2*sin (тета (t))

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

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

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

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

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

2t2 θ(t)=-ω02θ(t)diff (тета (t), t, 2) ==-omega_0^2*theta (t)

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

Решите уравнение eqnLinear при помощи dsolve. Задайте начальные условия в качестве второго аргумента. Упростите решение путем принятия ω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) = 

θ0cos(ω0t)+θt0sin(ω0t)ω0theta_0*cos (omega_0*t) + (theta_t0*sin (omega_0*t))/omega_0

Шаг 4: физическое значение ω0

Термин ω0t называется фазой. Косинусные и синусоидальные функции повторяют каждый 2π. Время должно было измениться ω0t 2π называется периодом времени.

T=2πω0=2πrg.

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

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

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

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

  • Гравитационное ускорение g=9.81m/s2

  • Длина маятника r=1m

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');

После нахождения решения для θ(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: определите нелинейное движение маятника Используя постоянные энергетические пути

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

E=12mr2(dθdt)2+mgr(1-cosθ)

Используйте тригонометрическую идентичность 1-cosθ=2sin2(θ/2) и отношение ω0=g/r переписать масштабированную энергию.

Emr2=12[(dθdt)2+(2ω0sinθ2)2]

Поскольку энергия сохраняется, движение маятника может быть описано постоянными энергетическими путями в фазовом пространстве. Фазовое пространство является абстрактным пробелом с координатами θ и dθ/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');

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

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

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

Шаг 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) = 

(t θ(t)=θt(t)t θt(t)=-ω02sin(θ(t)))[diff (тета (t), t) == theta_t (t); diff (theta_t (t), t) ==-omega_0^2*sin (тета (t))]

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

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

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

(1001)[sym (1), sym (0); sym (0), sym (1)]

F = 

(θt(t)-981sin(θ(t))100)[theta_t (t); - (981*sin (тета (t)))/100]

M и F обратитесь к этой форме.

M(t,x(t))dxdt=F(t,x(t)).

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

f = M\F
f = 

(θt(t)-981sin(θ(t))100)[theta_t (t); - (981*sin (тета (t)))/100]

Преобразуйте 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.

От энергетического контурного графика замкнутые контуры удовлетворяют условию θ0=0, θt0/2ω01. Сохраните начальные условия θ и dθ/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,:) представляет угловое смещение θ и sols.y(2,:) представляет скорость вращения 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. От энергетического контурного графика открытые контуры удовлетворяют условию θ0=0, θt0/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 проигрывать анимацию движения маятника.

Для просмотра документации необходимо авторизоваться на сайте