ODE с временным входом, как ускорить Up без использования интерполяции?

Question

Я пытаюсь решить систему ODE, и мое входное возбуждение является функцией времени.

Я использовал interp1 внутри функции интеграции, но это не похоже на очень эффективный способ сделать это. Я знаю, что это не так, потому что как только я изменю входное возбуждение на функцию sin, которая не требует вызова interp1 внутри функции, я получаю гораздо более быстрые результаты. Но делая интерполяцию, каждый шаг занимает примерно 10-20 раз дольше, чтобы сходиться. Итак, есть ли лучший способ решения ODE для произвольного зависящего от времени возбуждения, без необходимости делать интерполяцию или некоторые другие трюки для ускорения?

Я просто скопировать модифицированную версию simple example from The MathWorks здесь:

входного возбуждением является gradually increasing sin функцией, но после того, как через некоторое время он становится constant amplitude sin функции.

Dt = 0.01; % sampling time step 
Amp0 = 2;  % Final Amplitude of signal 
Dur_G = 10; % Duration of gradually increasing part of signal 
Dur_tot = 25; % Duration of total signal 
t_G = 0 : Dt : Dur_G; % time of gradual part 
A = linspace(0, Amp0, length(t_G)); 
carrier_1 = sin(5*t_G); % Unit Normal Signal 
carrier_A0 = Amp0*sin(5*t_G); 
out_G = A.*carrier_1; % Gradually Increasing Signal 
% Total Signal with Gradual Constant Amplitude Parts 
t_C = Dur_G+Dt:Dt:Dur_tot; % time of constant part 
out_C = Amp0*sin(5*t_C); % Signal of constant part 
ft = [t_G t_C]; % total time 
f = [out_G out_C]; % total signal 
figure; plot(ft, f, '-b'); % input excitation 


function dydt = myode(t,y,ft,f) 
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t 
g = 2;    % a constant 
dydt = -f.*y + g; % Evaluate ODE at time t 

tspan = [1 5]; ic = 1; 
opts = odeset('RelTol',1e-2,'AbsTol',1e-4); 
[t,y] = ode45(@(t,y) myode(t,y,ft,f), tspan, ic, opts); 

figure; 
plot(t,y); 

Обратите внимание, что я объяснил только первую часть моей проблемы выше, которая решает систему постепенно возрастающей sin функции. Во второй части мне нужно решить ее для произвольного входного возбуждения (например, входного сигнала ускорения заземления).

Answer 1

Для этого примера, вы можете использовать griddedInterpolant класс, чтобы получить немного ускорения:

ft = linspace(0,5,25); 
f = ft.^2 - ft - 3; 
Fp = griddedInterpolant(ft,f); 
gt = linspace(1,6,25); 
g = 3*sin(gt-0.25); 
Gp = griddedInterpolant(gt,g); 

tspan = [1 5]; 
ic = 1; 
opts = odeset('RelTol',1e-2,'AbsTol',1e-4); 
[t,y] = ode45(@(t,y)myode(t,y,Fp,Gp),tspan,ic,opts); 

figure; 
plot(t,y); 

Функция ОДА затем:

function dydt = myode(t,y,Fp,Gp) 
f = Fp(t);  % Interpolate the data set (ft,f) at time t 
g = Gp(t);  % Interpolate the data set (gt,g) at time t 
dydt = -f.*y + g; % Evaluate ODE at time t 

На моей системе с R2015b, то звонок в ode45 примерно в три раза быстрее (0,011 секунды против 0,035 с) для вашего примера. Вы можете получить немного большую скорость, переключившись на ode23. Вы можете узнать больше о классе griddedInterpolanthere.

Если ваша фактическая система, дискретно переключается между входами в определенные моменты времени, то вы, вероятно, должны решить проблему кусочно, интегрируя каждый случай отдельно. См. this question и this question. Если система переключается на основе значения переменных состояния, вы должны использовать event location (см. this question). Однако, если «решение ODE для случайного зависящего от времени возбуждения» означает, что вы добавляете случайный шум в систему, тогда у вас есть SDE, а не ODE, который является совершенно другим зверем.