Resolver un sistema de ecuación diferencial de retardo (DDE) restringido para proporcionar soluciones no negativas

Question

En MATLAB, ode45 tiene un parámetro llamado NonNegative que restringe las soluciones para que no sean negativas. They even wrote a paper about how this method works y cómo no es algo tan estúpido como simplemente establecer y_i en 0 cada vez que se vuelve negativo, ya que generalmente no funcionará.

Ahora, MATLAB también tiene dde23 para resolver ecuaciones diferenciales de retardo, pero no existe un parámetro NonNegative equivalente para este integrador.

Desafortunadamente, tengo la tarea de agregar un retraso a un sistema ODE existente que se soluciona usando ode45 con NonNegative activado.

¿Alguna idea de cómo debo proceder?

EDIT:

No estoy seguro de si esto es útil, pero ...

La parte DDE de mi sistema básicamente se parece a:

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

donde X (la capital variable de letra en la tercera ecuación) es la versión diferida de x. Luego, estoy vinculando este sistema DDE a un sistema ODE existente (y más grande) agregando un par de términos a las ecuaciones para x y z, y luego integrando el sistema combinado todos juntos.

Answer 1

Tiene un problema difícil y no estoy seguro de si hay una solución de un solo paso. Estaría más que contento de felicitar a cualquiera que esté dispuesto a proporcionar una respuesta alternativa.

Dependiendo de la duración de la demora, una opción sería ejecutar la ecuación varias veces, con cada iteración pasando los valores anteriores de x a la última actualización.

Por ejemplo, supongamos que su retraso es de una hora. En la primera hora, ejecute ode45 con NonNegative marcado. Almacene el valor en una matriz nueva junto con el parámetro de tiempo y vuelva a ejecutar el algoritmo. En esta ocasión, asegúrese de añadir dos parámetros de entrada: su matriz solución antigua y la matriz de antaño

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; 

Ahora lave, enjuague y repita. Tenga en cuenta que X ahora es un término cuasi dependiente del tiempo como se ve en el ejemplo 3 del ode45.

Answer 2

Recientemente tuve este problema con un código mío.La solución 'simple' es hacer lo siguiente, en primer lugar, una vez que la solución llegue a 0, hay que mantenerlo a 0, así cambiar

dx = 1/(1+y*z) - x; 

a (aviso en el que se evalúa el caso x == 0)

if x > 0 
    dx = 1/(1+y*z) - x; 
else % if x <= 0 
    dx = 0; 
end 

o tal vez a (dependiendo de por qué nunca puede ser 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 

Nota sin embargo que esto crea una discontinuidad de salto en la primera derivada, WHI ch cuando el solucionador DDE llega al punto de tener t - delay cerca de este punto, no lo resuelve de manera muy eficiente a menos que sepa el punto exacto donde se encuentra esta discontinuidad (generalmente tiene una opción adicional de decirle a Matlab dónde está, pero si siga los pasos a continuación, entonces eso no será necesario).

Para determinar el lugar de esta discontinuidad es necesario utilizar la opción DDE de events (baje hasta 'Evento ubicación Propiedad', también se puede ver en estos examples, uno de esos ejemplos en realidad se ocupa de un sistema similar en negativo los valores no están permitidos en un ODE - los eventos para ODE y DDE son casi idénticos). Básicamente, un evento es una función de Matlab con una salida vectorial, cada una de las entradas del vector es alguna u otra evaluación de sus variables; en cada paso Matlab comprueba si uno de ellos eqauls 0, cuando uno de ellos es igual a 0, el DDE se detiene y devuelve una solución hasta ese punto, desde el cual debe reiniciar el DDE con esa solución parcial como su historial, es decir, en lugar de ejecutar

sol = dde23(ddefun, lags, history, [t0 tEnd], options); 

que corren (aviso sol y t0 cambiado)

sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); 

En este caso una de las entradas del vector va a ser x (como se desea que el DDE para detener cuando x es igual a 0). Además, la línea del código elseif dxTmp <= 0 crea otra discontinuidad, por lo que necesita un evento para cuando dxTmp se convierta en 0, es decir, 1/(1+y*z) - x será otro componente de la salida del vector.

Ahora cuando reinicia el ODE, Matlab automáticamente asume que hay una discontinuidad en ese punto, por lo tanto, no tiene que preocuparse de decirle a Matlab que hay uno allí.

El siguiente problema es que Matlab nunca logra hacerlo bien, los valores de x, y, z y X será ligeramente negativo. Por lo tanto, si se va a crear un problema, tendrá que corregir el valor de x (y los demás valores de manera similar) con

if x < 0 
    x = 0; 
end 

antes de calcular los derivados. Pero esto solo cambia el valor de x localmente. Por lo tanto, es posible que desee cambiar todos los valores negativos de x en la solución final a 0 también. Sugiero que no intente cambiar sol antes de ingresarlo en sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); ya que hice varios intentos y no pude hacerlo funcionar.

Answer 3

Hay una respuesta más simple , use createOptimProblem para configurar una optimización. Debe incluir límites para cada parámetro con este método, pero resulta trivial forzar a un parámetro a permanecer positivo.

Detalles aquí: http://www.mathworks.com/help/gads/createoptimproblem.html