Time-varying coefficient in ODE.
6 views (last 30 days)
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.
Jan on 2 Oct 2021
Edited: Jan on 2 Oct 2021
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);
for k = 1:1e3
yi = Interp1_eqx(x, y, xi);
for k = 1:1e3
yi = interp1(x, y, xi, 'linear');
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
Check the speed locally. The timings in Matlab online are not accurate.