Why did i receive error with my for loop ode45 function?
4 views (last 30 days)
Show older comments
Hi matlabbers, I have been trying to figure out how to find thetadot2 and plot it on a thetadot1 vs time graph.
For context, The main issue is with the ode45 function. I did the dynamics for a double pendelum system symbolically and now im trying to convert it into numbers for the ode45 function but its not working. Any idea on how to fix this? I was thinking that in order for it to work, i only need one initial condition, in my case thetadot2.
Red text error message.
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
@(T,Y)[Y(2);-((SUBS(DL_DTHETA2,THETADOT1,THETADOT1_ITERATION)-SUBS(DL_DTHETADOT2,THETADOT1,THETADOT1_ITERATION)+Q)/JBLADE)] returns a vector of length 1, but the length of initial conditions vector is 2. the initial conditions vector must have the same number of elements.
Code
%{Assuming we have Parameters for
%L1, L2, m, Jblade, kt, M_Preload, Jblade, time_intervals,
%thetadot1_invervals, thetadot1&2, M_preload, M_centrifugal
%initial conditions
dL_dtheta2_numeric = subs(dL_dtheta2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade}); % Perform symbolic substitutions to get numeric versions of dL_dtheta2 and dt_dL_dthetadot2
dt_dL_dthetadot2_numeric = subs(dt_dL_dthetadot2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade});
theta2_solution_values_total = []; % Loop through each time interval and store thetadot2values
time_solution_values_total = [];
for i = 1:length(time_intervals) - 1
thetadot1_i_value = thetadot1_rad(i); % Get the current value of thetadot1 in radians
Q = M_Preload - M_centrifugal(thetadot1_i_value); % Compute the generalized force term Q
odeFunction = @(t, y) [y(2); -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric + Q) / Jblade]; % Define the ODE function with fully numeric values
tspan = [time_intervals(i), time_intervals(i + 1)]; % Time span for the current interval
[t_interval, sol_interval] = ode45(odeFunction, tspan, [theta2_initial; thetadot2_initial]); % Solve the ODE
theta2_initial = sol_interval(end, 1);% Update initial conditions for the next interval
thetadot2_initial = sol_interval(end, 2);
theta2_solution_values_total = [theta2_solution_values_total; sol_interval(:, 1)]; % Store solutions
time_solution_values_total = [time_solution_values_total; t_interval];
end
3 Comments
Walter Roberson
on 4 Nov 2024
dt_dL_dthetadot2_numeric = subs(dt_dL_dthetadot2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade});
That creates a symbolic value
odeFunction = @(t, y) [y(2); -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric + Q) / Jblade]; % Define the ODE function with fully numeric values
That creates an anonymous function that returns a symbolic expression
[t_interval, sol_interval] = ode45(odeFunction, tspan, [theta2_initial; thetadot2_initial]); % Solve the ODE
ode45 cannot use return values that are symbolic expressions.
Answers (2)
Walter Roberson
on 4 Nov 2024
Moved: Walter Roberson
on 4 Nov 2024
% Goal: Solve EOM using Lagrange for theta2 and plot theta2 vs time
% Clear old results
close all
clc
clear all
% Define symbols for DOF theta variables and parameters
syms L1 L2 mblade kt Jblade theta1(t) theta2(t)
% Define derivatives of theta
thetadot1 = diff(theta1, t);
thetadot2 = diff(theta2, t);
% Position of the second joint
x2 = L1*cos(theta1) + L2*cos(theta1 + theta2);
y2 = L1*sin(theta1) + L2*sin(theta1 + theta2);
% Velocities of the second joint
xdot2 = diff(x2, t);
ydot2 = diff(y2, t);
% Kinetic and Potential Energy
KE_Translation = (1/2)*mblade*(xdot2^2 + ydot2^2);
KE_Rotation = (1/2)*Jblade*(thetadot1 + thetadot2)^2;
PETorSpring = (1/2)*kt*theta2^2;
T = KE_Translation + KE_Rotation;
V = PETorSpring;
L = T - V; % Lagrangian
% Derivatives of the Lagrangian
dL_dtheta1 = diff(L, theta1);
dL_dtheta2 = diff(L, theta2);
dL_dthetadot1 = diff(L, thetadot1);
dL_dthetadot2 = diff(L, thetadot2);
dt_dL_dthetadot1 = diff(dL_dthetadot1, t);
dt_dL_dthetadot2 = diff(dL_dthetadot2, t);
% Define parameters (numeric values)
L1 = 0.60; % Length of first arm (m)
L2 = 0.60; % Length of second arm (m)
mblade = 2; % Mass of one blade (kg)
kt = 0.3; % Spring constant (Nm/rad)
M_Preload = 5; % Spring preload (Nm)
M_centrifugal = @(thetadot1) mblade * (L1 + L2 / 2) * thetadot1^2 * (L1 + L2 / 2);
Jblade = (1/12 * mblade * L2^2) + (mblade * (L1 + L2/2)^2); % ICR Moment of inertia of the blade
time_intervals = [0, 30, 60, 90, 120];
thetadot1_rpm = [0, 50, 100, 150, 0];
thetadot1_rad = thetadot1_rpm * (2 * pi / 60); % Convert RPM to rad/s
theta2_initial = 0
thetadot2_initial = 0
% Define the function for centrifugal moment
% Perform symbolic substitutions to get numeric versions of dL_dtheta2 and dt_dL_dthetadot2
dL_dtheta2_numeric = subs(dL_dtheta2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade});
dt_dL_dthetadot2_numeric = subs(dt_dL_dthetadot2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade});
% Loop through each time interval and thetadot1 values
theta2_solution_values_total = [];
time_solution_values_total = [];
for i = 1:length(time_intervals) - 1
% Get the current value of thetadot1 in radians
thetadot1_i_value = thetadot1_rad(i);
% Compute the generalized force term Q
Q = M_Preload - M_centrifugal(thetadot1_i_value);
% Define the ODE function with fully numeric values
temp = -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric + Q) / Jblade
whos temp
odeFunction = @(t, y) [y(2); -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric + Q) / Jblade];
% Time span for the current interval
tspan = [time_intervals(i), time_intervals(i + 1)];
% Solve the ODE
[t_interval, sol_interval] = ode45(odeFunction, tspan, [theta2_initial; thetadot2_initial]);
% Update initial conditions for the next interval
theta2_initial = sol_interval(end, 1);
thetadot2_initial = sol_interval(end, 2);
% Store solutions
theta2_solution_values_total = [theta2_solution_values_total; sol_interval(:, 1)];
time_solution_values_total = [time_solution_values_total; t_interval];
end
% Plot results
plot(time_solution_values_total, theta2_solution_values_total);
xlabel('Time (s)');
ylabel('\theta_2 (rad)');
title('\theta_2 vs Time with Generalized Forces');
grid on;
Your expression -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric + Q) / Jblade is a symfun. When [] together with something else, the result is a single symfun that returns the concatenation of the values. So your ode function is returning a single symfun, and ode45 is complaining about that.
1 Comment
Walter Roberson
on 4 Nov 2024
Your error is mixing symbolic values with numeric values. ode45() cannot handle symbolic results. You should be using matlabFunction or odeFunction to create your function handle.
Torsten
on 5 Nov 2024
Edited: Torsten
on 5 Nov 2024
I used numerical values for L1, L2, mblade, kt and Jblade right from the beginning and got expressions for dL_dtheta2_numeric and dt_dL_dthetadot2_numeric that depend on theta1, theta1', theta1'', theta2, theta2' and theta2''.
Therefore, I expected two second-order differential equations for theta1 and theta2, but cannot find them. Could you write down the mathematical system of ODEs you are trying to solve ?
% Goal: Solve EOM using Lagrange for theta2 and plot theta2 vs time
% Clear old results
close all
clc
clear all
% Define symbols for DOF theta variables and parameters
syms theta1(t) theta2(t)
% Define parameters (numeric values)
L1 = 0.60; % Length of first arm (m)
L2 = 0.60; % Length of second arm (m)
mblade = 2; % Mass of one blade (kg)
kt = 0.3; % Spring constant (Nm/rad)
M_Preload = 5; % Spring preload (Nm)
Jblade = (1/12 * mblade * L2^2) + (mblade * (L1 + L2/2)^2); % ICR Moment of inertia of the blade
% Define derivatives of theta
thetadot1 = diff(theta1, t);
thetadot2 = diff(theta2, t);
% Position of the second joint
x2 = L1*cos(theta1) + L2*cos(theta1 + theta2);
y2 = L1*sin(theta1) + L2*sin(theta1 + theta2);
% Velocities of the second joint
xdot2 = diff(x2, t);
ydot2 = diff(y2, t);
% Kinetic and Potential Energy
KE_Translation = (1/2)*mblade*(xdot2^2 + ydot2^2);
KE_Rotation = (1/2)*Jblade*(thetadot1 + thetadot2)^2;
PETorSpring = (1/2)*kt*theta2^2;
T = KE_Translation + KE_Rotation;
V = PETorSpring;
L = T - V; % Lagrangian
% Derivatives of the Lagrangian
dL_dtheta1 = diff(L, theta1);
dL_dtheta2 = diff(L, theta2);
dL_dthetadot1 = diff(L, thetadot1);
dL_dthetadot2 = diff(L, thetadot2);
dt_dL_dthetadot1 = diff(dL_dthetadot1, t);
dt_dL_dthetadot2 = diff(dL_dthetadot2, t);
M_centrifugal = @(thetadot1) mblade * (L1 + L2 / 2) * thetadot1^2 * (L1 + L2 / 2);
time_intervals = [0, 30, 60, 90, 120];
thetadot1_rpm = [0, 50, 100, 150, 0];
thetadot1_rad = thetadot1_rpm * (2 * pi / 60); % Convert RPM to rad/s
theta2_initial = 0
thetadot2_initial = 0
% Define the function for centrifugal moment
% Perform symbolic substitutions to get numeric versions of dL_dtheta2 and dt_dL_dthetadot2
dL_dtheta2_numeric = subs(dL_dtheta2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade})
dt_dL_dthetadot2_numeric = subs(dt_dL_dthetadot2, {L1, L2, mblade, kt, Jblade}, {L1, L2, mblade, kt, Jblade})
% Loop through each time interval and thetadot1 values
theta2_solution_values_total = [];
time_solution_values_total = [];
for i = 1:length(time_intervals) - 1
% Get the current value of thetadot1 in radians
thetadot1_i_value = thetadot1_rad(i);
% Compute the generalized force term Q
Q = M_Preload - M_centrifugal(thetadot1_i_value);
% Define the ODE function with fully numeric values
odeFunction = @(t, y) [y(2); -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric + Q) / Jblade];
% Time span for the current interval
tspan = [time_intervals(i), time_intervals(i + 1)];
% Solve the ODE
[t_interval, sol_interval] = ode45(odeFunction, tspan, [theta2_initial; thetadot2_initial]);
% Update initial conditions for the next interval
theta2_initial = sol_interval(end, 1);
thetadot2_initial = sol_interval(end, 2);
% Store solutions
theta2_solution_values_total = [theta2_solution_values_total; sol_interval(:, 1)];
time_solution_values_total = [time_solution_values_total; t_interval];
end
% Plot results
plot(time_solution_values_total, theta2_solution_values_total);
xlabel('Time (s)');
ylabel('\theta_2 (rad)');
title('\theta_2 vs Time with Generalized Forces');
grid on;
10 Comments
Torsten
on 7 Nov 2024
Edited: Torsten
on 8 Nov 2024
"solve" is used to get the equations that "ode45" has to compute afterwards. Thus "solve" is not a substitute for "ode45", but used as a preprocessing step.
Here is the complete solution with ode45.
Note that I used theta1 = 2, theta1dot = 0, theta2 = 0, thea2dot = 0 at t = 0 as initial values and worked without external forces. "odeToVectorField" does what "solve" did in the other code.
% Goal: Solve EOM using Lagrange for theta2 and plot theta2 vs time
% Clear old results
close all
clc
clear all
% Define symbols for DOF theta variables and parameters
syms theta1(t) theta2(t)
% Define parameters (numeric values)
L1 = 0.60; % Length of first arm (m)
L2 = 0.60; % Length of second arm (m)
mblade = 2; % Mass of one blade (kg)
kt = 0.3; % Spring constant (Nm/rad)
M_Preload = 5; % Spring preload (Nm)
Jblade = (1/12 * mblade * L2^2) + (mblade * (L1 + L2/2)^2); % ICR Moment of inertia of the blade
% Define derivatives of theta
thetadot1 = diff(theta1, t);
thetadot2 = diff(theta2, t);
% Position of the second joint
x2 = L1*cos(theta1) + L2*cos(theta1 + theta2);
y2 = L1*sin(theta1) + L2*sin(theta1 + theta2);
% Velocities of the second joint
xdot2 = diff(x2, t);
ydot2 = diff(y2, t);
% Kinetic and Potential Energy
KE_Translation = (1/2)*mblade*(xdot2^2 + ydot2^2);
KE_Rotation = (1/2)*Jblade*(thetadot1 + thetadot2)^2;
PETorSpring = (1/2)*kt*theta2^2;
T = KE_Translation + KE_Rotation;
V = PETorSpring;
L = T - V; % Lagrangian
% Derivatives of the Lagrangian
dL_dtheta1 = diff(L, theta1);
dL_dtheta2 = diff(L, theta2);
dL_dthetadot1 = diff(L, thetadot1);
dL_dthetadot2 = diff(L, thetadot2);
dt_dL_dthetadot1 = diff(dL_dthetadot1, t);
dt_dL_dthetadot2 = diff(dL_dthetadot2, t);
% Euler-Lagrange equations for theta1 and theta2
eq1 = simplify(dt_dL_dthetadot1 - dL_dtheta1) == 0; % Eq of motion for theta1
eq2 = simplify(dt_dL_dthetadot2 - dL_dtheta2) == 0;
V = odeToVectorField([eq1,eq2]);
M = matlabFunction(V,'vars',{'t','Y'});
interval = [0 20];
yInit = [2 0 0 0]; %[theta1,theta1dot,theta2,theta2dot] at t = 0
ySol = ode45(M,interval,yInit);
tValues = linspace(0,20,100);
yValues = deval(ySol,tValues);
plot(tValues,yValues(1,:))
See Also
Categories
Find more on Code Generation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!