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.
采纳的回答
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.
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
tic
for k = 1:1e3
yi = interp1(x, y, xi, 'linear');
end
toc
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.
2 个评论
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
2021-10-3
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
