Sistemas de equações diferenciais com atrasos fixos - SIRS

Anterior: Exemplo: equação diferencial com retardo | Índice: Ecologia | Próximo: Simulação e modelo de campo médio

O sistema

O exemplo a ser considerado é da biologia matemática, muitas ferramentas matemáticas utilizado na dinâmica de populações encontra-se também na dinâmica de epidemias. O objetivo deste tópico não é se aprofundar na modelagem em si, mas a utilizar como motivação para compreender melhor os sistemas de equações diferenciais com atraso baseado no artigo "Oscillations in SIRS model with distributed delays".

Uma doença do tipo SIRS é caracterizada por um tempo de infecção ${\tau }_i$ e um tempo de imunidade ${\tau }_r$. Isto é, um indivíduo que se infecte em um tempo $t$ vai deterministicamente se recuperar em um tempo $t+{\tau }_i$, tornando-se temporariamente imune por um tempo ${\tau }_r$. Então no instante $t+{\tau }_i+{\tau }_r\equiv t+{\tau }_0$ perde a imunidade tornando-se suscetível novamente. Este sistema pode ser representado por um conjunto de equações para a fração de população infecciosa (e infectada) $i\left(t\right)$ e suscetível $s\left(t\right)$:

$$\dot{s}\left(t\right)=-\beta s\left(t\right)i\left(t\right)+\beta s(t-{\tau }_0)i(t-{\tau }_0)$$

$$\dot{i}\left(t\right)=\beta s\left(t\right)i\left(t\right)-\beta s(t-{\tau }_i)i(t-{\tau }_i)$$

Onde $\beta$ é a taxa de contágio por indivíduo. A fração da subpopulação em recuperação (imune) é dada por $r\left(t\right)=1-s\left(t\right)-i\left(t\right)$. O primeiro termos em ambas as equações representa o contágio dos suscetíveis pelos infecciosos, daqui em diante muitas vezes será referido como: $c\left(t\right)=\beta s\left(t\right)i\left(t\right)$. Os segundos termos representam a recuperação ou a perda de imunidade após ${\tau }_i$ e ${\tau }_0$ decorrido desde a infecção, respectivamente. Isto é:

$$\dot{s}\left(t\right)=-c\left(t\right)+c(t-{\tau }_0)$$

$$\dot{i}\left(t\right)=c\left(t\right)-c(t-{\tau }_i)$$

Em outras palavras, uma fração $c\left(t\right)$ se contagia no instante $t$, no instante $t+{\tau }_i$ ela se recupera, diminuindo os infecciosos. No instante $t+{\tau }_i+{\tau }_r$ ela perde a imunidade, aumentando a quantidade de suscetíveis. Estas equações exigem condições iniciais, matematicamente usualmente se fornece funções arbitrárias $s\left(t\right)$ e $i\left(t\right)$ no intervalo $[-{\tau }_0,0]$. Porém de um ponto de vista epidemiológico, é razoável prover condições iniciais em $t=0$ e dinâmicas complementares nos intervalos $[0,{\tau }_i)$ e $[{\tau }_i,{\tau }_0)$ .

• $[0,{\tau }_i)\to$ somente contágio local. Utiliza-se as equações apenas com os primeiros termos, sem atraso.
• $[{\tau }_i,{\tau }_0)\to$Transição do estágio de infeccioso para em recuperação, ou seja, retira-se apenas o segundo termo na equação da dinâmica dos suscetíveis.

E com as condições iniciais:

$$s\left(0\right)=1-i_0,i\left(0\right)=i_0,r\left(0\right)=0$$

Ponto de equilíbrio

Para análise dos pontos de equilíbrio uma representação integral é uma forma melhor forma de representar o sistema, pois se fizermos algo análogo ao que foi feito em Equações diferenciais com atrasos, não será possível obter nenhuma informação sobre o ponto de equilíbrio, pois aparentemente qualquer par de valores poderia representar um ponto de equilíbrio. Para a dinâmica dos infecciosos é proposto então:

$$i\left(t\right)=c_1+\underset{t-{\tau }_i}{\overset{t}{\int }}c\left(u\right)du$$

A interpretação é direta, integra-se sobre sobre todos os indivíduos que se contagiaram entre o tempo $t-{\tau }_i$ e $t$. Estes serão os infecciosos no instante $t$, uma vez que todos infectados anteriormente a este tempo já estarão recuperados. A constante de integração $c_1$ a princípio é arbitrária, mas espera-se que $c_1=0$ pois não deve haver outras fontes de infecções adicionais. Complementariamente:

$$1-r\left(t\right)=c_2-\underset{t-{\tau }_0}{\overset{t-{\tau }_i}{\int }}c\left(u\right)du$$

Pelo lado esquerdo, pode-se perceber ver que isto se refere a população que não está se recuperando, isto é $1-r\left(t\right)=s\left(t\right)+i\left(t\right)$, ou seja, os infectados e suscetíveis. A integral cobre os contágios que ocorreram durante ${\tau }_r$ anterior. Desta forma, esta integral fornece o número de pessoas que estão no período de recuperação, logo é razoável supor que $c_2=1$ uma vez que também não há outras fontes de pessoas em recuperação. Para um estado de equilíbrio então que a fração de contagiados em um instante $t$ qualquer é constante, pode-se escrever $c\left(u\right)=c_0$ . Para $i\left(t_0\right)$:

$$\begin{array}{rl}i\left(t_0\right)& =c_0\underset{t_0-{\tau }_i}{\overset{t_0}{\int }}du\\ & =c_0(t_0-t_0+{\tau }_i)\\ & =c_0{\tau }_i\end{array}$$

E para $s\left(t_0\right)$, de maneira análoga:

$$\begin{array}{rl}s\left(t_0\right)& =1-\underset{t_0-{\tau }_0}{\overset{t_0-{\tau }_i}{\int }}c\left(u\right)du-i\left(t_0\right)\\ & =1-c_0\underset{t_0-{\tau }_0}{\overset{t_0-{\tau }_i}{\int }}du-c_0{\tau }_i\\ & =1-c_0(t_0-{\tau }_i-t_0+{\tau }_0)-c_0{\tau }_i\\ & =1-c_0(-{\tau }_i+{\tau }_i+{\tau }_r)-c_0{\tau }_i\\ & =1-c_0{\tau }_r-c_0{\tau }_i\\ & =1-c_0({\tau }_r+{\tau }_i)\\ & =1-c_0{\tau }_0\end{array}$$

Uma vez que $s\left(t\right)=1-r\left(t\right)-i\left(t\right)$. Lembrando que $c\left(t\right)=\beta s\left(t\right)i\left(t\right)$ , então no estado de equilíbrio $c_0=\beta s\left(t_0\right)i\left(t_0\right)$, denotando os pontos de equilíbrio apenas como $s\left(t_0\right)=s^{*}$ e $i\left(t_0\right)=i^{*}$, obtém-se do primeiro resultado:

$$\begin{array}{rl}{i}^{*}& =\beta s^{*}{i}^{*}{\tau }_i\\ s^{*}& =\frac{1}{\beta {\tau }_i}\end{array}$$

E do segundo:

$$\begin{array}{rl}1& =1-\beta s^{*}{i}^{*}{\tau }_0\\ 1& =\frac{1}{s^{*}}-\beta i^{*}{\tau }_0\\ \beta i^{*}{\tau }_0& =\beta {\tau }_i-1\\ i^{*}& =\frac{\beta {\tau }_i-1}{\beta {\tau }_0}\end{array}$$

Assim o ponto de equilíbrio é: $(s^{*},i^{*})=(\frac{1}{\beta {\tau }_i},\frac{\beta {\tau }_i-1}{\beta {\tau }_0})$.

Validade da aproximação e sistema de equações sem atraso

Uma atenção especial deve ser dada a consideração de que $c_1=0$ quando foi proposto a equação integral para $i\left(t\right)$. A formulação matemática do modelo sem atraso é usualmente escrito como:

$$\begin{array}{rl}\frac{ds\left(t\right)}{dt}=& -\beta s\left(t\right)i\left(t\right)+\frac{r\left(r\right)}{{\tau }_r}\\ \frac{di\left(t\right)}{dt}=& \beta s\left(t\right)i\left(t\right)-\frac{i\left(r\right)}{{\tau }_i}\\ \frac{dr\left(t\right)}{dt}=& \frac{i\left(r\right)}{{\tau }_i}-\frac{r\left(r\right)}{{\tau }_r}\end{array}$$

Integrando a equação da dinâmica de $i\left(t\right)$ do sistema sem atraso, de $t^{*}-{\tau }_i$ a $t^{*}$, obtém-se:

$$\begin{array}{rl}\frac{di\left(t\right)}{dt}& =\beta s\left(t\right)i\left(t\right)-\frac{i\left(t\right)}{{\tau }_i}\\ \underset{i(t^{*}-{\tau }_i)}{\overset{i\left(t^{*}\right)}{\int }}di& =\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt-\frac{1}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}i\left(t\right)dt\\ i\left(t^{*}\right)& =\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt+[i(t^{*}-{\tau }_i)-\frac{1}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}i\left(t\right)dt]eq:tri\end{array}$$

Considerando no equilíbrio que $i\left(t^{*}\right)=i(t^{*}-{\tau }_i)$:

$$\begin{array}{rl}0& =\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt-\frac{i^{*}}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}i\left(t\right)dt\\ \frac{i^{*}}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}i\left(t\right)dt& =\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt\end{array}$$

Sendo que no equilíbrio os valores são constantes $i\left(t\right)=i^{*}$ e $s\left(t\right)=s^{*}$:

$$\begin{array}{rl}\frac{i^{*}}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}dt& =\beta s^{*}{i}^{*}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}dt\\ i^{*}& =\beta s^{*}{i}^{*}{\tau }_{i}eq:siete\end{array}$$

De onde é possível obter também o mesmo ponto de equilíbrio obtido anteriormente:

$$\begin{array}{rl}{s}^{*}& =\frac{1}{\beta {\tau }_i}\end{array}$$

Outra formar de comparar as equações, é que para obter versão integral apresentada para $i\left(t\right)$ a partir de::

$$i\left(t^{*}\right)=\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt+[i(t^{*}-{\tau }_i)-\frac{1}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}i\left(t\right)dt]$$

É é necessário que o termo entre colchete seja zerado. Isto é, considerando que no equilíbrio $i(t^{*}-{\tau }_i)$ e $i(t^{*}-{\tau }_i)$ sejam constantes:

$$\begin{array}{rl}i(t^{*}-{\tau }_i)& =\frac{1}{{\tau }_i}\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}i\left(t\right)dt\\ i(t^{*}-{\tau }_i)& =i\left(t^{*}\right)\end{array}$$

Desta forma a condição imposta anteriormente é recuperada. Agora percebe-se que o ponto de equilíbrio é independente das condições iniciais, depende apenas dos parâmetros escolhidos $\beta ,{\tau }_i$ e ${\tau }_r$. Por exemplo, resolvendo numericamente o sistema de equações diferenciais proposto, com os parâmetros $\beta =0.4$, ${\tau }_i=5$ e ${\tau }_r=50$, e condições iniciais $(s_0,i_0)=(0.5,0.5)$, o sistema atinge o equilíbrio precisamente nos pontos calculados.

Solução numérica do sistema de equações diferenciais sem atraso obtido via Mathematica.

Além disto, pode-se verificar que $i^{*}\approx 0.45$ e $c^{*}=\beta s^{*}{i}^{*}{\tau }_i\approx 0.45$, concordando com nossos cálculos analíticos. Isto implica que no equilíbrio, a população total infecciosa em um instante $t$ é sempre dada inteiramente pelos que foram infectados durante um período anterior $t-{\tau }_i$, ou seja, a área dada pelo retângulo de largura ${\tau }_i$ e altura $\beta s^{*}{i}^{*}$, uma vez que no equilíbrio as frações $i\left(t\right)$ e $s\left(t\right)$ são constantes. Estes resultados concordam com a suposição de que $c_1=0$.

Porém tentando o mesmo procedimento para a equação com atraso, obtém-se:

$$\begin{array}{rl}\frac{di\left(t\right)}{dt}=& \beta s\left(t\right)i\left(t\right)-\beta s(t-{\tau }_i)i(t-{\tau }_i)\\ \underset{i(t^{*}-{\tau }_i)}{\overset{i\left(t^{*}\right)}{\int }}di=& \underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt-\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s(t-{\tau }_i)i(t-{\tau }_i)dt\\ i\left(t^{*}\right)=& \underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt+[i(t^{*}-{\tau }_i)-\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s(t-{\tau }_i)i(t-{\tau }_i)dt]eq:cat\end{array}$$

Utilizando a consideração $i\left(t^{*}\right)=i(t^{*}-{\tau }_i)$, é possível obter apenas $i\left(t^{*}\right)=i\left(t^{*}\right)$. Para recuperar a equação original seria necessário que, de modo análogo a equação sem atraso, a seguinte igualdade fosse válida:

$$i\left(t^{*}\right)=\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt$$

Isto é, de modo análogo ao caso anterior, que no equilíbrio a fração de infecciosos em um tempo $t$ seja dado inteiramente pelos que foram contagiados em um tempo $t-{\tau }_i$. Porém neste caso qualquer par de pontos $(i^{*},s^{*})$ poderia consistir em um ponto de equilíbrio para um conjunto de parâmetros $\{\beta ,{\tau }_i,{\tau }_r\}$ qualquer. A dependência se torna das condições iniciais. Utilizando por exemplo condições iniciais constantes como frequentemente é utilizado na literatura de equações com atraso, a própria condição inicial se torna um ponto de equilíbrio $i_0=i^{*}$. Pode-se pensar que desta fração $i^{*}$ de infecciosos , a cada instante $t$, uma quantidade de $c=\beta i_0{s}_0$ de infecciosos são curados, porém outra quantidade igual $c$ de suscetíveis são infectados, mantendo a taxa de variação zerada e a quantidade total de infecciosos constante. Assim sendo, pensando em situações reais, para um instante $t$, há quantidade $i_c=i_0-c$ da fração que permanece infecciosa mesmo após ${\tau }_i$, que soma-se aos novos infectados, mantendo a fração constante. Isto é, para a equação:

$$i\left(t\right)=c_1+\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt$$

Então $c_1\ne 0$. Para retirar esta constante, pode-se pensar que é como se no equilíbrio precisasse de $$n=\frac{i_0}{c}=\frac{i_0}{\beta i_0{s}_0}=\frac{1}{\beta s_0}$$ tempo para que uma pessoa infecciosa se curasse, ao invés de ${\tau }_i$. Pois após $n$ tempo toda a população $i_0$ terá sido substituída por uma nova população infecciosa isto é $i_0-nc=0$, uma vez que a cada instante $t=n$ tem-se que uma fração $c$ da população suscetível sendo infectada, e outra fração $c$ da população infecciosa que se recuperou pois foi infectada em um instante anterior $t-{\tau }_i$. . Reescrevendo a integral com este novo limite

$$i\left(t\right)=\underset{t^{*}-n}{\overset{t^{*}}{\int }}\beta s\left(t\right)i\left(t\right)dt=\beta s^{*}{i}^{*}n=\frac{\beta s^{*}{i}^{*}}{\beta s_0}=i^{*}$$

Uma vez que $s_0=s^{*}$ para estas condições iniciais. Além disto, a expressão para $n$ é análoga a equação encontrada para o ponto de equilíbrio do sistema sem atraso, isto é ${\tau }_i=1/\beta s^{*}$. Toda discussão foi feita para condições iniciais constantes, quando se utiliza o método proposto, a análise torna-se mais complicada, mas ainda é notável que o ponto de equilíbrio depende das condições iniciais.Utilizando o conjunto de parâmetros $\beta =0.4$, ${\tau }_i=5$ e ${\tau }_r=50$, para diferentes valores iniciais é possível obter:

Solução em equilíbrio do sistema de equações diferenciais com atraso para diferentes valores iniciais de infecciosos.

Pode-se notar que aproximadamente:

$i_0=$	${\tau }_{i}c\approx$	$i^{*}\approx$
$0.20$	$0.06$	$0.26$
$0.15$	$0.06$	$0.21$
$0.10$	$0.06$	$0.16$
$0.05$	$0.06$	$0.11$

Lembrando da integral proposta para a condição de equilíbrio condição de equilíbrio:

$$\begin{array}{rl}i\left(t\right)& =c_1+\underset{t^{*}-{\tau }_i}{\overset{t^{*}}{\int }}c\left(t\right)dteq:veinte\end{array}$$

Novamente $c_1\ne 0$, mais especificamente $i_0=c_1$. Este resultado faz sentido, pois de acordo com o procedimento proposto, não há momento em que os infecciosos iniciais $i_0$ se recuperam. Durante o tempo $0\le t<{\tau }_i$ uma fração da população se infecta, e então esta mesma fração começa a se recuperar após $t>{\tau }_i$, mas os infecciosos iniciais, nunca se recuperam, funcionando como fontes permanentes de infecção. Outra forma de comparar, é que olhando o conjunto de equações diferenciais originais o termo responsável pela recuperação dos infecciosos no instante $t$ depende da quantidade total de infecciosos no próprio instante $t$. No conjunto de equações com atraso, é a fração da população da população que foi contagiada em um instante anterior. Ou seja, uma fração da população que se infectou em um instante $t$, é a mesma que se recuperara em um instante posterior $t+{\tau }_i$, porém a quantidade inicial de infecciosos não se contagiou em nenhum momento desta ’linha do tempo', então também não se recuperam em nenhum momento.

Porém ainda que $c_1\ne 0$, quanto mais próximo a população inicial de infecciosos for de $0$, melhor será aproximação, Por isso a análise em torno do ponto de equilíbrio, e o próprio ponto de equilíbrio se torna uma aproximação válida apenas quando é adotado valores próximos de $0$ para $i_0$.

Análise do ponto de equilíbrio

Obtendo a equação transcendental

Aplicando então uma perturbação nos pontos de equilíbrios, isto é: $s\left(t\right)=s^{*}+x\left(t\right)$ e $i$$\left(t\right)=i^{*}+y\left(t\right)$ pode-se obter uma aproximação linear em torno do ponto de equilíbrio (uma observação é que o teste utilizando os limites leva a uma indeterminação, uma vez que não há termo linear). Então analisando o termo que chamamos de contágios:

$$\begin{array}{rl}c\left(t\right)& =\beta s\left(t\right)i\left(t\right)\\ \frac{c\left(t\right)}{\beta }& =(s^{*}+x\left(t\right))(i^{*}+y\left(t\right))\\ & =s^{*}{i}^{*}+i^{*}x\left(t\right)+s^{*}y\left(t\right)+x\left(t\right)y\left(t\right)\end{array}$$

Ignorando os termos não lineares e denotando $f\left(t\right)=f$ e $f(t-\tau )=f_{\tau }$:

$$\begin{array}{rl}\frac{c}{\beta }& =s^{*}{i}^{*}+i^{*}x+s^{*}y\end{array}$$

Tem-se então:

$$\begin{array}{rl}\frac{\dot{x}}{\beta }=& -(s^{*}{i}^{*}+i^{*}x+s^{*}y)+(s^{*}{i}^{*}+i^{*}{x}_{{\tau }_0}+s^{*}{y}_{{\tau }_0})\\ & =-(i^{*}x+s^{*}y)+(i^{*}{x}_{{\tau }_0}+s^{*}{y}_{{\tau }_0})\end{array}$$

E:

$$\begin{array}{rl}\frac{\dot{y}}{\beta }=& -(s^{*}{i}^{*}+i^{*}{x}_{{\tau }_i}+s^{*}{y}_{{\tau }_i})+(s^{*}{i}^{*}+i^{*}x+s^{*}y)\\ & =(i^{*}x+s^{*}y)-(i^{*}{x}_{{\tau }_i}+s^{*}{y}_{{\tau }_i})\end{array}$$

O sistema linearizado em torno do ponto de equilíbrio é então:

$$\begin{array}{rl}\frac{\dot{x}}{\beta }& =-(i^{*}x+s^{*}y)+(i^{*}{x}_{{\tau }_0}+s^{*}{y}_{{\tau }_0})\\ \frac{\dot{y}}{\beta }& =+(i^{*}x+s^{*}y)-(i^{*}{x}_{{\tau }_i}+s^{*}{y}_{{\tau }_i})\end{array}$$Reescrevendo matricialmente:$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\beta \left(\begin{array}{cc}-i^{*}& -s^{*}\\ i^{*}& s^{*}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\beta \left(\begin{array}{cc}{i}^{*}& s^{*}\\ 0& 0\end{array}\right)\left(\begin{array}{c}{x}_{{\tau }_0}\\ y_{{\tau }_0}\end{array}\right)+\beta \left(\begin{array}{cc}0& 0\\ -i^{*}& -s^{*}\end{array}\right)\left(\begin{array}{c}{x}_{{\tau }_i}\\ y_{{\tau }_i}\end{array}\right)$$Ou ainda:$$\dot{\mathit{x}}\left(t\right)=A_0\mathit{x}+A_1\mathit{x}(t-{\tau }_0)+A_2\mathit{x}(t-{\tau }_i)$$

A equação característica para um sistema do tipo:

$$\dot{\mathit{x}}\left(t\right)=A_0\mathit{x}\left(t\right)+\sum _{i=1}^k{A}_i\mathit{x}(t-T_i)$$ É dado por:

$$\det \left[\mathrm{\Delta }\left(\lambda \right)\right]=\det [\lambda \mathbb{𝕀}-A_0-\sum _{i=1}^k{A}_i{e}^{-\lambda T_i}]$$

A função $\mathrm{\Delta }\left(\lambda \right)$ é chamada de quase-polinômio característico. Como de costume, se para para todo $\lambda \in \mathbb{ℂ}$, for $Re\left(\lambda \right)<0$ o sistema é assintoticamente estável.

Detalhes:

• Se não houver atraso, a equação é reduzida a $\det [\lambda \mathbb{𝕀}-A_0]$. Como para achar as raízes deve-se fazer $\det \left[\mathrm{\Delta }\left(\lambda \right)\right]=0$, obtém-se o resultado tradicional.
• Se houver apenas uma equação, e não um sistema, então os termos $A_j$ serão constantes e não uma matrizes. Substituindo $x\left(t\right)=c{e}^{\lambda t}$:

$$\begin{array}{rl}\lambda c{e}^{\lambda t}& =A_{0}c{e}^{\lambda t}+\sum _{i=1}^k{A}_{i}c{e}^{\lambda t}{e}^{-\lambda T_i}\\ \lambda & =A_0+\sum _{i=1}^k{A}_i{e}^{-\lambda T_i}\end{array}$$

Então:

$$\lambda -A_0-\sum _{i=1}^k{A}_i{e}^{-\lambda T_i}=0$$

Obtém-se o quase-polinômio característico proposto pro sistema conforme. Dessa forma, a equação transcendental característica para o sistema que estamos interessados é:

$$\begin{array}{rl}\det \left[\mathrm{\Delta }\left(\lambda \right)\right]=& \det [\lambda \mathbb{𝕀}-A_0-A_1{e}^{-\lambda {\tau }_0}-A_2{e}^{-\lambda {\tau }_i}]\\ =& \det [\left(\begin{array}{cc}\lambda & 0\\ 0& \lambda \end{array}\right)-\beta \left(\begin{array}{cc}-i^{*}& -s^{*}\\ i^{*}& s^{*}\end{array}\right)-\beta \left(\begin{array}{cc}{i}^{*}& s^{*}\\ 0& 0\end{array}\right)e^{-\lambda {\tau }_0}+\beta \left(\begin{array}{cc}0& 0\\ i^{*}& s^{*}\end{array}\right)e^{-\lambda {\tau }_i}]\\ =& \det \left[\left(\begin{array}{cc}\lambda +\beta i^{*}(1-e^{-\lambda {\tau }_0})& \beta s^{*}(1-e^{-\lambda {\tau }_0})\\ \beta i^{*}(e^{-\lambda {\tau }_i}-1)& \lambda +\beta (e^{-\lambda {\tau }_i}-1)s^{*}\end{array}\right)\right]\\ =& [\lambda +\beta i^{*}(1-e^{-\lambda {\tau }_0})][\lambda +\beta (e^{-\lambda {\tau }_i}-1)s^{*}]-i^{*}\beta (e^{-\lambda {\tau }_i}-1)s^{*}\beta (1-e^{-\lambda {\tau }_0})\\ =& {\lambda }^2+\lambda \beta [s^{*}(e^{-\lambda {\tau }_i}-1)-i^{*}(e^{-\lambda {\tau }_0}-1)]+{\beta }^2{s}^{*}{i}^{*}(e^{-\lambda {\tau }_i}-1)(1-e^{-\lambda {\tau }_0})-{\beta }^2{s}^{*}{i}^{*}(e^{-\lambda {\tau }_i}-1)(1-e^{-\lambda {\tau }_0})\\ =& {\lambda }^2+\lambda \beta [s^{*}(e^{-\lambda {\tau }_i}-1)-i^{*}(e^{-\lambda {\tau }_0}-1)]\end{array}$$

Para encontrar as raízes

$$\begin{array}{rl}{\lambda }^2+\lambda \beta [s^{*}(e^{-\lambda {\tau }_i}-1)-i^{*}(e^{-\lambda {\tau }_0}-1)]& =0\end{array}$$

Esta equação é chamada de equação característica transcendental. Uma raiz é $\lambda =0$, e a outra pode ser obtida a partir de:

$$\begin{array}{rl}\lambda +\beta [s^{*}(e^{-\lambda {\tau }_i}-1)-i^{*}(e^{-\lambda {\tau }_0}-1)]& =0\end{array}$$

Como uma raiz é ${\lambda }_1=0$, se ${\lambda }_2>0$ há um equilíbrio instável próximo ao ponto de equilíbrio.

• Observações:
  • Equação transcendental: equação que contém uma função transcendental da variável que deve ser resolvida.
  • Função transcendental: função que não satisfaz uma equação polinomial, em contraste com uma função algébrica, isto é, transcende a álgebra”. Um exemplo é a função exponencial.
  • Função algébrica: função que pode ser definida como a raiz de uma equação polinomial.

Soluções aproximadas são obtidas por métodos gráficos conforme feito anteriormente para soluções reais em Equações diferenciais com atrasos, ou até simplesmente observando o gráfico. Mas neste caso será utilizado algoritmos.

Obtendo as raízes da equação transcendental

Método de Newton

Também conhecido como método de Newton-Raphson, aproxima a função com uma linha. Esta linha atravessa o ponto $(x_i,f\left(x_i\right))$ e tem inclinação igual a derivada da própria função. Matematicamente isto é:

$$\begin{array}{rl}{f}^{\prime }\left(x_i\right)=& \frac{y_{i+1}-f\left(x_i\right)}{x_{i+1}-x_i}\\ x_{i+1}=& x_i+\frac{(y_{i+1}-f\left(x_i\right))}{f^{\prime }\left(x_i\right)}\end{array}$$

Então encontra-se o $x$ que atravessa o eixo e o utiliza como uma nova tentativa. Ou seja o segundo ponto sobre a reta é $(x_{i+1},0)$, então:

$$\begin{array}{rl}{x}_{i+1}=& x_i-\frac{f\left(x_i\right)}{f^{\prime }\left(x_i\right)}\end{array}$$

Para estender o método para o plano complexo, basta substituir variável real $x$ por uma variável complexa $z=x+yi$. A partir disto, o loop é repetido.

Aplicando o método

Então sendo:

$$\begin{array}{rl}f\left(\lambda \right)& =\lambda +\beta [s^{*}(e^{-\lambda {\tau }_i}-1)-i^{*}(e^{-\lambda {\tau }_0}-1)]\\ f^{\prime }\left(\lambda \right)& =1+\beta [{\tau }_0{i}^{*}{e}^{-\lambda {\tau }_0}-s^{*}{\tau }_i{e}^{-\lambda {\tau }_i}]\end{array}$$

Foi utilizado um algoritmo em Python para buscar as raízes complexas. O gráfico original é $R_0\times \frac{{\tau }_r}{{\tau }_i}$. Onde $R_0=\beta {\tau }_i$ é chamado de taxa de reprodução básica. Para facilitar foi escolhido denotar $R_1=\frac{{\tau }_r}{{\tau }_i}$. Com poucas manipulações é possível mostrar que:

$$\begin{array}{rl}{\tau }_i=& \frac{R_0}{\beta }\\ {\tau }_r=& \frac{R_1{R}_0}{\beta }\end{array}$$

Assim para um $\beta$ constante pode-se variar apenas $R_0$ e $R_1$. Deste modo o algoritmo proposto busca a existência de raízes com parte real positiva, isto é, o conjunto de equações apresenta comportamento instável na proximidade do ponto de equilíbrio. Esta busca ocorre para cada par de valores $(R_1,R_0)$ na malha formada por $R_0\times R_1$, onde $1<R_0<10$ e $0\le R_1<3.5$, com intervalos de $\mathrm{\Delta }=0.1$. Especificamente em cada par de valores é executado o algoritmo de Newton-Raphson com o chute inicial partindo de cada ponto possível do plano complexo dentro da região delimitada por $-1\le \mathbb{ℝ}\left(\lambda \right),\mathbb{𝕀}\left(\lambda \right)<1$ , com espaçamento entre os pontos de $\mathrm{\Delta }=0.1$. Isto ocorre enquanto nenhuma raiz com parte real positiva seja encontrada, além disso, é realizada no máximo $N=10^7$ tentativas de aproximação, e é considerado $f\left(x\right)=0$ se $\left|f\left(x\right)\right|<10^{-12}$.

A esquerda a imagem original retirada do artigo. A região preta são soluções sem oscilação e a colorida onde há a presença de oscilações. A direita a imagem obtida através do algoritmo proposto onde a região verde indica que foi encontrado autovalores positivos.

Algoritmos

O conjunto de equações diferenciais sem atraso pode ser solucionado via Mathematica:

tmax = 400; b = 0.4; ti = 5; tr = 5;
sol = NDSolve[{
    s'[t] == -b*s[t]*i[t] + r[t]/tr,
    i'[t] == +b*s[t]*i[t] - i[t]/ti,
    r'[t] == i[t]/ti - r[t]/tr,
    s[0] == 0.5, i[0] == 0.5, r[0] == 0},
   {s, i, r}, {t, 0, 500}];
Plot[{i[t] /. sol, s[t]} /. sol, {t, 0, 500}, PlotRange -> {0., 1.}]

O sistema de equações diferenciais com atraso pode ser pode ser solucionado via Python utilizando o método de Euler conforme proposto :

import matplotlib.pyplot as plt
import numpy as np
#Função para resolver o sistema de equações diferenciais
def sistema(R0,R1):
    #Parâmetro R0=b.ti
    #Parâmetro R1=tr/ti
    #Parâmetros da dinâmica
    b =4
    ti=R0/b
    tr=R1*R0/b
    to=ti+tr
    # Listas para guardar a evolução do sistema
    s=[]
    i=[]
    d=0.001 #Passo para o método de Euler
    #Primeira parte:
    N1=int(ti/d)     #Quantidade de passos
    i.append(1e-16)  #Condição inicial de inectado i0
    s.append(1-i[0]) #Condição inicial de suscetíveis s0=1-i0
    #Resolve o sistema Usando o método de Euler
    for k in range(N1):
        s.append(s[k]+d*(-b*s[k]*i[k]))
        i.append(i[k]+d*(b*s[k]*i[k]))
    #Segunda parte
    N2=int(to/d)
    for k in range(N1,N2):
        s.append(s[k]+d*(-b*s[k]*i[k]))
        i.append(i[k]+d*(b*s[k]*i[k]-b*s[k-int(ti/d)]*i[k-int(ti/d)]))
    #Terceira parte
    N3=int(1800/d)
    for k in range(N2,N3):
        s.append(s[k]+d*(-b*s[k]*i[k]+b*s[k-int(to/d)]*i[k-int(to/d)]))
        i.append(i[k]+d*(+b*s[k]*i[k]-b*s[k-int(ti/d)]*i[k-int(ti/d)]))

O sistema também pode ser resolvido via Mathematica:

b = 0.8; ti = 5; to = 20;
ci1 = NDSolve[{
    s1'[t] == -b*s1[t]*i1[t],
    i1'[t] == +b*s1[t]*i1[t],
    s1[0] == 0.99, i1[0] == 0.01},
   {s1, i1}, {t, 0, ti}];
ci2 = NDSolve[{
    s2'[t] == -b*s2[t]*i2[t],
    i2'[t] == +b*s2[t]*i2[t] - b*s2[t - ti]*i2[t - ti],
    s2[t /; t <= ti] == s1[t] /. ci1, 
    i2[t /; t <= ti] == i1[t] /. ci1},
   {s2, i2}, {t, 0, to}];
sol = NDSolve[{
    s'[t] == -b*s[t]*i[t] + b*s[t - to]*i[t - to],
    i'[t] == +b*s[t]*i[t] - b*s[t - ti]*i[t - ti],
    s[t /; t <= to] == s2[t] /. ci2, i[t /; t <= to] == i2[t] /. ci2},
   {s, i}, {t, 0, 300}];
Plot[{i[t] /. sol, s[t] /. sol}, {t, 0, 300}, PlotRange -> {0., 1.}]

O código abaixo foi escrito em Python e é responsável por buscas raízes complexas com a parte real positiva na equação transcendental discutida anteriormente.

#Código para resolver a equação transcendental
def raizes():
    sol=np.zeros((35,90)) # Matriz para guardar as raízes
    err=1e-12             # Erro admitido   
    sy=0
    for R0 in np.arange(1,10,0.1):      #Percorrer os valores RO=b.
        print(str(100*(sy+1)/90)+"%")
        sx=0
        for R1 in np.arange(0,3.5,0.1): #Percorrer os valores R1=tr/ti
            #Parâmetros da dinâmica
            b =4
            ti=R0/b
            tr=R1*R0/b
            to=ti+tr
            #Pontos de equilíbrio
            io=(b*ti-1)/(b*to)
            so=1/(b*ti)
            #Matriz das raízes para os parâmetros atuais
            r=np.zeros((20, 20))
            #Os chutes iniciais serão pontos dentro da malha [-1,1] em 2D
            m=0
            for a in np.arange(-1,1,0.1):
                n=0
                for y in np.arange(-1,1,0.1):
                    i=y*1j
                    x=a+i  # Ponto inicial
                    N=10000000 # Quantidade máxima de aproximações
                    #Valor da função para o chute inicial
                    f=x+b*(so*(np.power(np.e,-x*ti)-1)-io*(np.power(np.e,-x*to)-1))
                    for k in range(N):
                        if (k==N-2):
                            print("!") #Indicando que chegou no último passo sem achar a raízs
                        if (abs(f)<err):
                            break        # Se chou a raíz, do loop
                        try:
                            f =x+b*(so*(np.power(np.e,-x*ti)-1)-io*(np.power(np.e,-x*to)-1))
                            df=1+b*(to*io*np.power(np.e,-x*to)-ti*so*np.power(np.e,-x*ti))
                            nx=x-f/df        # Próximo chute
                            x=nx
                        except:
                            print(str(R0)+","+str(R1))
                            break
                    if(abs(np.real(x))<err):   # A raíz é zero
                        r[m,n]=0
                    elif (np.real(x)>0):        # Raíz positiva
                        r[m,n]=1
                        break
                    else:
                        r[m,n]=-1               # Raíz negativa
                    n=n+1
                else:        # Quando acabar o loop interno
                    m=m+1
                    continue # Continua
                break        # Se o loop interno foi encerrado antes da hora, encerra o externo
            res=[]        
            for i in range(20):
                res.append(max(r[i]))
            sol[sx,sy]=max(res)                 # Se houve raíz positiva, obtém.
            sx=sx+1
        sy=sy+1
    #Registra
    f = open("raizes.dat", "w")
    for i in sol:
         for j in i:
             f.write(str(j)+"	")
         f.write("\n")
    f.close()
raizes()

A partir dos dados gerados com o código anterior, é possível visualizá-los em R com o seguinte código:

library('tseries')      #Bibliotea pra ler matriz
library('plot.matrix')  #Biblioteca pra plotar matriz

m <- read.matrix( "Raizes.dat") #Os dados são lidos no formato de matriz
m<-t(m)                         # São feitas algumas correções devido à ordem que os dados são gravados
m <-m[90:2, 1:35]
rownames(m) <- (99:11)/10       #Ajusta-se os nomes das linhas e colunas
colnames(m) <- (0:34)/10
#E plota-se
plot(m, col=c('black', 'white', 'green'), breaks=c(-2, -0.25,0.25, 2), xlab='R1', ylab='R0',main='Raizes',border=NA)   

Principais materiais utilizados

1. Oscillations in SIRS model with distributed delays (S. Gonçalves e outros, The European Physical Journal B)
2. Root finding algorithms (Eugeniy E. Mikhailov, Faculdade de William e Mary)
3. Stability and stabilization of time-delay systems ( Gerhard Manfred Schoen ,Instituto Federal de Tecnologia de Zurique)

Anterior: Exemplo: equação diferencial com retardo | Índice: Ecologia | Próximo: Simulação e modelo de campo médio