Why did i receive error with my for loop ode45 function?

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.

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error using odearguments (line 98)

@(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

[Torsten]

I don't see the line for which the error is reported in your code:

@(T,Y)[Y(2);-((SUBS(DL_DTHETA2,THETADOT1,THETADOT1_ITERATION)-SUBS(DL_DTHETADOT2,THETADOT1,THETADOT1_ITERATION)+Q)/JBLADE)]

Please supply executable code that reproduces the error you get.

[Walter Roberson]

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.

You need to use matlabFunction or odeFunction to generate the function handle.

[Miles]

Hi Torsten, this is the full code. I wanted to only put the for loop to not overwhelm anyone. If you run this code, there should be the error as shown above with ode45

Full 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 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

theta2_initial = 0

thetadot2_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

Error using odearguments (line 98)
@(T,Y)[Y(2);-(DL_DTHETA2_NUMERIC-DT_DL_DTHETADOT2_NUMERIC+Q)/JBLADE] returns a vector of length 1, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[Y(2);-(DL_DTHETA2_NUMERIC-DT_DL_DTHETADOT2_NUMERIC+Q)/JBLADE] and the initial conditions vector must have the same number of elements.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% 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;

Answer 1

[Walter Roberson]

% 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

theta2_initial = 0

thetadot2_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

  Name      Size            Bytes  Class     Attributes

  temp      1x1                 8  symfun              

Error using odearguments (line 98)
@(T,Y)[Y(2);-(DL_DTHETA2_NUMERIC-DT_DL_DTHETADOT2_NUMERIC+Q)/JBLADE] returns a vector of length 1, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[Y(2);-(DL_DTHETA2_NUMERIC-DT_DL_DTHETADOT2_NUMERIC+Q)/JBLADE] and the initial conditions vector must have the same number of elements.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% 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.

[Walter Roberson]

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.

Answer 2

[Torsten]

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

theta2_initial = 0

thetadot2_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

Error using odearguments (line 98)
@(T,Y)[Y(2);-(DL_DTHETA2_NUMERIC-DT_DL_DTHETADOT2_NUMERIC+Q)/JBLADE] returns a vector of length 1, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[Y(2);-(DL_DTHETA2_NUMERIC-DT_DL_DTHETADOT2_NUMERIC+Q)/JBLADE] and the initial conditions vector must have the same number of elements.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% 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;

[Walter Roberson]

This code continues to have the problem that odeFunction returns a scalar symfun, and symfun are of length 1 which is the wrong size for ode45.

[Torsten]

@Walter Roberson

This code continues to have the problem that odeFunction returns a scalar symfun, and symfun are of length 1 which is the wrong size for ode45.

This is a technical problem. The main problem is to get a proper mathematical formulation of the ODEs.

[Miles]

So only diff equation i should have is EOM2 since i already am given values for theta1dot. We get 1 second order diff equation from lagrange, solve them on matlab and then I only plot the theta 2 equation. I made new assumtions that should make solving this easier being that the lagrange equation equals to zero and preload force is defined and is accounted for in the potential energy. I also tried using matlab function but now get this error.

Array indices must be positive integers or logical values.

Error in sym/matlabFunction>@(t,theta2_initial,theta2,thetadot2,theta1,thetadot1,L1,L2,mblade,kt,Jblade)-kt.*theta2_initial+kt.*theta2(t)+Jblade.*diff(theta1(t),t,t)+Jblade.*diff(theta2(t),t,t)+L2.^2.*mblade.*diff(theta1(t),t,t)+L2.^2.*mblade.*diff(theta2(t),t,t)+L1.*L2.*mblade.*sin(theta2(t)).*diff(theta1(t),t).^2+L1.*L2.*mblade.*cos(theta2(t)).*diff(theta1(t),t,t)

Error in 413Design2>@(t,y)[y(2);EOM2_matlabFunction(t,y(1),y(2),theta1_val,thetadot1_val,L1,L2,mblade,kt,Jblade)] (line 77)

EOM2_matlabFunction(t, y(1), y(2), theta1_val, thetadot1_val, L1, L2, mblade, kt, Jblade)];

%EOM2 =  dt_dL_dthetadot2 - dL_dtheta2 == 0; %Lagrange second order diff equations
%odeFunction = @(t, y) [y(2); -(dL_dtheta2_numeric - dt_dL_dthetadot2_numeric) / Jblade]; % Previous First order converted from EOM

%Full code #2
clear; close all; clc;

% Step 1: Define symbolic variables
syms theta1(t) theta2(t) L1 L2 mblade kt Jblade theta2_initial

% Define time derivatives of theta variables
thetadot1 = diff(theta1, t);
thetadot2 = diff(theta2, t);

% Position of the second joint (assuming planar mechanism)
x2 = L1*cos(theta1) + L2*cos(theta1 + theta2);
y2 = L1*sin(theta1) + L2*sin(theta1 + theta2);

% Step 2: Kinetic and Potential Energy
% Velocity components of the second joint
xdot2 = diff(x2, t);
ydot2 = diff(y2, t);

% Kinetic Energy (Translation + Rotation)
KE_Translation = (1/2) * mblade * (xdot2^2 + ydot2^2);
KE_Rotation = (1/2) * Jblade * (thetadot1 + thetadot2)^2;

% Potential Energy (Spring term)
PETorSpring = (1/2) * kt * (theta2 - theta2_initial)^2;

% Lagrangian
T = KE_Translation + KE_Rotation;
V = PETorSpring;
L = T - V; % Lagrangian

% Derivatives of the Lagrangian for theta2 only
dL_dtheta2 = diff(L, theta2);
dL_dthetadot2 = diff(L, thetadot2);
dt_dL_dthetadot2 = diff(dL_dthetadot2, t);

% Step 3: Define the equation of motion expression
EOM2_expr = simplify(dt_dL_dthetadot2 - dL_dtheta2);

% Ensure that each input variable in 'Vars' is defined as a symbolic variable
vars = [sym('t'), sym('theta2_initial') sym('theta2'), sym('thetadot2'), sym('theta1'), sym('thetadot1'), sym('L1'), sym('L2'), sym('mblade'), sym('kt'), sym('Jblade')];

% Create the MATLAB function
EOM2_matlabFunction = matlabFunction(EOM2_expr, 'Vars', vars);

% Step 5: Define parameters (numeric values) and initial conditions
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)
Jblade = (1/12 * mblade * L2^2) + (mblade * (L1 + L2/2)^2); % Moment of inertia of the blade

% Convert RPM to rad/s for thetadot1 values (in each time interval)
time_intervals = [0, 30, 60, 90, 120];
thetadot1_rpm = [0, 50, 100, 150, 0];
thetadot1_rad = thetadot1_rpm * (2 * pi / 60);

% Initial conditions for theta2
theta2_initial_val = 30 * pi / 180;  % Initial angle in radians
thetadot2_initial_val = 0;           % Initial angular velocity

% Step 6: Solve the ODE using ode45 over multiple intervals
theta2_solution_values_total = [];
time_solution_values_total = [];

for i = 1:length(time_intervals) - 1
    % Update theta1 and thetadot1 in radians/s for the current interval
    theta1_val = 0;  % Fixed or known value
    thetadot1_val = thetadot1_rad(i);
    
    % Define the numeric ODE function for ode45
    odeFunction_i = @(t, y) [y(2); ...
                             EOM2_matlabFunction(t, y(1), y(2), theta1_val, thetadot1_val, L1, L2, mblade, kt, Jblade)];
    
    % Define the time span for the current interval
    tspan = [time_intervals(i), time_intervals(i + 1)];
    
    % Initial conditions for theta2 and thetadot2
    initial_conditions = [theta2_initial_val; thetadot2_initial_val];
    
    % Solve the ODE
    [t_interval, sol_interval] = ode45(odeFunction_i, tspan, initial_conditions);
    
    % Update initial conditions for next interval
    theta2_initial_val = sol_interval(end, 1);
    thetadot2_initial_val = sol_interval(end, 2);
    
    % Store solutions for plotting
    theta2_solution_values_total = [theta2_solution_values_total; sol_interval(:, 1)];
    time_solution_values_total = [time_solution_values_total; t_interval];
end

% Step 8: 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;

[Torsten]

Your expression EOM2_expr still contains second-order derivatives for theta1 and theta2. Thus you need 4 first-order differential equations for theta1, theta1', theta2 and theta2'.

[Miles]

Yes, you're right about the expression. That makes things tricky. so you are saying from that expression i need to convert it into 4 different first order diff eqs for theta 1, thetadot1, theta2, thetadot2? I tred odetovectorfield, but that didnt work.

However, I already have values for theta dot i want to input [0,50,100,150,0]. For theta, could i just assume radian values and cut it down to 2 first order diff equations or this would not work?

[Torsten]

That makes things tricky. so you are saying from that expression i need to convert it into 4 different first order diff eqs for theta 1, thetadot1, theta2, thetadot2?

I don't know the sense of this expression. Do you mean you have the equation "expression = 0" as part of your ODE system ?

However, I already have values for theta dot i want to input [0,50,100,150,0]. For theta, could i just assume radian values and cut it down to 2 first order diff equations or this would not work?

I don't know what you mean by "for theta dot i want to input [0,50,100,150,0]".

You must prescribe initial conditions at t = 0 for theta1, theta1dot, theta2 and theta2dot - thus four values. Further you need four first-order differential equations for theta1, theta1dot, theta2 and theta2dot to solve.

[Miles]

Im not sure honestly. That seems very complicated. I wanted to create a plot that shows the change of the angle theta 2 over time. I want to do this on a thetadot1 vs time graph, but i think in order to do that i need to find theta dot 2.

I made some changes to simplify the equations and try to use algebra to solve for the double derivative portion of theta2 in the second order diff equation, but even this isnt possible unless you assume the double derivative of theta 1 is equal to zero.

The lagrange equation for dynamics i learned was

d/dt(Dl/dqdotx)-dL/dqx = 0 for conservative systems with no applied forces

% Derivatives of the Lagrangian for theta1 and theta2
dL_dtheta1 = diff(L, theta1);
dL_dthetadot1 = simplify(diff(L, thetadot1));
dL_dtheta2 = simplify(diff(L, theta2));
dL_dthetadot2 = simplify(diff(L, thetadot2));

% Euler-Lagrange equations for theta1 and theta2
eq1 = simplify(diff(dL_dthetadot1, t) - dL_dtheta1); % Eq of motion for theta1
eq2 = simplify(diff(dL_dthetadot2, t) - dL_dtheta2);
% Eq of motion for theta2
d2theta1_dt2 = 0
d2theta2_dt2 = (- (kt * theta2(t) - kt * theta2_initial) - ...
                (Jblade + L2^2 * mblade + L1 * L2 * mblade * cos(theta2(t))) * d2theta1_dt2 - ...
                L1 * L2 * mblade * sin(theta2(t)) * (thetadot1)^2) / (Jblade + L2^2 * mblade);

d2theta2_dt2simple = simplify(d2theta2_dt2)

%Convert Equation 1 to first order diff equation

[V] = odeToVectorField(d2theta2_dt2)

%Convert Symbolic functions into a matlab function with all the variables
M = matlabFunction(V,'vars', {'t','Y'})

% Display the result
disp('The equation of motion for d^2 theta2/dt^2:');
disp(ddtheta2);

The equation you get when you run the code for

[Torsten]

The below code gives explicit equations for d^2theta1/dt^2 and d^2theta2/dt^2. I don't know how extra external forces are included.

Without the centrifugal force, your system of equation is

dtheta1/dt = theta1dot

dtheta1dot/dt = a

dtheta2/dt = theta2dot

dtheta2dot/dt = b

with a, b as below.

% 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)
syms a b
eq1 = subs(eq1,[diff(theta1,t,2),diff(theta2,t,2)],[a,b])
eq2 = subs(eq2,[diff(theta1,t,2),diff(theta2,t,2)],[a,b])
% a = d^2theta1/dt^2, b = d^2theta2/dt^2
sol = solve([eq1,eq2],[a b])

sol = struct with fields:
    a: (35*theta2(t) + 15*cos(theta1(t) + theta2(t))^2*theta2(t) + 15*sin(theta1(t) + theta2(t))^2*theta2(t) + 15*cos(theta1(t))*cos(theta1(t) + theta2(t))*theta2(t) + 15*sin(t...
    b: -(35*theta2(t) + 15*sin(theta1(t))^2*theta2(t) + 15*cos(theta1(t) + theta2(t))^2*theta2(t) + 15*sin(theta1(t) + theta2(t))^2*theta2(t) + 15*cos(theta1(t))^2*theta2(t) + ...

disp(char(sol.a))

(35*theta2(t) + 15*cos(theta1(t) + theta2(t))^2*theta2(t) + 15*sin(theta1(t) + theta2(t))^2*theta2(t) + 15*cos(theta1(t))*cos(theta1(t) + theta2(t))*theta2(t) + 15*sin(theta1(t))*sin(theta1(t) + theta2(t))*theta2(t) + 84*cos(theta1(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 - 84*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 84*cos(theta1(t))*sin(theta1(t) + theta2(t))*diff(theta2(t), t)^2 - 84*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta2(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t) + theta2(t))^3*diff(theta1(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))^3*diff(theta1(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t) + theta2(t))^3*diff(theta2(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))^3*diff(theta2(t), t)^2 - 36*cos(theta1(t))*sin(theta1(t))*cos(theta1(t) + theta2(t))^2*diff(theta1(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)^2 + 72*cos(theta1(t))*sin(theta1(t) + theta2(t))^3*diff(theta1(t), t)*diff(theta2(t), t) - 72*sin(theta1(t))*cos(theta1(t) + theta2(t))^3*diff(theta1(t), t)*diff(theta2(t), t) + 36*cos(theta1(t))*cos(theta1(t) + theta2(t))^2*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 36*cos(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 36*cos(theta1(t))*cos(theta1(t) + theta2(t))^2*sin(theta1(t) + theta2(t))*diff(theta2(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)^2 - 36*sin(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))^2*diff(theta2(t), t)^2 + 168*cos(theta1(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) - 168*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) + 72*cos(theta1(t))*cos(theta1(t) + theta2(t))^2*sin(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) - 72*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)*diff(theta2(t), t))/(12*(3*cos(theta1(t))^2*sin(theta1(t) + theta2(t))^2 + 3*sin(theta1(t))^2*cos(theta1(t) + theta2(t))^2 + 7*cos(theta1(t))^2 + 7*sin(theta1(t))^2 - 6*cos(theta1(t))*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))))

disp(char(sol.b))

-(35*theta2(t) + 15*sin(theta1(t))^2*theta2(t) + 15*cos(theta1(t) + theta2(t))^2*theta2(t) + 15*sin(theta1(t) + theta2(t))^2*theta2(t) + 15*cos(theta1(t))^2*theta2(t) + 30*cos(theta1(t))*cos(theta1(t) + theta2(t))*theta2(t) + 30*sin(theta1(t))*sin(theta1(t) + theta2(t))*theta2(t) + 84*cos(theta1(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 - 84*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 84*cos(theta1(t))*sin(theta1(t) + theta2(t))*diff(theta2(t), t)^2 - 84*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta2(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t) + theta2(t))^3*diff(theta1(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))^3*diff(theta1(t), t)^2 + 36*cos(theta1(t))^3*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 - 36*sin(theta1(t))^3*cos(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t) + theta2(t))^3*diff(theta2(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))^3*diff(theta2(t), t)^2 - 72*cos(theta1(t))*sin(theta1(t))*cos(theta1(t) + theta2(t))^2*diff(theta1(t), t)^2 - 36*cos(theta1(t))^2*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta1(t), t)^2 - 36*cos(theta1(t))*sin(theta1(t))*cos(theta1(t) + theta2(t))^2*diff(theta2(t), t)^2 + 72*cos(theta1(t))*sin(theta1(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t))^2*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 36*cos(theta1(t))*sin(theta1(t))*sin(theta1(t) + theta2(t))^2*diff(theta2(t), t)^2 + 72*cos(theta1(t))*sin(theta1(t) + theta2(t))^3*diff(theta1(t), t)*diff(theta2(t), t) - 72*sin(theta1(t))*cos(theta1(t) + theta2(t))^3*diff(theta1(t), t)*diff(theta2(t), t) + 36*cos(theta1(t))*cos(theta1(t) + theta2(t))^2*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 72*cos(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 + 36*cos(theta1(t))*cos(theta1(t) + theta2(t))^2*sin(theta1(t) + theta2(t))*diff(theta2(t), t)^2 + 36*cos(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta2(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)^2 - 72*sin(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)^2 - 36*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))^2*diff(theta2(t), t)^2 - 36*sin(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta2(t), t)^2 + 168*cos(theta1(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) - 168*sin(theta1(t))*cos(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) - 72*cos(theta1(t))*sin(theta1(t))*cos(theta1(t) + theta2(t))^2*diff(theta1(t), t)*diff(theta2(t), t) + 72*cos(theta1(t))*sin(theta1(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)*diff(theta2(t), t) + 72*cos(theta1(t))*cos(theta1(t) + theta2(t))^2*sin(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) + 72*cos(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t) - 72*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))^2*diff(theta1(t), t)*diff(theta2(t), t) - 72*sin(theta1(t))^2*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))*diff(theta1(t), t)*diff(theta2(t), t))/(12*(3*cos(theta1(t))^2*sin(theta1(t) + theta2(t))^2 + 3*sin(theta1(t))^2*cos(theta1(t) + theta2(t))^2 + 7*cos(theta1(t))^2 + 7*sin(theta1(t))^2 - 6*cos(theta1(t))*sin(theta1(t))*cos(theta1(t) + theta2(t))*sin(theta1(t) + theta2(t))))

[Miles]

I tested this out and it actually gives an output! Thank you so much for making this work! May I ask why solve works better than using ode45?

[Torsten]

"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,:))