Time-varying coefficient in ODE.
조회 수: 4 (최근 30일)
이전 댓글 표시
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.
댓글 수: 0
채택된 답변
Jan
2021년 10월 2일
편집: Jan
2021년 10월 2일
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
Jan
2021년 10월 2일
편집: Jan
2021년 10월 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
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Mathematics and Optimization에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!