Stabilità delle traiettorie, modelli di popolazione e criterio di Routh

L'idea

Spesso non serve sapere come evolve un sistema istante per istante: basta sapere se resta fermo, se torna dove era o se si allontana. Per i sistemi lineari la risposta sta nel segno della parte reale degli autovalori, e il criterio di Routh la dà senza calcolarli.

Stabilità di una traiettoria

Con un ingresso che varia nel tempo non c'è più un punto fermo: c'è una traiettoria. Si fissa un ingresso $u_e(t)$ e una traiettoria di riferimento $x_e(t)$ che lo segue:

$$\dot x_e(t)=f\big(x_e(t),u_e(t)\big)$$

La traiettoria parte da $x_{e,0}=x_e(0)$; l'evoluzione $x(t)$, con lo stesso ingresso, parte da $x_0$. La domanda è la stessa degli equilibri: se si parte un po' fuori, si resta vicini alla traiettoria? La traiettoria è stabile se per ogni $\varepsilon>0$ esiste $\delta>0$ tale che

$$\begin{aligned} &\|x_0-x_{e,0}\|<\delta\ \Longrightarrow\\ &\|x(t)-x_e(t)\|<\varepsilon\quad\text{per ogni } t\ge0 \end{aligned}$$

È asintoticamente stabile se in più $\|x(t)-x_e(t)\|\to0$, ed è instabile se non è stabile. L'aggettivo è lo stesso degli equilibri, cambia il sostantivo.

Caso lineare. Con lo stesso ingresso, le due evoluzioni hanno lo stesso integrale di convoluzione, che nella differenza si cancella:

$$\begin{aligned} x(t)&=\Phi(t)x_0+\int_0^t H(t-\tau)u_e(\tau)\,d\tau\\ x(t)-x_e(t)&=\Phi(t)\big(x_0-x_{e,0}\big) \end{aligned}$$

Quindi le condizioni sono quelle degli equilibri: $\Phi(t)$ limitata per la stabilità, $\Phi(t)\to0$ per quella asintotica. Un sistema lineare stabile ha equilibri e traiettorie stabili, e lo stesso vale per l'asintotica.

Due evoluzioni con lo stesso ingresso

$t$

• riferimento $x_e(t)$, da $0$
• $x(t)$, da $0,5$

$a$

−0,3

Per un sistema non lineare il risultato non si estende con l'approssimazione lineare: questa vale vicino a un punto, mentre una traiettoria si allontana da ogni punto. Tutto ciò riguarda lo stato, cioè la stabilità interna.

Diffusione di un'influenza

$S$ sono i sani, $I$ gli infetti, $R$ i guariti. Un sano si ammala incontrando un infetto, con velocità proporzionale a $SI$ secondo $\beta$; entrano $N$ nuovi individui per unità di tempo; ogni classe muore con tasso $m$; gli infetti guariscono con tasso $\gamma$:

$$\begin{aligned} \dot S&=-\beta SI+N-mS\\ \dot I&=\beta SI-(\gamma+m)I\\ \dot R&=\gamma I-mR \end{aligned}$$

Equilibri. La seconda equazione è $I(\beta S-\gamma-m)=0$. Con $I=0$ anche $R=0$ e $S=N/m$: è la popolazione senza infezione,

$$x_{e1}=\begin{pmatrix}N/m\\0\\0\end{pmatrix}$$

Con $S=(\gamma+m)/\beta$ si trova un secondo equilibrio, con infetti presenti:

$$x_{e2}=\begin{pmatrix}\dfrac{\gamma+m}{\beta}\\[1ex]\dfrac{N}{\gamma+m}-\dfrac{m}{\beta}\\[1ex]\dfrac{\gamma}{m}\Big(\dfrac{N}{\gamma+m}-\dfrac{m}{\beta}\Big)\end{pmatrix}$$

Ha senso solo se la sua $I$ è positiva.

Arriva un infetto. Partire vicino a $x_{e1}$ vuol dire aggiungere qualche infetto. La jacobiana in $x_{e1}$ è

$$A=\begin{pmatrix}-m&-\beta N/m&0\\0&\beta N/m-(\gamma+m)&0\\0&\gamma&-m\end{pmatrix}$$

con autovalori $\lambda_1=-m$, $\lambda_2=-m$ e $\lambda_3=\beta N/m-(\gamma+m)$. Decide $\lambda_3$: serve $\beta N/m<\gamma+m$, cioè un tasso di contagio minore di quello con cui gli infetti guariscono o muoiono. L'equilibrio senza infezione è quindi localmente asintoticamente stabile se

$$R_0=\frac{\beta N}{m(\gamma+m)}<1$$

ed è instabile con $R_0>1$. $R_0$ è l'indice di riproduzione. Con $R_0=1$ l'approssimazione lineare non decide.

Gli infetti vicino all'equilibrio senza infezione

$t$

• $I(t)/I(0)$, con $\gamma +m=0,5$

$R_0$

0,7

Il risultato dice su che cosa agire: ridurre i contatti riduce $\beta$, e quindi $R_0$.

Prede e predatori

$x_1$ sono le prede (per esempio pesci), che crescono da sole con tasso $K$; $x_2$ i predatori (squali), che senza prede diminuiscono con tasso $h$. Gli incontri tolgono prede (coefficiente $a$) e nutrono predatori (coefficiente $b$). Il modello di Volterra è

$$\begin{aligned} \dot x_1&=Kx_1-ax_1x_2\\ \dot x_2&=-hx_2+bx_1x_2 \end{aligned}$$

Da $x_1(K-ax_2)=0$ e $x_2(-h+bx_1)=0$ gli equilibri sono l'origine, senza nessuna delle due specie, e la coesistenza in $(h/b,\ K/a)$.

Un trattamento che uccide entrambe le specie in proporzione al loro numero, con tassi $c$ e $d$, aggiunge $-cx_1$ e $-dx_2$:

$$\begin{aligned} \dot x_1&=Kx_1-ax_1x_2-cx_1\\ \dot x_2&=-hx_2+bx_1x_2-dx_2 \end{aligned}$$

L'equilibrio di coesistenza diventa

$$\left(\frac{h+d}{b},\ \frac{K-c}{a}\right)$$

Per ogni $c$ e $d$ positivi (con $c$ minore di $K$) all'equilibrio ci sono più prede e meno predatori. Se le prede sono un parassita e i predatori il suo antagonista, il trattamento fa crescere il parassita.

L'equilibrio si sposta

primadopo$x_1$, prede$x_2$, predatori

$c$

0,3

$d$

0,4

La figura usa $h=b=K=a=1$.

Condizione necessaria sui coefficienti

Per la stabilità asintotica tutte le radici del polinomio caratteristico

$$p(\lambda)=a_n\lambda^n+a_{n-1}\lambda^{n-1}+\dots+a_0$$

devono avere parte reale negativa. Allora ogni fattore è $\lambda+|\lambda_i|$, o $(\lambda+|\alpha|)^2+\omega^2$ per una coppia complessa, con coefficienti positivi. Moltiplicando, ogni coefficiente è una somma di prodotti positivi; $a_n$ davanti cambia, se è negativo, tutti i segni insieme, come in $-(\lambda+1)(\lambda+2)=-\lambda^2-3\lambda-2$:

se tutte le radici hanno parte reale negativa, i coefficienti sono tutti non nulli e dello stesso segno.

La condizione è solo necessaria. $\lambda^7+1$ la viola, quindi non ha tutte le radici a parte reale negativa. Se invece vale, si passa alla tabella.

La tabella di Routh

Le righe sono $n,n-1,\dots,0$. Le prime due hanno i coefficienti alterni, a partire da $a_n$ e da $a_{n-1}$. Ogni riga dopo si calcola dalle due precedenti:

$$b_1=\frac{\begin{vmatrix}a_n&a_{n-2}\\a_{n-1}&a_{n-3}\end{vmatrix}}{-a_{n-1}},\qquad b_2=\frac{\begin{vmatrix}a_n&a_{n-4}\\a_{n-1}&a_{n-5}\end{vmatrix}}{-a_{n-1}}$$

La prima colonna del determinante resta fissa, la seconda scorre a destra. Un elemento che manca vale zero nel calcolo ma non si scrive. Le righe $1$ e $0$ hanno sempre un solo elemento.

Lettura. Nella prima colonna ogni variazione di segno fra elementi consecutivi è una radice a parte reale positiva, ogni permanenza una a parte reale negativa. Tutte permanenze: sistema asintoticamente stabile.

Per $\lambda^4+\lambda^3+2\lambda^2+\lambda+2$:

$$\begin{array}{c|ccc} 4&1&2&2\\ 3&1&1&\\ 2&1&2&\\ 1&-1&&\\ 0&2&& \end{array}$$

La prima colonna $1,1,1,-1,2$ ha due variazioni: due radici a parte reale positiva, sistema instabile.

Uno zero nella prima colonna

Se il primo elemento di una riga è zero, la riga dopo non si calcola. Il sospetto è una radice sull'asse immaginario, che le variazioni non sanno contare; in ogni caso non c'è stabilità asintotica.

• Riga non tutta nulla. Se i primi $j$ elementi sono nulli, si somma alla riga la stessa riga spostata a sinistra di $j$ posti e moltiplicata per $(-1)^j$. Per esempio $(0,0,1,2)$ con $j=2$ diventa $(1,2,1,2)$. Le variazioni contano ancora le radici a parte reale positiva; con sole permanenze le radici hanno parte reale $\le0$.
• Riga tutta nulla. Se si annulla la riga $i$, che è sempre dispari, con gli elementi della riga $i+1$ si scrive il polinomio ausiliario $p_2(\lambda)$, con le sole potenze pari: per esempio, dalla riga $4$ con $c_1,c_2,c_3$ viene $p_2=c_1\lambda^4+c_2\lambda^2+c_3$. Allora $p(\lambda)=p_1(\lambda)p_2(\lambda)$, le radici di $p_1$ seguono la tabella fino alla riga $i+1$ e quelle di $p_2$ sono simmetriche rispetto all'origine. Si sostituisce la riga nulla con i coefficienti di $p_2'(\lambda)$ e si prosegue; con sole permanenze le radici hanno parte reale $\le0$.

Per $\lambda^4+\lambda^3+2\lambda^2+\lambda+1$ la riga $1$ si annulla. La riga $2$ è $(1,1)$, quindi $p_2=\lambda^2+1$ e $p_2'=2\lambda$:

$$\begin{array}{c|ccc} 4&1&2&1\\ 3&1&1&\\ 2&1&1&\\ 1&2&&\\ 0&1&& \end{array}$$

Tutte permanenze, ma la tabella è stata modificata: le radici hanno parte reale $\le0$ e $p_2$ dà $\pm j$. Il sistema è stabile, non asintoticamente.

Le radici di ${\lambda }^4+{\lambda }^3+2{\lambda }^2+\lambda +1$

$j$$-j$ReIm

Un parametro nella tabella

Per quali $k$ ha tutte le radici a parte reale negativa questo polinomio?

$$p(\lambda)=\lambda^4+(k+1)\lambda^3+k\lambda^2+\lambda+1$$

La condizione necessaria chiede $k>0$. La tabella è

$$\begin{array}{c|ccc} 4&1&k&1\\ 3&k+1&1&\\ 2&\dfrac{k^2+k-1}{k+1}&1&\\ 1&-\dfrac{k+2}{k^2+k-1}&&\\ 0&1&& \end{array}$$

Con $k>0$ il terzo elemento è positivo solo se $k^2+k-1>0$, il quarto solo se $k^2+k-1<0$: non possono valere insieme. Nessun $k$ rende il polinomio asintoticamente stabile.

Se un primo elemento si annulla per qualche valore del parametro, la tabella non vale per quel valore: lo si sostituisce nelle righe già calcolate e si prosegue con i numeri. Per $k=0$ il polinomio è $\lambda^4+\lambda^3+\lambda+1$ e la prima colonna è $1,1,-1,2,1$: due radici a parte reale positiva.