How to Output k_mb for Plotting in ODE45?

  • Context: MATLAB 
  • Thread starter Thread starter Sirsh
  • Start date Start date
  • Tags Tags
    Ode45 Variables
Click For Summary
SUMMARY

The discussion focuses on how to output the variable "k_mb" from the ODE solver in MATLAB using the ode45 function for vibration analysis. The user provided a code snippet that defines the state space equations and global variables, including k_mb, which represents torsional stiffness. The main script initializes the global variables and sets up the ODE solver with specific tolerances. The solution is intended to visualize the behavior of k_mb during integration.

PREREQUISITES
  • Understanding of MATLAB programming, specifically with ODE solvers.
  • Familiarity with global variables in MATLAB.
  • Knowledge of state space representation in dynamic systems.
  • Basic concepts of vibration analysis and torsional stiffness.
NEXT STEPS
  • Learn how to visualize ODE solver outputs in MATLAB using plots.
  • Explore the use of global variables in MATLAB for dynamic simulations.
  • Investigate the impact of varying torsional stiffness on system behavior.
  • Study advanced options in the ode45 function for improved accuracy and performance.
USEFUL FOR

Engineers, researchers, and students involved in mechanical systems analysis, particularly those working with MATLAB for solving ordinary differential equations and analyzing vibration dynamics.

Sirsh
Messages
262
Reaction score
10
Hi all,

I have the following ode code, and I want to know how I could output the variable "k_mb" within my main script so that I can plot it to see if it's actually active when the solver is integrating.
odescript
Matlab:
function zdot = odenesttest(t,z)
global r_p r_g I_p I_g m_p m_g V k_mb
A_or_B = [5,10];
V = [0,5,15,20,30,35,45,50,60];
%% state space equations are defined below:
zdot(2) = ((r_p*k_mb)/I_p)*(-r_p*z(1)+r_g*z(3)+z(5)-z(7));  % Eq 1ss
zdot(6) = (k_mb/m_p)*(r_p*z(1)-r_g*z(3)-z(5)+z(7));         % Eq 2ss
zdot(4) = ((r_g*k_mb)/I_g)*(-r_g*z(3)+r_p*z(1)+z(7)-z(5));  % Eq 3ss
zdot(8) = (k_mb/m_g)*(-r_p*z(1)+r_g*z(3)-z(7)+z(5));        % Eq 4ss

zdot(1) = z(2); % Eq 5ss
zdot(3) = z(4); % Eq 6ss
zdot(5) = z(6); % Eq 7ss
zdot(7) = z(8); % Eq 8ss

zdot(9) = z(1); % Pinion angular displacement (control angle).
torsional_stiffness_changer
   function torsional_stiffness_changer
       if 0 <= zdot(9)
           lower_index = find(V < zdot(9), 1, 'last');
           even_or_odd = mod(lower_index,2);
           k_mb = A_or_B(2 - even_or_odd);
       end
   end
% Outputs
zdot = zdot';
end

Mainscript
Matlab:
% Mainscript for vibration analysis (19/01/2017).
clear all;
clc;
%% Defining global variables
global r_p r_g I_p I_g m_p m_g k_mb
k_mb = 50;  % torsional stiffness.
r_p = 10;   % pinion pitch(?) radius
r_g = 15;   % gear pitch(?) radius
I_p = 100;  % pinion inertia.
I_g = 200;  % gear inertia.
m_p = 50;   % pinion mass
m_g = 100;  % gear mass
%==========================================================================
%% initial coditions for ode solver
%==========================================================================
z0 = [1,0,0,0,0,0,0,0,0]; % zero initial conditions.
t_step = [0:1:10];        % timestep 0 to 5 seconds.
opts = odeset('Reltol', 1e-3, 'AbsTol', 1e-3);  % conservative tolerances.
[t,z] = ode45(@odenesttest, t_step, z0, opts);
%==========================================================================
%% Screen Output and Plotting
%==========================================================================
[t,z(:,:)]

Thanks a lot!
 
Physics news on Phys.org
if it is a global, you should be able to see it all the time, write it to a file or vector and update the plot...
 
Dr Transport said:
if it is a global, you should be able to see it all the time, write it to a file or vector and update the plot...

Sorry, in the mainscript code that I listed I meant to take out the "k_mb = 50;". I have edited it and put it below.

Updated mainscript
Code:
% Mainscript for vibration analysis (19/01/2017).
clear all;
clc;
%% Defining global variables
global r_p r_g I_p I_g m_p m_g k_mb
r_p = 10;   % pinion pitch(?) radius
r_g = 15;   % gear pitch(?) radius
I_p = 100;  % pinion inertia.
I_g = 200;  % gear inertia.
m_p = 50;   % pinion mass
m_g = 100;  % gear mass
%==========================================================================
%% initial coditions for ode solver
%==========================================================================
z0 = [1,0,0,0,0,0,0,0,0]; % zero initial conditions.
t_step = [0:1:10];        % timestep 0 to 5 seconds.
opts = odeset('Reltol', 1e-3, 'AbsTol', 1e-3);  % conservative tolerances.
[t,z] = ode45(@odenesttest, t_step, z0, opts);
%==========================================================================
%% Screen Output and Plotting
%==========================================================================
[t,z(:,:)]
 
Last edited:

Similar threads

  • · Replies 6 ·
Replies
6
Views
4K
  • · Replies 2 ·
Replies
2
Views
4K
  • · Replies 4 ·
Replies
4
Views
2K
  • · Replies 5 ·
Replies
5
Views
5K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 5 ·
Replies
5
Views
9K
  • · Replies 2 ·
Replies
2
Views
7K
  • · Replies 34 ·
2
Replies
34
Views
4K
  • · Replies 5 ·
Replies
5
Views
6K
Replies
1
Views
2K