Efficient MATLAB Repeating Function for User-Defined Rotations

[Sirsh]

Hey all,

I need your help writing the meat of a function that I need for a MATLAB script I'm writing.

The jist of it is, starting from 0°, when between 0.0° - 5.0° apply a certain value, when between 5.0° - 15.0° apply another value. However, the relative amount of rotation can be user defined, so I cannot make predefined matrices. The upper and lower limits are absolute, but the intervals within them should be able to be any size.

I thought it would be possible through a if else statement, but I cannot figure out how to make the condition of the statement applicable to what I've stated above. As if the user input rotation is 360°, value 'A' would be applied at 0°-5°, 15°-20° 30°-35° and value 'B' would be applied at 5°-15°, 20°-30°, 35°-45° etc.

Any help would be appreciated, it's probably cake for you seasoned programmers!

 

[Mark44]

Hey all,

I need your help writing the meat of a function that I need for a MATLAB script I'm writing.

The jist of it is, starting from 0°, when between 0.0° - 5.0° apply a certain value, when between 5.0° - 15.0° apply another value. However, the relative amount of rotation can be user defined, so I cannot make predefined matrices. The upper and lower limits are absolute, but the intervals within them should be able to be any size.

I thought it would be possible through a if else statement, but I cannot figure out how to make the condition of the statement applicable to what I've stated above. As if the user input rotation is 360°, value 'A' would be applied at 0°-5°, 15°-20° 30°-35° and value 'B' would be applied at 5°-15°, 20°-30°, 35°-45° etc.

Any help would be appreciated, it's probably cake for you seasoned programmers!

I think what would work is a cascade of if...else statements. This or something close to it would work, I think. I don't have matlab, so you will need to verify that my syntax is correct.

if (0.0 < angle && angle < 5.0)
   % do something
elseif (15.0 < angle <  20.0)
   % do something else
.
.
.
etc.

 

[Sirsh]

I think what would work is a cascade of if...else statements. This or something close to it would work, I think. I don't have matlab, so you will need to verify that my syntax is correct.

if (0.0 < angle && angle < 5.0)
   % do something
elseif (15.0 < angle <  20.0)
   % do something else
.
.
.
etc.

Do you mean for each interval there would be an if statement associated with it? This is the nature of the application that I require for an example of 60°:
https://postimg.org/image/xcz38o687/
See how the function needs to be able to periodically apply the values dependent on the total angle defined?

 

[Mark44]

Do you mean for each interval there would be an if statement associated with it? This is the nature of the application that I require for an example of 60°:
https://postimg.org/image/xcz38o687/

From your image, 60° seems to be the total angle, which is broken up into subintervals. What I showed above would work for this.

Does you code also need to work for different situations (like a different range of angles, with possibly different subintervals)?

See how the function needs to be able to periodically apply the values dependent on the total angle defined?

It isn't clear to me what you mean here.

 

[FactChecker]

If you want the values to alternate as your example shows, you can use the 'find' and 'mod' functions.

Suppose your breakpoint values are in vector
V = [0, 5, 15, 20, 30, 35, 45, 50, 60];
And your A/B values are in the vector
A_or_B = [5, 10];
Then
lower_index = find( V < angle, 1); % will return the first index that is less than angle.
even_or_odd = mod( lower_index, 2); % will return a 0 if lower_index is even and 1 if it is odd.
a_or_b_for_angle = A_or_B(even_or_odd); % will have the appropriate value of A or B

You may have to work with this some. I can't test it and I am sure that there are some corrections needed.

 

[Sirsh]

It isn't clear to me what you mean here.

Maybe it'll be more clear if I explain the application. I am trying to make a function that allows me to apply a torsional stiffness for gear teeth in mesh, as the gear rotates the gear (in this case) will be in either singular contact or double contact. Where I'm assigning 5° for single and 10° for double, so if I want to simulate 10 rotations (3600°) then I need the function to be able to apply the torsional stiffness in 5° and 10° intervals.

I'll be running and ODE solver which will be coupled with this function I need to make so it can evaluate the torsional stiffness at any point and be able to it apply to the dynamic equations representing the physical system as they're solved.

You may have to work with this some. I can't test it and I am sure that there are some corrections needed.

I've chucked the code you wrote into my MATLAB editor, and associated an arbitrary value for the 'angle', set at 60°. The outputs don't make much sense to me, can you explain in laymen terms what the code is doing so I can have a play with it some more? Thanks.

 

[FactChecker]

I've chucked the code you wrote into my MATLAB editor, and associated an arbitrary value for the 'angle', set at 60°. The outputs don't make much sense to me, can you explain in laymen terms what the code is doing so I can have a play with it some more?

I tried to do that in the comments in the right part of each line, but maybe they were not clear. Here are some added comments. I tried to correct a mistake in the last line but I still can't test it.

% Initialize vectors V with the values where the A/B switch should occur:
V = [0, 5, 15, 20, 30, 35, 45, 50, 60];
% Initialize vector A_or_B with your values of A and B
A_or_B = [5, 10];

% Find the index of the value in V that is just below your angle
lower_index = find( V < angle, 1); % will return the first index that is less than angle.

% Determine if the index is even or odd.
even_or_odd = mod( lower_index, 2); % will return a 0 if lower_index is even and 1 if it is odd.

% Odd indices begin the A value section and even indices begin the B value section.
a_or_b_for_angle = A_or_B( 2 - even_or_odd); % will have the appropriate value of A or B

 

[Sirsh]

I tried to do that in the comments in the right part of each line, but maybe they were not clear. Here are some added comments. I tried to correct a mistake in the last line but I still can't test it.

Again, many thanks. I'm quite sure I know what it is trying to achieve.

If I run:

angle = 4;
V = [0,5,15,20,30,35,45,50,60];
A_or_B = [5,10];

lower_index = find(V < angle, 1);
even_or_odd = mod(lower_index,2);
a_or_b_for_angle = A_or_B(2 - even_or_odd);

The values of each vector are:
a_or_b_for_angle = [5]
angle = [4]
even_or_odd = [1]
lower_index = [1]

However, if I run the same code except:
angle = 12;

The values for each vector are:
a_or_b_for_angle = [5]
angle = [12]
even_or_odd = [1]
lower_index = [1]

I am assuming that for the case of angle = 12, that the a_or_b_for_angle should be equal to variable B, which is 10?

 

[FactChecker]

Again, many thanks. I'm quite sure I know what it is trying to achieve.

If I run:

angle = 4;
V = [0,5,15,20,30,35,45,50,60];
A_or_B = [5,10];

lower_index = find(V < angle, 1);
even_or_odd = mod(lower_index,2);
a_or_b_for_angle = A_or_B(2 - even_or_odd);

The values of each vector are:
a_or_b_for_angle = [5]
angle = [4]
even_or_odd = [1]
lower_index = [1]

However, if I run the same code except:
angle = 12;

The values for each vector are:
a_or_b_for_angle = [5]
angle = [12]
even_or_odd = [1]
lower_index = [1]

I am assuming that for the case of angle = 12, that the a_or_b_for_angle should be equal to variable B, which is 10?

Yes, that is the desired result of this code. I don't know if that is what you want. Obviously there is a bug. If you want to debug this code, look at why lower_index is 1 in the second case. It should have been 2.

PS. This should work better. It will find the largest index whose V value is less than angle:
lower_index = find(V < angle, 1, 'last');

 

[Sirsh]

Yes, that is the desired result of this code. I don't know if that is what you want. Obviously there is a bug. If you want to debug this code, look at why lower_index is 1 in the second case. It should have been 2.

PS. This should work better. It will find the largest index whose V value is less than angle:
lower_index = find(V < angle, 1, 'last');

The addition of 'last' in the lower_index, fixed the problem instantly. Thus, for:
angle = 4 -> a_or_b_for_angle = 5
angle = 6 -> a_or_b_for_angle = 10
angle = 46 -> a_or_b_for_angle = 5
angle = 59 -> a_or_b_for_angle = 10

Thank you very much for your help, I'll use this as my base for my ODE solvers to use.

 

[FactChecker]

The addition of 'last' in the lower_index, fixed the problem instantly. Thus, for:
angle = 4 -> a_or_b_for_angle = 5
angle = 6 -> a_or_b_for_angle = 10
angle = 46 -> a_or_b_for_angle = 5
angle = 59 -> a_or_b_for_angle = 10

Thank you very much for your help, I'll use this as my base for my ODE solvers to use.

If speed is an issue, you might want to modify it to find the first V value greater than angle; then subtract 1. The 'last' forces it to look through the entire V array (unless it is smart enough to go backward; that would be fast with no changes).

 

[Sirsh]

If speed is an issue, you might want to modify it to find the first V value greater than angle; then subtract 1. The 'last' forces it to look through the entire V array (unless it is smart enough to go backward; that would be fast with no changes).

Once I start testing the ODE solver, I'll be sure to keep that in mind if it becomes an issue.

One further thing, see how the array, V = [0,5,15,20,30,35,45,50,60]; (pre-determined and put in). I want to write something that allows the user to input the angle and then the array is generated. So far this is the beginning of what i'll use:

t_rotation = input('Enter the integer rotation: ');
if isempty(t_rotation)
t_rotation = 1;
end
t_rot_degree = t_rotation*360; % converts the integer rotation to degrees.

I've been toying with some sort of series summation commands that MATLAB has, but I'm not too sure. I know that V must start at 0, always, and from here the series is 5,10,5,10,...,etc. If the user specified 60 degrees, then there will be four 5's and four 10's in the series, I am just not sure how to make something add them in series to create the array V like you've shown.

 

[FactChecker]

It's not clear to me what is fixed and what varies. Does 60° change? Does the 5,10 alternating pattern change? It's not clear what angle the user is inputting and how it effects those numbers.

 

[Sirsh]

It's not clear to me what is fixed and what varies. Does 60° change? Does the 5,10 alternating pattern change? It's not clear what angle the user is inputting and how it effects those numbers.

Sorry for the confusion, let me explain. The purpose of this is to simulate the revolution of a pair of gears in mesh, so if the user wants to simulate 10 revolutions the angle will be 3600°, if they want to rotate the gears 1 revolution then the angle will be 360°, etc. So the angle is the only variable that is user dependent. However, the 5,10 alternating pattern doesn't change at all.

So for example, if 0.5 revolutions is desired, then the array that I need to somehow produce should be:
V = [0,5,15,20,30,35,45,50,60,65,75,80,90,95,105,110,120,125,135,140,150,155,165,170,180]

As doing this manually for 10+ revolutions will be extremely tedious.

I thought that it would be possible to make an array of 5's, 10's, and an array with the first index equal to zero, then use some sort of summation command to produce the above array, but I'm not too sure at the moment.

Thanks again for your help!

 

[FactChecker]

Thanks. That is much clearer. Try something like this. (Again, I can not test the code, so you will probably have to work on it some.)

% set number of revolutions to simulate and the corresponding maximum angle
num_revolutions = 10.5;
max_angle = num_revolutions * 360;

% initialize the angle and index to start at the first V entry
current_angle = 0;
current_index = 1;
V(1) = 0;

% loop through angles till max_angle is exceeded
while current_angle le max_angle
   if mod( current_index, 2) eq 1
      % for odd indices the next angle should be 5 larger
      V(current_index+1) = V(current_index) + 5;
   else
      % for even indices the next angle should be 10 larger
      V(current_index+1) = V(current_index) + 10;
   end
   % move the index and angle up for the next loop
   current_index = current_index + 1;
   current_angle = V(current_index);
end

PS. The Math Software and LaTeX forum section under Mathematics is a better place to ask MATLAB questions.

 

[Sirsh]

Thanks. That is much clearer. Try something like this. (Again, I can not test the code, so you will probably have to work on it some.)

PS. The Math Software and LaTeX forum section under Mathematics is a better place to ask MATLAB questions.

This is great, works perfectly with the addition of this at the end:

if V(current_index) ~= max_angle
    V(current_index) = [];
else
    V(current_index);
end

Otherwise the last entry in the V array is either over by 5 or 10. But with that, it works nicely.

I'll be sure to check out those sub-forums, it never crossed my mind to post this in there actually!

 

[Sirsh]

PS. The Math Software and LaTeX forum section under Mathematics is a better place to ask MATLAB questions.

Just one more thing (I didn't think it'd be worthwhile beginning a new thread for this when you've responded in here). I'm now trying to incorporate the repeating function into my ODE solving function which obviously requires a bit of modification to the code we have in the previous posts. I'm currently using ODE45 in the form of:

function zdot = odess(t,z)

where my 'z' is a state space variable. Now, the problem is, that I need the stiffness varied as the angle changes as we have achieved previously, however if my initial conditions for the ode solver are zero and the solver begins iterating the code, the stiffness switching needs to be able to know what angle the ode solver is at and which variable value to apply. As the 'max_angle' is not defined by the user in this case, the angle is whatever the solver gets to when the time step has finished (I think).

If the state space variable z(2) is my angle, then I think 'current_angle' can be changed to this, but I'm not too sure how to dictate or indicate the max_angle based on the solvers dynamics.

Below is my current ODE solver function code, implanted with the code that is in the previous post (obvious will be changed). Any thoughts to this?

function zdot = odess(t,z)
global k_mb r_p r_g I_p I_g m_p m_g % defining global variables

% initialize the angle and index to start at the first V entry
current_angle = 0;
current_index = 1;
V(1) = 0;

% loop through angles till max_angle is exceeded
while current_angle <= max_angle
   if mod(current_index, 2) == 1
      % for odd indices the next angle should be 5 larger (single contact)
      V(current_index+1) = V(current_index) + 5;    % 5 can be user specified in later versions.
   else
      % for even indices the next angle should be 10 larger (double
      % contact)
      V(current_index+1) = V(current_index) + 10;   % 10 can be user specified in later versions.
   end
   % move the index and angle up for the next loop
   current_index = current_index + 1;
   current_angle = V(current_index);
end

if V(current_index) ~= max_angle % Truncates the last entry if not equal.
    V(current_index) = [];
else
    V(current_index);   % If equal, does nothing and ends.
end

%% Adding in the next step - extracted from 'sfunc.m'. 
A_or_B = [5,10];    % array holding the variable values (i.e. stiffnesses).

lower_index = find(V < current_angle, 1, 'last');   % array V, generated in prev section.
even_or_odd = mod(lower_index,2);
a_or_b_for_angle = A_or_B(2 - even_or_odd);

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 = zdot';

 

[Sirsh]

One easy way to solve this I think is to make a preset array like 'V' extremely long and then change the current angle to z(2).

 

[FactChecker]

One easy way to solve this I think is to make a preset array like 'V' extremely long and then change the current angle to z(2).

I'm not sure that I understand what you want. If you need to pass the value of max_angle in, you can add it to the input parameter list. If you need to pass the value out, you can return it from the function: function [zdot, max_angle] = odess( ...

 

[Sirsh]

I'm not sure that I understand what you want. If you need to pass the value of max_angle in, you can add it to the input parameter list. If you need to pass the value out, you can return it from the function: function [zdot, max_angle] = odess( ...

Sorry, I'm still an intermediate at using the ode solvers with matlab. I'll list below what I have written (excluding what you have shown me previously).

Main script:

% Mainscript for vibration analysis (16/01/2017).
clear all;
clc;
%% Defining global variables
global k_mb r_p r_g I_p I_g m_p m_g
k_mb = 10;  % 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 (using ode45 initially)
%==========================================================================
z0 = [0,0,0,0,0,0,0,0]; % zero initial conditions.
t_step = [0 10];        % timestep 0 to 10 seconds.
opts = odeset('Reltol', 1e-3, 'AbsTol', 1e-3);  % conservative tolerances.
[t,z] = ode45(@odess, t_step, z0, opts);

ODE function:

% This script is for the 'guts' of the ode function, where state space
% equations live.

function zdot = odess(t,z)
% defining global variables
global k_mb r_p r_g I_p I_g m_p m_g

%% 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 = zdot';

As you can see in the main script, k_mb is the stiffness (what I've been trying to create a repeating function for). However, I've just set it to a fixed value for the time being. As I previously stated the stiffness needs to change as the gears rotate between 5 and 10 etc. However, that is based on a maximum angle which is all good and well, except in the ODE solver script the maximum angle isn't defined, the timestep is, so I'm trying to figure out how I can use the angular displacement (z(2)) to make the stiffness applied. If that makes sense? Cheers.