У меня недавно была эта проблема с некоторым моим кодом. «Самое простое» решение состоит в том, чтобы сделать следующее: во-первых, как только решение достигнет 0, вы должны оставить его равным 0,
таким образом изменить
dx = 1/(1+y*z) - x;
до (обратите внимание, где оценивается дело x == 0)
if x > 0
dx = 1/(1+y*z) - x;
else % if x <= 0
dx = 0;
end
или, может быть (в зависимости от того, почему оно никогда не может быть 0)
dxTmp = 1/(1+y*z) - x;
if x > 0
dx = dxTmp;
elseif dxTmp > 0
% x becomes positive again
dx = dxTmp;
else
dx = 0;
end
Заметим, однако, что это создает скачок скачка в первой производной, который, когда решатель DDE достигает точки, имеющей t - delay вблизи этой точки, не решает ее очень эффективно, если не знает точного места, где находится этот разрыв (обычно вы используете дополнительную опцию, чтобы указать Matlab, где он находится, но если вы выполните следующие шаги, это не понадобится).
Чтобы определить место этого разрыва, вам нужно использовать опцию DDE events (прокрутите вниз до «Property Location Property», вы также можете посмотреть на эти examples , один из этих примеров фактически имеет дело с подобной системой, где отрицательные значения недопустимы в ODE - события для ODE и DDE почти идентичны). По сути, событие - это функция Matlab с выходным вектором, каждая из записей вектора является той или иной оценкой ваших переменных; на каждом шаге Matlab проверяет, равен ли один из них 0, когда один из них равен 0, DDE останавливается и возвращает решение до той точки, с которой вы должны перезапустить DDE с этим частичным решением в качестве вашей истории, т.е. вместо выполнения
sol = dde23(ddefun, lags, history, [t0 tEnd], options);
вы запускаете (уведомление sol и t0 изменено)
sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options);
В этом случае одна из записей вектора будет x (так как вы хотите, чтобы DDE остановился, когда x равно 0). Также строка кода elseif dxTmp <= 0 создает еще один разрыв, поэтому вам нужно событие, когда dxTmp становится 0, т. Е. 1/(1+y*z) - x будет еще одним компонентом векторного вывода.
Теперь, когда вы перезапускаете ODE, Matlab автоматически предполагает, что в этот момент есть разрыв, поэтому вам не нужно беспокоиться о том, чтобы сообщить Matlab, что он там есть.
Следующая проблема состоит в том, что Matlab никогда не достигает этого правильно, значения x, y, z и X будут слегка отрицательными. Таким образом, если это создаст проблему, вам нужно исправить значение x (и другие значения аналогично) с помощью
if x < 0
x = 0;
end
до расчета производных. Но это только меняет значение x локально. Поэтому вы можете изменить все отрицательные значения x в решении final на 0. Я предлагаю вам не пытаться изменить sol перед вводом его в sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options);, так как я сделал несколько попыток и не смог заставить его работать.