ODE45 for Equations of Motion

Hi, I am hoping someone can help me with the following issue. I have tried this problem multiple times and still can't seem to write the correct code (I am new to MATLAB).

The following three equations of motion are to be solved:

theta dot out is the input velocity into the system, it is time dependant and is characterised by the equation:

All other variables are constant (e.g. theta_dot_zero,J1,c1,k1,f (this is not frequency) etc.).

I believe you can use ode45 to solve this however I have had no luck.

Initial conditions can be taken as zero, a timespan of 0-20seconds can be used as an example.

The aim is to then plot velocity, acceleration and displacement graphs for the unknown theta variables of subscript 1,2,n and out.

I hope someone can help with this! I would really appreciate a draft of some code.

Thanks!

[Steven Lord]

This sounds like a homework assignment. If it is, show us the code you've written to try to solve the problem and ask a specific question about where you're having difficulty and we may be able to provide some guidance.

If you aren't sure where to start because you're not familiar with how to write MATLAB code, I suggest you start with the free MATLAB Onramp tutorial to quickly learn the essentials of MATLAB.

If you aren't sure where to start because you're not familiar with the mathematics you'll need to solve the problem, I recommend asking your professor and/or teaching assistant for help.

[Sam Butler]

Hi Steven, this is not a homework problem. It is to produce theoretical models based to compare to experimental data in the future.

I have started by just using the 1st equation with only theta1 and thetaout as the unknows.

My function file:

function dx = ujRig(t,Y,param)

theta_out = Y(1); 
theta_out_dot = Vin+(0.02*vin*cos(72*pi*f_mech*param.t))+(0.005*vin*cos(144*pi*f_mech*param.t));
theta_one = Y(2); theta_one_dot = Y(3);
theta_one_dotdot = -1/param.J1((param.c1*(theta_one_dot-theta_out_dot))+(param.k1(theta_one-theta_out)));

dx = [theta_out;theta_one_dot;theta_one_dotdot];

end

My script file:

param.vin = 600*(2*pi/60); %% 600 is the input velocity in rpm
param.f_mech = param.vin/(2*pi);
param.J1 = 2.08e-4;
param.c1 = 0.15;
param.k1 = 270;


% initial condition
Y0 = [0;0;0];

tMax =50; tMin = 0; % time

% solving the ODE options
options = odeset('reltol',1e-6,'abstol',1e-6);

%% time data
numFrequenciesUsed = 2e-4; %% Time step (data points)
timelist = linspace(tMin,tMax,numFrequenciesUsed)'; %% define time span

% initial condition
Y0 = [0;0;0];

for i = 1:length(timelist) %%from time at the start to time at the end
    
    i
    
    param.t = timelist(i);
    
    %% omega list
    numCyclesOfSimulation = 50;
    max_input_velocity =1200; %% 
    numPointsPerCycle = 50;
    input_span = linspace(0,max_input_velocity,numPointsPerCycle*numCyclesOfSimulation)';

    [tlist,Ylist] = ode45(@ujRig,input_span,Y0,options,param);
    
    % resetting initial condition
    Y0 = Ylist(end,:)';
    

    
    
end

The idea is to plot graphs for 600rpm-1200rpm based on what i input. For now, just starting with 600rpm. Have i defined the variables and the equations correctly in my function file? does the for loop give the program enough information to then plot the graphs of theta_out_dot and theta_one_dot?

Thank you for your help.

[Torsten]

@Sam Butler

Vin and vin are not interpreted as the same variable.

And to use ode45 to solve system (6), you have to rewrite it as a system of first-order differential equations. That means you will have 6 variables to solve for, not 3.

[Sam Butler]

Ah thank you, I did not spot this.

Sorry I don't understand what you mean by that, I have split them up into their first order equations. Would you mind altering the code to explain further please?

[Torsten]

I mean that Y0 and dx must be 6x1 vectors, not 3x1 vectors. They should be

Y0 = [theta_1,theta_1_dot,theta_2,theta_2_dot,theta_n,theta_n_dot] and

dx = [theta_1_dot,theta_1_dotdot,theta_2_dot,theta_2_dotdot,theta_n_dot,theta_n_dotdot].

You don't need to solve an additional differential equation for theta_out. If theta_0_dot is a constant and you know theta_0 and theta_out at t=0, it should be no problem for you to integrate the equation theta_out_dot = ... analytically to get theta_out.

[Sam Butler]

Hi Torsten, thank you for the help so far, I appreciate your support.

For now, to keep it simple, I am just using the input velocity (theta_out_dot) equation and theta_one_dotdot and just not use the other terms in the equation (e.g. 2 and n), as I just want to make sure I can do this on a simpler scale first.

I have made a few changes :

function dx = ujRig(t,Y,param)

theta_out = Y(1); 
theta_out_dot = param.vin+(0.02*param.vin*cos(72*pi*param.f_mech*t))+(0.005*param.vin*cos(144*pi*param.f_mech*t));
theta_one = Y(2); theta_one_dot = Y(3);
theta_one_dotdot = -1/param.J1*((param.c1*(theta_one_dot-theta_out_dot))+(param.k1*(theta_one-theta_out)));

dx = [theta_out_dot;theta_one_dot;theta_one_dotdot];


end

param.vin = 600*(2*pi/60); %% 600 is the input velocity in rpm
param.f_mech = param.vin/(2*pi);
param.J1 = 2.08e-4;
param.c1 = 0.15;
param.k1 = 270;


% initial condition
Y0 = [0;0;0];

tMax =5; tMin = 0; % time

% solving the ODE options
options = odeset('reltol',1e-6,'abstol',1e-6);

dt = 2e-4;
fs=1/dt;

%% time data
numFrequenciesUsed = fs; %% Time step (data points)
timelist = linspace(tMin,tMax,fs)'; %% define time span
OscAmpl_list_fwd = zeros(size(timelist));

% initial condition
Y0 = [0;0;0];

for i = 1:length(timelist) %%from time at the start to time at the end
    
    i
    
    param.t = timelist(i);
    
    %% omega list
    numCyclesOfSimulation = 50;
    max_input_velocity =1200; %% 
    numPointsPerCycle = 50;
    input_span = linspace(0,max_input_velocity,numPointsPerCycle*numCyclesOfSimulation)';

    [tlist,Ylist] = ode45(@ujRig,input_span,Y0,options,param);
      OscAmpl_list_fwd(i) = max(Ylist(end-5*numPointsPerCycle:end,1)); % steady state amplitude
      OscAmpl_list_fwd_vel(i) = max(Ylist(end-5*numPointsPerCycle:end,2)); % steady state amplitude
    
 OscAmpl_list_fwd(i) = OscAmpl_list_fwd(i);

    % resetting initial condition
    Y0 = Ylist(end,:)';
    


    
end
    figure
    plot(timelist,OscAmpl_list_fwd_vel)

I get no errors but it seems to take a really long time to compute (to the point where I cancel the simulation). So I dont know if i am stuck in an infinite loop or if it is to do with my line:

    timelist = linspace(tMin,tMax,fs)'; %% define time span

Hopefully this helps with the explanation.

[Torsten]

Do you see the problem why the integration takes so long and is almost impossible to perform ?

param.vin = 600*(2*pi/60); %% 600 is the input velocity in rpm
param.f_mech = param.vin/(2*pi);
theta_out_dot = @(t)param.vin+(0.02*param.vin*cos(72*pi*param.f_mech*t))+(0.005*param.vin*cos(144*pi*param.f_mech*t));
t=0:0.000001:0.1;
plot(t,theta_out_dot(t))

[Shree Charan]

@Sam Butler In the question you mention that "theta_dot_out" is the input velocity into the system and "theta_dot_zero" is a constant.

However in the line

theta_out_dot = param.vin+(0.02*param.vin*cos(72*pi*param.f_mech*t))+(0.005*param.vin*cos(144*pi*param.f_mech*t));

theta_dot_zero is replaced with vin which is the input velocity.

Can you please clarify if theta_out_dot = vin as mentioned in the question or if theta_dot_zero = vin.