# Help with ODE tolerance options

4 views (last 30 days)
Siddharth Jain on 15 Mar 2023
Edited: Siddharth Jain on 16 Mar 2023
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
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')
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')
grid on
nexttile
plot(t,y(:,6))
xlabel('Time')
title('Driver gear angular velocity vs time')
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')
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')
grid on
nexttile
plot(t,y(:,12))
xlabel('Time')
title('Driven gear angular velocity vs time')
grid on
% nexttile
% plot(t,T_a(1:60001))
% ylabel('Torque (Nm)')
% xlabel('Time (sec)')
% title('Torque vs time')
% grid on
Jan on 15 Mar 2023
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 on 15 Mar 2023
I didn't notice you had already suggested that @Jan - sorry. That's pretty much all my answer below is about.

Mike Croucher on 15 Mar 2023
Edited: Mike Croucher on 15 Mar 2023
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.
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
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')
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')
grid on
nexttile
plot(t,y(:,6))
xlabel('Time')
title('Driver gear angular velocity vs time')
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')
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')
grid on
nexttile
plot(t,y(:,12))
xlabel('Time')
title('Driven gear angular velocity vs time')
grid on
Mike Croucher on 15 Mar 2023
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 on 16 Mar 2023
Edited: Siddharth Jain on 16 Mar 2023
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

### Categories

Find more on Gears in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!