Решение системы дифференциальных уравнений с запаздыванием (DDE), ограниченной для получения неотрицательных решений

Question

В MATLAB, ode45 имеет параметр с именем NonNegative что ограничивает решения быть неотрицательными. Они даже написали статью о том, как работает этот метод, и о том, что он не такой глупый, как просто установка y_i в 0, когда он становится отрицательным, так как это обычно не работает.

Теперь, MATLAB также имеет dde23 для решения дифференциальных уравнений с запаздыванием, но нет эквивалента NonNegative Параметр для этого интегратора.

К сожалению, мне поручено добавить задержку в существующую систему ODE, которая решается с помощью ode45 с NonNegative включен.

Есть идеи, как мне поступить?

РЕДАКТИРОВАТЬ:

Я не уверен, что это полезно, но...

Часть DDE моей системы в основном выглядит так:

 dx = 1/(1+y*z) - x;
 dy = (y*z)^2/(1+(y*z)^2) - y;
 dz = X - z;

где X (переменная заглавной буквы в третьем уравнении) является отложенной версией x, Затем я связываю эту систему DDE с существующей (и более крупной) системой ODE, добавляя пару терминов в уравнения для x а также z и затем объединяя объединенную систему все вместе.

Answer 1

У меня недавно была эта проблема с некоторым моим кодом. "Самое простое" решение состоит в том, чтобы сделать следующее: во-первых, как только решение достигнет 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 (прокрутите вниз до "Свойства местоположения события", вы также можете посмотреть на эти примеры, один из этих примеров фактически имеет дело с аналогичной системой, где отрицательные значения недопустимы в 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 в окончательном решении до 0, а также. Я предлагаю вам не пытаться изменить sol перед вводом в sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); так как я сделал несколько попыток, я не смог заставить его работать.

Answer 2

У вас сложная проблема, и я не уверен, есть ли одношаговое решение. Я был бы более чем рад отдать должное любому, кто хочет дать альтернативный ответ.

В зависимости от длительности задержки, один из вариантов будет запускать уравнение несколько раз, при этом каждая итерация передает старые значения x до последнего обновления.

Например, скажем, ваша задержка составляет один час. В первый час запустите ode45 с пометкой "Неотрицательный". Сохраните значение в новую матрицу вместе с параметром Time и снова запустите алгоритм. На этот раз убедитесь, что вы добавили два входных параметра: матрицу старого решения и матрицу старого времени.

dx = 1/(1+y*z) - x; 
dy = (y*z)^2/(1+(y*z)^2) - y;
tindex = find(told>t,1) -1 % find the upper index which best approximates t
X = xold(tindex) + (xold(tindex+1)-xold(tindex))*(t-told(tindex))/(told(tindex+1)-told(tindex)) % or interpolation method of your choosing
dz = X - z;

Теперь вымойте, сполосните и повторите. Обратите внимание, что теперь X является квази-зависящим от времени термином, как видно в примере 3 из ode45.

Answer 3

Там гораздо более простой ответ, используйте createOptimProblem настроить оптимизацию. Вы должны включить границы для каждого параметра, использующего этот метод, но становится тривиально заставить параметр оставаться положительным.

Подробности здесь: http://www.mathworks.com/help/gads/createoptimproblem.html