ODE45 for Equations of Motion

Question

Sam Butler 2023-3-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932750-ode45-for-equations-of-motion

评论： Shree Charan 2023-5-4

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!

8 个评论
显示 6更早的评论隐藏 6更早的评论

Sam Butler 2023-3-21

编辑：Sam Butler 2023-3-21

在 MATLAB Online 中打开

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.

Sam Butler 2023-3-22

在 MATLAB Online 中打开

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 2023-3-22

编辑：Torsten 2023-3-22

在 MATLAB Online 中打开

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 2023-5-4

在 MATLAB Online 中打开

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

请先登录，再进行评论。

请先登录，再回答此问题。

ODE45 for Equations of Motion

8 个评论
显示 6更早的评论隐藏 6更早的评论

回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

ODE45 for Equations of Motion

8 个评论 显示 6更早的评论隐藏 6更早的评论

回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论