Trying to solve systems of ODEs and graph answer
1 次查看(过去 30 天)
显示 更早的评论
Trying to solve system of ODEs using other working code but getting error on line with "matlabFunction" I beleave the error has somethig to do with 'vars' but I can not figure out what. ODEs equation of motrions are correct and code works, error is just with graphing.
Original Code
clear;clc;
syms theta_1(t) theta_2(t) L_1 L_2 m_1 m_2 g
x_1 = L_1*sin(theta_1);
y_1 = -L_1*cos(theta_1);
x_2 = x_1 + L_2*sin(theta_2);
y_2 = y_1 - L_2*cos(theta_2);
vx_1 = diff(x_1);
vy_1 = diff(y_1);
vx_2 = diff(x_2);
vy_2 = diff(y_2);
ax_1 = diff(vx_1);
ay_1 = diff(vy_1);
ax_2 = diff(vx_2);
ay_2 = diff(vy_2);
syms T_1 T_2
eqx_1 = m_1*ax_1(t) == -T_1*sin(theta_1(t)) + T_2*sin(theta_2(t))
eqy_1 = m_1*ay_1(t) == T_1*cos(theta_1(t)) - T_2*cos(theta_2(t)) - m_1*g
eqx_2 = m_2*ax_2(t) == -T_2*sin(theta_2(t))
eqy_2 = m_2*ay_2(t) == T_2*cos(theta_2(t)) - m_2*g
Tension = solve([eqx_1 eqy_1],[T_1 T_2]);
eqRed_1 = subs(eqx_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);
eqRed_2 = subs(eqy_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);
L_1 = 1;
L_2 = 1.5;
m_1 = 2;
m_2 = 1;
g = 9.8;
eqn_1 = subs(eqRed_1)
eqn_2 = subs(eqRed_2)
[V,S] = odeToVectorField(eqn_1,eqn_2)
S
M = matlabFunction(V,'vars',{'t','Y'})
initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);
plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')
New codewith matlabFunction error
clear; clc;
% Initalize Variables
syms theta1(t) theta2(t) r1 r2 m1 m2 g
% Initalize Kinematic Constraints
x1 = r1*sin(theta1);
y1 = -r1*cos(theta1);
x2 = x1+r2*sin(theta2+theta1);
y2 = y1-r2*cos(theta2+theta1);
% Initalize Velocity Equations
x1dot=diff(x1);
y1dot=diff(y1);
x2dot=diff(x2);
y2dot=diff(y2);
% Solve for Potential Energy
pe=m1*g*y1+m2*g*y2;
%Solve for Kinetic Energy
ke=1/2*m1*(x1dot^2+y1dot^2)+1/2*m2*(x2dot^2+y2dot^2);
simplify(ke);
% Solve for L=KE-PE
l=ke-pe;
simplify(l);
% Initalize Variables
syms theta1dot theta2dot
% Solve for partial derivitive of L by theta1dot
eqn1=subs(diff(subs(l(t), diff(theta1(t)), theta1dot), theta1dot), theta1dot, diff(theta1(t)));
eqn1=simplify(eqn1);
% Solve for the time derivative of eqn1
eqn1=diff(eqn1);
eqn1=simplify(eqn1);
% Solve for partial derivitive of L by theta1
a=simplify(diff(l,theta1));
% Solve for first equation of motion
eqn1=simplify(eqn1-a)
% Solve for partial derivitive of L by theta2dot
eqn2=subs(diff(subs(l(t), diff(theta2(t)), theta2dot), theta2dot), theta2dot, diff(theta2(t)));
eqn2=simplify(eqn2);
% Solve for the time derivative of eqn2
eqn2=diff(eqn2);
eqn2=simplify(eqn2);
% Solve for partial derivitive of L by theta2
b=simplify(diff(l,theta2));
% Solve for second equation of motion
eqn2=simplify(eqn2-b)
r1 = 1;
r2 = 1.5;
m1 = 2;
m2 = 1;
g = 9.8;
[V,S] = odeToVectorField(eqn1,eqn2)
M = matlabFunction(V,'vars',{'t','Y'})
initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);
plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')
0 个评论
回答(1 个)
Walter Roberson
2022-12-1
If you define a symbolic variable and use it in an expression, and later assign a numeric value to the same name, then that does not change the symbolic variable that is in the expression. The situation is much like how if you
A = 3
B = A*2 + 1
A = 9
then B does not suddenly change from 3*2+1 to instead evaluate to 9*2+1 . The value of A as of the time that B was defined is what is used in B, and changes to A after that do not affect B. The same is true for if it had been syms A instead of assigning a numeric value to A -- changing A afterwards would not change the place the symbolic version was already used.
% Initalize Variables
syms theta1(t) theta2(t) r1 r2 m1 m2 g
% Initalize Kinematic Constraints
x1 = r1*sin(theta1);
y1 = -r1*cos(theta1);
x2 = x1+r2*sin(theta2+theta1);
y2 = y1-r2*cos(theta2+theta1);
% Initalize Velocity Equations
x1dot=diff(x1);
y1dot=diff(y1);
x2dot=diff(x2);
y2dot=diff(y2);
% Solve for Potential Energy
pe=m1*g*y1+m2*g*y2;
%Solve for Kinetic Energy
ke=1/2*m1*(x1dot^2+y1dot^2)+1/2*m2*(x2dot^2+y2dot^2);
simplify(ke);
% Solve for L=KE-PE
l=ke-pe;
simplify(l);
% Initalize Variables
syms theta1dot theta2dot
% Solve for partial derivitive of L by theta1dot
eqn1=subs(diff(subs(l(t), diff(theta1(t)), theta1dot), theta1dot), theta1dot, diff(theta1(t)));
eqn1=simplify(eqn1);
% Solve for the time derivative of eqn1
eqn1=diff(eqn1);
eqn1=simplify(eqn1);
% Solve for partial derivitive of L by theta1
a=simplify(diff(l,theta1));
% Solve for first equation of motion
eqn1=simplify(eqn1-a)
% Solve for partial derivitive of L by theta2dot
eqn2=subs(diff(subs(l(t), diff(theta2(t)), theta2dot), theta2dot), theta2dot, diff(theta2(t)));
eqn2=simplify(eqn2);
% Solve for the time derivative of eqn2
eqn2=diff(eqn2);
eqn2=simplify(eqn2);
% Solve for partial derivitive of L by theta2
b=simplify(diff(l,theta2));
% Solve for second equation of motion
eqn2=simplify(eqn2-b)
r1 = 1;
r2 = 1.5;
m1 = 2;
m2 = 1;
g = 9.8;
[V,S] = odeToVectorField(subs(eqn1),subs(eqn2))
M = matlabFunction(V,'vars',{'t','Y'})
initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);
plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!