В этом примере показано, как симулировать движение математического маятника с помощью Symbolic Math Toolbox™. Выведите уравнение движения маятника, затем решите уравнение аналитически для маленьких углов и численно для любого угла.
Маятник является простой механической системой, которая следует за дифференциальным уравнением. Маятник находится первоначально в покое в вертикальном положении. Когда маятник перемещен углом и выпущенный, сила тяжести задерживает его к своей позиции отдыха. Его импульс заставляет его промахиваться и прибывать в угол (при отсутствии фрикционных сил), и так далее. Сила восстановления вдоль движения маятника из-за силы тяжести . Таким образом, согласно второму закону Ньютона, массовые времена ускорение должно равняться .
syms m a g theta(t) eqn = m*a == -m*g*sin(theta)
eqn(t) =
Для маятника с длиной , ускорение боба маятника равно угловым ускоряющим временам .
.
Займите место при помощи subs
.
syms r
eqn = subs(eqn,a,r*diff(theta,2))
eqn(t) =
Изолируйте угловое ускорение в eqn
при помощи isolate
.
eqn = isolate(eqn,diff(theta,2))
eqn =
Соберите константы и в один параметр, который также известен как собственную частоту.
.
syms omega_0
eqn = subs(eqn,g/r,omega_0^2)
eqn =
Уравнение движения нелинейно, таким образом, это затрудняет, чтобы решить аналитически. Примите, что углы малы и линеаризуют уравнение при помощи Разложения Тейлора .
syms x approx = taylor(sin(x),x,'Order',2); approx = subs(approx,x,theta(t))
approx =
Уравнение движения становится линейным уравнением.
eqnLinear = subs(eqn,sin(theta(t)),approx)
eqnLinear =
Решите уравнение eqnLinear
при помощи dsolve
. Задайте начальные условия в качестве второго аргумента. Упростите решение путем принятия действительное использование 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) =
Термин называется фазой. Косинусные и синусоидальные функции повторяют каждый . Время должно было измениться называется периодом времени.
.
Период времени пропорционально квадратному корню из длины маятника, и это не зависит от массы. Для линейного уравнения движения период времени не зависит от начальных условий.
Постройте движение маятника для малого углового приближения.
Задайте физические параметры:
Гравитационное ускорение
Длина маятника
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');
После нахождения решения для , визуализируйте движение маятника.
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
проигрывать анимацию движения маятника.
Чтобы изучить нелинейное движение маятника, визуализируйте путь к маятнику при помощи уравнения для полной энергии. Полная энергия сохраняется.
Используйте тригонометрическую идентичность и отношение переписать масштабированную энергию.
Поскольку энергия сохраняется, движение маятника может быть описано постоянными энергетическими путями в фазовом пространстве. Фазовое пространство является абстрактным пробелом с координатами и . Визуализируйте эти пути с помощью 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');
Постоянные энергетические контуры симметричны о ось и ось, и является периодической вперед ось. Рисунок показывает две области отличного поведения.
Более низкие энергии контурного графика закрываются на себя. Маятник качается назад и вперед между двумя максимальными углами и скоростями. Кинетической энергии маятника недостаточно, чтобы преодолеть гравитационную энергию и позволить маятнику сделать полный цикл.
Более высокие энергии контурного графика не закрываются на себя. Маятник всегда перемещается в одно угловое направление. Кинетической энергии маятника достаточно, чтобы преодолеть гравитационную энергию и позволить маятнику сделать полный цикл.
Нелинейные уравнения движения являются дифференциальными уравнениями второго порядка. Численно решите эти уравнения при помощи 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) =
eqs = subs(eqs,omega_0,omega_0Value); vars = [theta, theta_t];
Найдите большую матрицу M
из системы и правых сторон уравнений F
.
[M,F] = massMatrixForm(eqs,vars)
M =
F =
M
и F
обратитесь к этой форме.
Чтобы упростить дальнейшие расчеты, перепишите систему в форме .
f = M\F
f =
Преобразуйте 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)]
Решите ОДУ для закрытых энергетических контуров при помощи ode45
.
От энергетического контурного графика замкнутые контуры удовлетворяют условию , . Сохраните начальные условия и в переменной 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,:)
представляет скорость вращения .
Постройте решение замкнутого пути.
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
проигрывать анимацию движения маятника.
Решите ОДУ для открытых энергетических контуров при помощи ode45
. От энергетического контурного графика открытые контуры удовлетворяют условию , .
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
проигрывать анимацию движения маятника.