Help with ODE tolerance options

I am modelling a single stage helical gear transmission, and want to adjust the ODE tolerance options to allow for faster simulation times. Is there a way to figure out the tolerance values accurate for my model. I have tried various values for REltol and Abstol, but cannot get accurate results with that.

My main code-

function [yp] = reduced2(t,y,T_a)  
beta = 13*(pi/180);                      %Helix angle (rad)
P = 20*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 

Freq = 1000*20/60;

R_a = 51.19e-3;                             %Radius
R_p = 139.763e-3;                           %Radius

% Common parameters
e_a = (pi*2*R_a*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
e_p = (pi*2*R_p*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
% e = 0;
e = (pi*2*(R_a+R_p)*tan(beta))/(4*cos(P));

z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
Cs = 0.1 + 0.0001*sin(2*pi*Freq*t);                     %Shear damping


% Driver gear  
m_a = 0.5;                      %mass of driver gear (kg)
c_ay=500;                       %Damping of driver gear in y direction (Ns/m)
c_az=500;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 8e7;                      %Stiffness of driver gear in y direction (N/m)
k_az= 5e7;                      %Stiffness of driver gear in z direction (N/m)
I_a = 0.5*m_a*(R_a^2);          %Moment of inertia of driver gear (Kg.m3)

% y_ac= e_a + theta_a_vec*R_a;    %circumferential displacement at the meshing point on the driver gear (m)
% y_pc = e_p - theta_a_vec*R_a;   %circumferential displacement at the meshing point on the driven gear (m)

z_a = e_a*tan(beta);  
z_p = e_p*tan(beta);  
% z_ac = (z_a-y_ac)*tan(beta);    %axial displacements of the meshing point on the driving gear  (m)
% z_pc = (z_p+y_pc)*tan(beta);    %axial displacements of the meshing point on the driven gear   (m)

% Driven gear  
m_p = 0.5;                                  %mass of driven gear (kg)
c_py=500;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=500;                                   %Damping of driven gear in z direction (Ns/m)
k_py= 9.5e7;                                %Stiffness of driven gear in y direction (N/m)
k_pz =9e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 0.5*m_p*(R_p^2);                      %Moment of inertia of driven gear (Kg.m3)
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio

% % Excitation forces 
% Fy=ks*cos(beta)*(y_ac-y_pc-e) + Cs*cos(beta)*2*R_a*speed*2*pi;                %axial dynamic excitation force at the meshing point (N)
% Fz=ks*sin(beta)*(z_ac-z_pc-z) - Cs*sin(beta)*2*tan(beta)*R_a*speed*2*pi;      %circumferential dynamic excitation force at the meshing point (N)


%Time interpolation  
time = 0:0.00001*2:10;
% Torque = interp1(time,T_a,t)/1000;

% torque_range = T_a;  
% torque_interp = griddedInterpolant(time, torque_range);
% Torque = torque_interp(t) / 1000; 

Torque = griddedInterpolant(time, T_a);
Torque = Torque(t)/1000;


% torque_table = [time', T_a];
% torque = interp1(torque_table(:,1), torque_table(:,2), t, 'linear', 'extrap') / 1000;
Opp_Torque = 68.853/1000;                        % Average torque value

%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
ks = 0.9e8 + 20000*sin(2*pi*Freq*t);                   %Shear stiffness 
TE = y(11)*R_p - y(5)*R_a;                   %transmission error
b = 20e-6;                                   %Backlash

if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end

   
yp(1) = y(2); 
yp(2) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Cs*cos(beta)*R_a*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_a*(R_a*y(5)+R_p*y(11)) -KS*cos(beta)*R_a*(e_a-e_p-e))/I_a; %angular acceleration


%Driven gear equations
yp(7) = y(8);  
yp(8) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (Opp_Torque*i - Cs*cos(beta)*R_p*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_p*(R_a*y(5)+ R_p*y(11)) -KS*cos(beta)*R_p*(e_a-e_p-e))/I_p;                % angular accelration (rad/s2) 
end

My command window code-

tic
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:2:1000001);                  %Torque (Nm) chnages with time step

speed = 1000/60; 
tspan=0:0.00001*2:10;    %can try 4 or 5

y0 = [0;0;0;0;0;-104.719;0;0;0;0;0;104.719/3]; 
opts = odeset('RelTol',1e-1,'AbsTol',1e-1);
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);  
toc

% Driver gear graphs
nexttile  
plot(t,y(:,2))
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in y direction vs time')
axis padded  
grid on

nexttile  
plot(t,y(:,4))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in z direction vs time')  
axis padded  
grid on

nexttile  
plot(t,y(:,6))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driver gear angular velocity vs time')  
axis padded  
grid on

% Driven gear graphs
nexttile  
plot(t,y(:,8))  
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in y direction vs time')  
axis padded  
grid on

nexttile  
plot(t,y(:,10))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in z direction vs time')  
axis padded  
grid on

nexttile  
plot(t,y(:,12))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driven gear angular velocity vs time')  
axis padded  
grid on

% nexttile  
% plot(t,T_a(1:60001))  
% ylabel('Torque (Nm)')
% xlabel('Time (sec)')
% title('Torque vs time')  
% axis padded  
% grid on

My torque file- torque_for_Sid_60s - Copy.mat

[Jan]

You are still using a linear interpolation and hard steps for h and KS. This will attack the stepsize control and can reduce the accuracy of the calculations massively. I've told you the problem before already. It would not be efficient, to repeat this advice.

[Siddharth Jain]

Yes, I have tried spline. But it takes longer. If I can just ignore accuracy for the moment and want to get results faster, how can I Taylor my ode tolerance for that?

[Jan]

If run time matters, move the creation of the griddedInterpolant out of the loop, as I've suggested in your other question.

A fast dirty solution is not useful, because the applied method of integrating a non-smooth function causes numerical troubles. Starting with the fine-tuning by choosing optimal tolerances is not the point: The approach is not accurate.

[Mike Croucher]

I didn't notice you had already suggested that @Jan - sorry. That's pretty much all my answer below is about.

Answer 1

[Mike Croucher]

Thanks so much for including all of your code and data. It allowed me to run everything on my machine and use the profiler to figure out what was going on. In your function reduced2, you do this

time = 0:0.00001*2:10;
Torque = griddedInterpolant(time, T_a);

These are both constants for every call to reduced2 that ode45 makes. ode45 called reduced2 over 10,000 times and every time it recomputes these constant expressions. So, let's lift them out and just compute them once.

Instead of doing

tic
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);  
oldtime = toc

I do

tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*2:10;
Torque = griddedInterpolant(time, T_a);
[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0,opts); 
newtime = toc

I've created a new function reduced3 that doesn't take T_a. The only reason you needed T_a was to create Torque. I've done that outside the function and sent in in as an argument. I now get the following times.

oldtime =
31.9675
newtime =

0.1238
Speedup is 258.130906x

I've also checked that all results are equal. Here is my reduced3

function [yp] = reduced3(t,y,Torque)  
beta = 13*(pi/180);                      %Helix angle (rad)
P = 20*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;
R_a = 51.19e-3;                             %Radius
R_p = 139.763e-3;                           %Radius
% Common parameters
e_a = (pi*2*R_a*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
e_p = (pi*2*R_p*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
% e = 0;
e = (pi*2*(R_a+R_p)*tan(beta))/(4*cos(P));
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
Cs = 0.1 + 0.0001*sin(2*pi*Freq*t);                     %Shear damping
% Driver gear  
m_a = 0.5;                      %mass of driver gear (kg)
c_ay=500;                       %Damping of driver gear in y direction (Ns/m)
c_az=500;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 8e7;                      %Stiffness of driver gear in y direction (N/m)
k_az= 5e7;                      %Stiffness of driver gear in z direction (N/m)
I_a = 0.5*m_a*(R_a^2);          %Moment of inertia of driver gear (Kg.m3)
% y_ac= e_a + theta_a_vec*R_a;    %circumferential displacement at the meshing point on the driver gear (m)
% y_pc = e_p - theta_a_vec*R_a;   %circumferential displacement at the meshing point on the driven gear (m)
z_a = e_a*tan(beta);  
z_p = e_p*tan(beta);  
% z_ac = (z_a-y_ac)*tan(beta);    %axial displacements of the meshing point on the driving gear  (m)
% z_pc = (z_p+y_pc)*tan(beta);    %axial displacements of the meshing point on the driven gear   (m)
% Driven gear  
m_p = 0.5;                                  %mass of driven gear (kg)
c_py=500;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=500;                                   %Damping of driven gear in z direction (Ns/m)
k_py= 9.5e7;                                %Stiffness of driven gear in y direction (N/m)
k_pz =9e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 0.5*m_p*(R_p^2);                      %Moment of inertia of driven gear (Kg.m3)
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio
% % Excitation forces 
% Fy=ks*cos(beta)*(y_ac-y_pc-e) + Cs*cos(beta)*2*R_a*speed*2*pi;                %axial dynamic excitation force at the meshing point (N)
% Fz=ks*sin(beta)*(z_ac-z_pc-z) - Cs*sin(beta)*2*tan(beta)*R_a*speed*2*pi;      %circumferential dynamic excitation force at the meshing point (N)
%Time interpolation  
Torque = Torque(t)/1000;
% torque_table = [time', T_a];
% torque = interp1(torque_table(:,1), torque_table(:,2), t, 'linear', 'extrap') / 1000;
Opp_Torque = 68.853/1000;                        % Average torque value
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
ks = 0.9e8 + 20000*sin(2*pi*Freq*t);                   %Shear stiffness 
TE = y(11)*R_p - y(5)*R_a;                   %transmission error
b = 20e-6;                                   %Backlash
if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end
   
yp(1) = y(2); 
yp(2) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Cs*cos(beta)*R_a*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_a*(R_a*y(5)+R_p*y(11)) -KS*cos(beta)*R_a*(e_a-e_p-e))/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (Opp_Torque*i - Cs*cos(beta)*R_p*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_p*(R_a*y(5)+ R_p*y(11)) -KS*cos(beta)*R_p*(e_a-e_p-e))/I_p;                % angular accelration (rad/s2) 
end

and here is how we use it with your original version included for comparison

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:2:1000001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*2:10;    %can try 4 or 5
y0 = [0;0;0;0;0;-104.719;0;0;0;0;0;104.719/3]; 
opts = odeset('RelTol',1e-1,'AbsTol',1e-1);
tic
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);  
oldtime = toc


tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*2:10;
Torque = griddedInterpolant(time, T_a);
[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0,opts); 
newtime = toc
fprintf("Speedup is %fx\n",oldtime/newtime)

fprintf("Are all results equal?\n")
all(tnew==t)
all(ynew==y,'all')

% Driver gear graphs
nexttile  
plot(t,y(:,2))
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in y direction vs time')
axis padded  
grid on
nexttile  
plot(t,y(:,4))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,6))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driver gear angular velocity vs time')  
axis padded  
grid on
% Driven gear graphs
nexttile  
plot(t,y(:,8))  
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in y direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,10))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,12))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driven gear angular velocity vs time')  
axis padded  
grid on

[Siddharth Jain]

Hi Mike, thank you so much for spending time to solve my problem. I am trying to run the corrrections you suggested but I am getting the follwing errors. Aplogies if I am being dumb, matlab is not my forte.

Error using griddedInterpolant
Sample points vector corresponding to grid dimension 1 must contain 500001 elements.

Error in reduced2 (line 61)
Torque = griddedInterpolant(time, T_a);

Error in figinput_reduced3>@(t,y)reduced2(t,y,T_a) (line 9)
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in figinput_reduced3 (line 9)
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);
 

[Mike Croucher]

In the first instance, I suggest using the exact things I posted. I have a new function, reduced3. It needs to be called differently from your reduced2. Copy and paste the full functions I posted.

[Siddharth Jain]

Apologies, I had made some changes to reduced2 which caused the error, however I do not want to use the opts to chnage the tolerance since your modifications improve the speed.

The actual graphs look like this and the ode options make it something different.

[Siddharth Jain]

Also, is it possible to run the function for the improved time calculation only, so running only reduced 3

[Mike Croucher]

Sure. Just don't use your reduced2 and use reduced3 only. Can delete opts or change it to whatever you want.

[Siddharth Jain]

Thanks for that. However, I am getting different results from yourself as I removed ode options of reltol and abstol as they did not help in getting the correct results.

Please let me clarify my query. Is it possible to run the modifications you provided without ode opts and without comparing the simulation times between reduced2 and reduced 3 so that only the improved version of reduced 3 gives the results (i.e if reduced 2 and reduced 3 are run together as your code suggests, I still have to wait a long time)

[Mike Croucher]

Yes.

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:2:1000001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*2:10;    %can try 4 or 5
y0 = [0;0;0;0;0;-104.719;0;0;0;0;0;104.719/3]; 



tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*2:10;
Torque = griddedInterpolant(time, T_a);
[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0); 
newtime = toc

% Driver gear graphs
nexttile  
plot(t,y(:,2))
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in y direction vs time')
axis padded  
grid on
nexttile  
plot(t,y(:,4))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,6))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driver gear angular velocity vs time')  
axis padded  
grid on
% Driven gear graphs
nexttile  
plot(t,y(:,8))  
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in y direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,10))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,12))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driven gear angular velocity vs time')  
axis padded  
grid on

[Siddharth Jain]

sorry for troubling you again, but when I run this with reduced 3, I obtain the follwing error with a blank graph-

newtime =

33.6263

Unrecognized function or variable 't'.

Error in figinput_reduced3 (line 20)

plot(t,y(:,2))

[Mike Croucher]

Sorry.

Change

[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0); 

to

[t,y] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0); 

[Siddharth Jain]

Hi @Mike Croucher, thank you for all your help. Can you also help me with plotting KS vs time in my code

ks = 0.9e8 + 20000*sin(2*pi*Freq*t);                   %Shear stiffness 
TE = y(11)*R_p - y(5)*R_a;                   %transmission error
b = 20e-6;                                   %Backlash
if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end