MATLAB How to Output k_mb for Plotting in ODE45?

  • Thread starter Thread starter Sirsh
  • Start date Start date
  • Tags Tags
    Ode45 Variables
AI Thread Summary
The discussion focuses on how to output the variable "k_mb" from an ODE45 solver for plotting purposes. The user is trying to determine if "k_mb" is active during the integration process. It is clarified that since "k_mb" is defined as a global variable, it should be accessible throughout the script. Suggestions include writing "k_mb" to a file or vector for real-time updates in the plot. The user also corrected the main script by removing an initial assignment to "k_mb" to ensure it reflects the dynamics defined in the ODE function.
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

Back
Top