Пользовательская периодическая функция без счетчика

Это распространенный тип вопроса. Пользователи часто хотят периодически изменять переменную в функции интегрирования, используемой решателями ODE Matlab. К сожалению, обычно это плохая идея. Эти решатели предполагают, что дифференциальные уравнения являются гладкими и непрерывными. В лучшем случае ваш код будет медленнее. В худшем случае у вас будут большие ошибки или даже совершенно неправильные результаты. Использование жесткого решателя, такого как ode15s, может немного помочь, но это слишком предполагает преемственность. Поскольку ваша система, как указано, имеет аналитическое решение, самый простой способ «симулировать» это может состоять в том, чтобы пропустить последовательность импульсов через эту функцию времени.

Случай, когда параметр изменяется скачкообразно во времени, т. е. по независимой переменной, решается легко. Просто интегрируйте по каждому временному интервалу необходимое количество периодов:

t1 = 600;
t2 = 120;
TSpan = [0 t1]; % Initial time vector
NPeriods = 5;   % Number of periods
C0 = 400;       % Initial condition
Q1 = ...
Q2 = ...
V = 100;
C_i = 400;
S = 100;
u_vent = @(t,C,Q)(Q*(C_i-C(1))+S)/V; % Integration function as anonymous function
Cout = C0;       % Array to store solution states
tout = Tspan(1); % Array to store solution times
for i = 1:NPeriods
    [t,C] = ode45(@(t,C)u_vent(t,C,Q1), TSpan, C0);
    tout = [tout;t(2:end)];   % Append time, 2:end used to avoid avoid duplicates
    Cout = [Cout;C(2:end,:)]; % Append state
    TSpan = [0 t2]+t(end);    % New time vector
    C0 = C(end);              % New initial conditions
    [t,C] = ode45(@(t,C)u_vent(t,C,Q2), TSpan, C0);
    tout = [tout;t(2:end)];
    Cout = [Cout;C(2:end,:)];
    TSpan = [0 t1]+t(end);
    C0 = C(end);
end

Это позволяет ode45 интегрировать ОДУ за указанное время из набора начальных условий. Затем интегрирование перезапускается в момент окончания предыдущего интегрирования с новыми разрывными начальными условиями или другими параметрами. Это приводит к переходному процессу. Вам может не понравиться внешний вид кода, но вот как это делается. Еще сложнее, если параметры изменяются дискретно по отношению к переменным состояния.

Необязательный. Каждый раз, когда вы вызываете/перезапускаете ode45, функция должна определить начальный размер шага для использования. Это основные (минимальные) затраты/накладные расходы на перезапуск решателя. В некоторых случаях вы можете помочь решателю, указав 'InitialStep'< /a> параметр, основанный на последнем размере шага, использованном в предыдущем вызове. См. демонстрацию ballode для получения дополнительной информации об этом, набрав edit ballode в командном окне.

Если у вас есть набор инструментов обработки сигналов, вы можете использовать что-то вроде:

>> t = 0:3600;
>> y = pulstran(t,[660:720:3600],'rectpuls',120);
>> plot(t,y)
>> ylim([-0.1 1.1])

что дает следующее (в Octave должно быть то же самое в MATLAB):

Затем вам нужно масштабировать y, чтобы оно находилось между Q1 и Q2 вместо 0 и 1.

Не обязательно лучший подход, потому что предположения о непрерывности остаются нарушенными, но способ генерировать прямоугольную последовательность импульсов без цепочки if:

Q = Q2 + (Q1 - Q2) * (mod(t, period) < t_rise);

который работает как в скалярном, так и в векторном контексте.