# Help with ODE tolerance options

4 views (last 30 days)

Show older comments

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

##### 4 Comments

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

### Answers (1)

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.

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

##### 10 Comments

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);

### See Also

### Categories

### Products

### Community Treasure Hunt

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

Start Hunting!