Time-varying coefficient in ODE.

Hello everyone! I am developing a nonlinear dynamic model of a gearbox. The gear mesh stiffness between the gears is a function of shaft position. So, on each iteration step, I must evaluate the value of gear mesh stiffness for the ODE solver.

P.s. I used two approaches:

1. To calculate the gear mesh stiffness and save the values in a look-up table. Then on each iteration step in the ODE solver simply interpolate values from the look-up table. As I have found this approach is not the best, because the interp1 function is not optimal, and slows down the calculation process significantly.

2. Another approach is to make symbolic Fourier series outside of the ODE solver and represent this series as a function handle. Then this function handle is declared as a global variable. So, on each iteration step, the gear mesh stiffness is evaluated in the ODE solver. In my understanding, this approach should be faster, but it is not.

Both methods which I am using require a lot of time for simulation. So, I am trying to find the optimal way, how to simulate my dynamic model. Any suggestions highly appreciated.

Answer 1

[Jan]

Which ODE solver do you use? Remember that Matlab's builtin ODE solvers like ODE45 ar designed for smooth functions only. A parameter, which is obtained by a linear interpolation is not smooth, such that the stepsize controller of the integrator has an undefined behavior. Do not use this for scientific work.

Global variables are a mess. Use anonymous functions instead: Anonymous function for parameters

The computing time critically depends on the tolerances. Did you check with the profiler, if the interpolation is really the bottleneck? If this is the case:

You can use a faster implementation of the interpolation, e.g. https://www.mathworks.com/matlabcentral/fileexchange/25463-scaletime or:

x  = linspace(0, 2*pi, 100);
y  = sin(x);
xi = linspace(0, 2*pi, 1000);

tic
for k = 1:1e3
yi = Interp1_eqx(x, y, xi);
end
toc

Elapsed time is 0.064923 seconds.

tic
for k = 1:1e3
yi = interp1(x, y, xi, 'linear');
end
toc

Elapsed time is 0.249932 seconds.

function F = Interp1_eqx(x, y, xi)
% Linear interpolation.
% INPUT: x, y, xi: Vectors. x must have equidistant steps.
%        xi must inside the range of x.
% (C) Jan, Heidelberg, CC BY-SA 3.0
ny = numel(y);
u  = 1 + (xi - x(1)) / (x(ny) - x(1)) * (ny - 1);  % Normalize

v    = floor(u);
s    = u - v;         % Interpolation parameters
d    = (v == ny);     % Elements on the boundary
v(d) = v(d) - 1;
s(d) = 1;

F = y(v) .* (1 - s) + y(v + 1) .* s;  % Interpolate
end

Check the speed locally. The timings in Matlab online are not accurate.

[Jan]

An approach, which accepts matrices as input y:

function F = Interp1_eqx2D(x, y, xi)
% Linear interpolation.
% INPUT: x, xi: Vectors. x must have equidistant steps.
%               xi must inside the range of x.
%        y: column vector or matrix, operate on 1st dim!
% (C) Jan, Heidelberg, CC BY-SA 3.0

ny = size(y, 1);
u  = 1 + (xi(:) - x(1)) / (x(ny) - x(1)) * (ny - 1);  % Normalize

% #Uncomment if xi can be outside x:
% uout    = (u < 1) | (u > nrows);  % Outside the range of x
% u(uout) = 1;

v    = floor(u);
s    = u - v;         % Interpolation parameters

d    = (v == ny);     % Elements on the boundary
v(d) = v(d) - 1;
s(d) = 1;

F = y(v, :) .* (1 - s) + y(v + 1, :) .* s;  % Interpolate

% #Uncomment if xi can be outside x:
% F(uout) = NaN;        % Set out of range values to NaN
end

[Iurii Storozhenko]

Jan, thank you very much for your prompt response and valuable suggestions. I have tried different ODE solvers, but the fastest is ODE15s for my case. I will play around with the profiler and will get back to you! Thanks a lot for your help!