Giampaolo Mele  (estratto dalla tesi di laurea triennale in Matematica, Università di Pisa 2011)

Le equazioni differenziali con ritardo intervengono in quei modelli matematici in cui la reazione di un sistema sottoposto ad un certo stimolo non è sempre immediata ma può avvenire con un certo ritardo.   A volte tale ritardo è trascurabile, mentre in altri casi trascurarlo può portare ad errori molto grossi.   In questo articolo presentiamo una introduzione al problema generale svolta attraverso esempi ed applicazioni.  Successivamente diamo un cenno alla trattamento numerico di tali problemi toccando gli aspetti computazionali.

1 – Presentazione del problema

Generalmente i problemi del mondo reale in cui intervengono quantità che evolvono nel tempo, o che variano in dipendenza di altre variabili di natura diversa, sono descritti da una o più funzioni di queste variabili.   Le equazioni che caratterizzano il modello spesso legano la funzione incognita, ad esempio $y(t)$, con la sua derivata prima $y’(t)$ (o le sue derivate fino a un certo ordine).   Danno quindi origine a (sistemi di) equazioni differenziali ordinarie.   Nella forma più generale un sistema di equazioni differenziali ordinarie ai valori iniziali assume la forma \[\begin{array}{ll} y' (t) = f (t, y(t)),&\qquad t_0\lt t\le T\\ y(t_0)=y_0 \end{array} \] dove $y(t):[t_0,T]\to\mathbb R^n$ è la funzione incognita che è supposta nota in un punto $t_0$, detto istante iniziale ed $f(t,y):[t_0,T]\times \mathbb R^n\to \mathbb R^n$ è una opportuna funzione di $t$ e $y$ che caratterizza il problema. Questa classe di equazioni ammette un’unica soluzione se $f(t,y)$ è una funzione lipschitziana rispetto al secondo argomento nell’intervallo $[t_0 , T ]$, ovvero esiste un $K>0$ tale che $\|f (t, v) – f (t, w)\|\le K\| v – w\|$ per ogni $v, w \in \mathbb R^n$ e per ogni $t\in[t_0,T]$, dove $\|\cdot\|$ è una norma su $\mathbb R^n$.

        Un semplice esempio è dato dall’equazione \[ \begin{array}{ll} y'(t)=\lambda y(t)&\qquad t\in[t_0,T]\\ y(t_0)=y_0 \end{array} \] che si ottiene dal caso generale con $n=1$ e $f(t,y)=\lambda y$. È evidente che $f$ è lipschitziana in $y$ e che la soluzione è $y(t)=y_0e^{\lambda t}$.   Questa semplice equazione differenziale modella l’evolvere di un sistema la cui variazione nel tempo, espressa dalla derivata prima, è proporzionale alla funzione stessa. Questo accade ad esempio in un capitale finanziario che frutta nel tempo con tasso costante e i cui frutti vengono immediatamente re-investiti.   La stessa equazione descrive un modello semplificato di sviluppo di una popolazione, vista come funzione continua del tempo, in cui non ci sono decessi e la crescita dovuta alle nuove nascite è proporzionale al numero di individui e ogni nuovo nato è immediatamente in grado di riprodursi.

        I problemi in questa classe sono detti problemi ai valori iniziali o anche IVP (Initial Value Problems) Le ipotesi sull’esistenza e l’unicità della soluzione di un IVP si possono indebolire.   Per approfondire questo aspetto è sufficiente consultare qualsiasi libro di analisi, per esempio si veda il capitolo 6 di [4].

        Un esempio più significativo è dato di seguito.

Esempio 1.1.  (Equazione logistica).   Questo è un modello matematico proposto da Verhulst nel 1938 per descrivere la crescita delle popolazioni \[\left\{\begin{array}{ll} \frac{dN}{dt} = rN (t) \left(1 - \frac{N (t)}{K}\right)&\qquad t_0\lt t\le t_f\\ N(t_0)=N_0 \end{array}\right. \]

• $N$ è il numero di individui della specie
• $r$ è il tasso di crescita
• $K$ è la capacità di carico della popolazione, cioè il numero massimo di individui che possono trovarsi in un determinato ambiente.

In questo modello il numero $N(t)$ degli individui di una popolazione al tempo $t$ viene trattato come se fosse una funzione continua e derivabile del tempo, per cui la variazione di una popolazione (accrescimento o diminuzione del numero degli individui) viene espressa dalla derivata di $N(t)$ rispetto al tempo.   L’equazione esprime il fatto che la variazione della popolazione è proporzionale al numero totale degli individui e alla distanza tra il valore corrente della popolazione e il valore massimo raggiungibile. La proporzionalità rispetto al numero degli individui è naturale visto che ogni individuo può riprodursi o estinguersi.   La proporzionalità con la differenza $K-N(t)$ esprime la capacità dell’ambiente di ricevere nuovi individui.   Se le risorse disponibili, ad esempio il cibo e l’acqua, non sono sufficienti ad alimentare la specie è ragionevole supporre che la specie abbia meno capacità di sviluppo e questo accade in misura maggiore quanto più la popolazione è vicina al valore limite.

        È possibile dare una forma esplicita alla soluzione di tale equazione nel modo seguente \[ N (t) = \frac{KN_0e^{rt}}{K+N_0(e^{rt}-1)} \] Nella figura che segue si riporta il grafico della soluzione che rappresenta la crescita del Paramecium caudatum, assieme al grafico relativo all’osservazione del fenomeno reale.

1.2 – Equazioni differenziali con ritardo

I problemi ai valori iniziali sono impostati sull’idea intuitiva che interagire con un sistema provoca in esso una “immediata” variazione, tale idea è sensata dato che se si colpisce una palla questa immediatamente si sposta oppure se si lascia cadere un oggetto pesante questo non fluttuerà nell’aria prima di cadere. Daltronde non è sempre così e ci sono vari esempi in cui la reazione non è istantanea.   Si pensi alla doccia.   Nel fare la doccia, quando si sposta la leva del rubinetto termostatico verso la posizione di acqua calda, la temperatura dell’acqua che esce dal soffione non aumenta immediatamente ma occorrerà aspettare qualche istante affinch´ l’acqua calda percorra il tragitto che separa il rubinetto dal soffione.   Stessa situazione si incontra nel caso in cui l’acqua che esce dalla doccia è troppo calda e vogliamo ridurne la temperatura spostando la leva del rubinetto verso la posizione di acqua fredda.

     In generale un’equazione differenziale con ritardo, denotata con DDE (Delay Differential Equation) è un problema della forma \[ \begin{array}{ll} y'(t)=f(t,y(t),y(t-\tau))&\qquad t_0 \lt t \le T \\ y(t)=\phi(t) &\qquad t \le t_0 \end{array} \] Dove $\tau=\tau(t,y(t))$, che in generale è a sua volta una funzione, è detto ritardo, spesso chiamato anche lag, mentre $\phi$ è detta storia.   Si osserva subito che non è sufficiente conoscere un valore iniziale della soluzione come nei problemi ai valori iniziali, ma occorre conoscere tutta la storia della soluzione antecedente l’istante iniziale.   Infatti si possono costruire esempi di equazioni differenziali con ritardo che hanno soluzioni diverse tra loro ma con lo stesso valore iniziale. In analogia con il caso degli IVP, si ha che la soluzione di questi problemi esiste ed è unica se

• la funzione $f$ è continua e lipschitziana sugli ultimi due argomenti
• la funzione $\tau(t,y(t))$ è non negativa, continua e lipschitziana sull’ultimo argomento
• la funzione $\phi$ è continua e lipschitziana

La dimostrazione di queste proprietà è decisamente più complessa rispetto al caso degli IVP senza ritardo.   Per dettagli sulla dimostrazione si veda per esempio il capitolo 2 di [2]. I problemi si incontrano nel caso in cui $\tau$ si annulla.   Infatti se $\tau > 0$ allora la dimostrazione segue dai metodi risolutivi che si mostreranno nella prossima sezione.

1.3 – Equazione logistica con ritardo

Nel caso dell’equazione logistica data nell’esempio 1.1 si è mostrato un modello per studiare la dinamica delle popolazioni in cui l’ambiente, dotato di un numero limitato di risorse, condiziona la crescita della popolazione.  Tale modello infatti funziona abbastanza bene, ma dà per scontato che se l’ambiente cambia, ad esempio la popolazione cresce e consuma più risorse, istantaneamente gli individui della specie se ne accorgono.   Ma nella realtà non avviene esattamente così.

        È sensato in queste circostanze pensare che l’influenza del cambiamento sulla popolazione avvenga con un certo ritardo, cioè che gli individui si accorgano solo dopo qualche tempo dell’avvenuto cambiamento dell’ambiente.

        Un’altra conferma della necessità di introdurre modelli con ritardo sta nel fatto che il modello non è in accordo con i dati sperimentali. Ad esempio, nel caso dell’equazione logistica, la soluzione è una funzione crescente e limitata da $K$ mentre dalle osservazioni di alcuni fenomeni in biologia si è trovato che in alcune specie a volte $N(t)$ supera $K$ e in alcuni casi $N(t)$ decresce in certi intervalli di tempo. Quindi il modello dell’equazione logistica senza ritardo non è sempre adeguato a descrivere la realtà.

        Pertanto è conveniente modificare il modello in modo tale da imporre che gli individui della specie si accorgano con ritardo che l’ambiente che li ospita è cambiato.   In questo modo si arriva a un modello diverso di equazione logistica, proposto da Hutchinson nel 1948, noto come equazione logistica con ritardo. \[ \begin{array}{ll} \displaystyle \frac{d N}{d t} = r N(t) \left(1 - \frac{N(t-\tau)}{K} \right)& t_0 \lt t \le t_f \\ N(t)=\phi(t) & t \le t_0 \end{array} \] In questa formulazione, mentre rimane istantaneo il cambiamento dovuto alle nascite degli individui (proporzionalità rispetto a $N(t)$), l’influenza dell’ambiente (proporzionalità rispetto a $K-N(t)$) avviene ora con un certo ritardo ottenuto imponendo la proporzionalità rispetto al fattore $K-N(t-\tau)$.

        È interessante vedere come cambia la soluzione dell’equazione logistica nel modello così modificato.   Si supponga per semplicità $\tau$ costante, $\phi(t)=3$ e $K=100$.   Si ha che

• se $\tau=5$ allora il ritardo è trascurabile, infatti la soluzione cresce fino a $K$ e poi si stabilizza su $K$ come se il ritardo non ci fosse;
• se $\tau=10$ Il ritardo non è più trascurabile, ci sono delle oscillazioni iniziali dove $N$ supera $K$ ma dopo un certo tempo la soluzione si stabilizza su $K$;
• se $\tau=20$ si ottiene un caso ben più interessante: si ha che $N$ supera in continuazione $K$, la soluzione è periodica e ci sono dei veri e propri cicli.   Questi casi sono presenti in natura e per molto tempo sono rimasti inspiegati, ad esempio i lemming (roditori artici) seguono questo andamento, infatti in passato per giustificare questi andamenti che si discostano dal modello senza ritardo, si pensava ad un fenomeno di suicidi di massa.
• Se $\tau=27$ si ottiene un caso degenere.   Infatti anche se la soluzione è periodica, dopo il primo picco, per un certo tempo vale $N=0$, quindi è avvenuta l’estinzione della specie.

Grafici al variare di $\tau$

tau = 5	tau=10
tau = 20	tau = 27

I dettagli su questo esempio si trovano in [1] sezione 1.4. Molti altri esempi di DDE si trovano in [2].

2 – Metodi di risoluzione: il metodo dei passi

Per risolvere le DDE il metodo più comodo è quello dei passi che permette di ridurre una DDE ad una successione di IVP.   In questo modo la teoria degli IVP può essere applicata alle DDE. Inoltre i metodi numerici disponibili per risolvere gli IVP, quali i metodi di Runge-Kutta possono essere efficacemente utilizzati per la risoluzione numerica delle DDE. Diamo una breve descrizione di come sia possibile ricondurre una DDE ad una successione di IVP.

2.1 – Ritardo costante

Per maggiore chiarezza consideriamo prima il caso in cui il ritardo $\tau$ è una funzione costante. Si consideri quindi la DDE \[ \begin{array}{ll} y'(t) = f(t,y(t),y(t- \tau)) &\qquad t_0 \lt t \le T \\ y(t)=\phi(t)&\qquad t \le t_0 \end{array} \] Focalizziamo la nostra attenzione sulla soluzione ristretta all’intervallo $[t_0,t_0+\tau]$.   Dato che per $t$ appartenente a questo intervallo si ha che $t-\tau \in [t_0 - \tau, t_0]$, il valore di $y(t-\tau)$ risulta noto per tutti i $t\in[t_0,t_0+\tau]$ coincidendo con la storia passata $\phi(t-\tau)$ della soluzione.   Quindi è possibile trasformare il problema sostituendo $y(t-\tau)$ con $\phi(t-\tau)$ ottenendo il seguente IVP \[ \begin{array}{ll} y_1'(t) = f(t,y_1(t),\phi(t- \tau))&\qquad t_0 \lt t \le t_0+ \tau \\ y_1(t_0)=\phi(t_0) \end{array} \]      Risolto questo primo IVP possiamo procedere in modo analogo e determinare la soluzione nell’intervallo $[t_0 + \tau, t_0 + 2 \tau]$.   Infatti conoscendo ora $y_1(t)$ sull’intervallo $[t_0,t_0+\tau]$, possiamo formulare la DDE sull’intervallo $[t_0+\tau,t_0+2\tau]$ come un nuovo IVP \[ \begin{array}{ll} y_2'(t)=f(t,y_2(t),y_1(t-\tau))&\qquad t_0+\tau \lt t \le t_0 + 2\tau \\ y_2(t_0+\tau)= y_1(t_0+\tau) \end{array} \] Chiaramente la soluzione $y(t)$ in $[t_0,t_0+2\tau]$ sarà l’incollamento di $y_1(t)$ e $y_2(t)$, in generale per avere $y(t)$ nell’intervallo $[t_0,t_0 + N \tau]$ si dovranno risolvere $N$ IVP dati da \[ 1 \le i \le N \hspace{1cm} \begin{array}{ll} y_i'(t)=f(t,y_i(t),y_{i-1}(t-\tau))&\qquad t_0+(i-1)\tau \lt t \le t_0 + i\tau \\ y_i(t_0+(i-1)\tau)= y_{i-1}(t_0+(i-1)\tau) \end{array} \] dove si pone $y_0(t):=\phi(t)$.   Quindi $y(t)$ sarà l’incollamento di $y_1(t), \dots, y_N(t)$, ovvero \[ y(t)=y_i(t) \hspace{1cm} \mbox{se} \hspace{1cm} (i-1) \tau \le t \le i \tau \]
  Questo modo di ricondurre la risoluzione di una DDE alla risoluzione di una successione di IVP è noto come metodo dei passi.

Esempio 2.1.  Si consideri la seguente DDE \[ \begin{array}{ll} y'(t)=-y(t-1) &\qquad 0 \lt t \le T \\ y(t)=1 &\qquad t \le 0 \end{array} \] Usando il metodo dei passi, si determina la soluzione in $[0,1]$ , in questo caso si ha l’IVP \[ \begin{cases} y_1'(t)=-1 \hspace{1cm} 0 \lt t \le 1 \\ y_1(0)=1 \end{cases} \] che ha come soluzione $y_1(t)=1-t$.   Nell’intervallo $[1,2]$ si ha \[ \begin{cases} y_2'(t)=t-1 \hspace{1cm} 1 \lt t \le 2 \\ y_2(1)=0 \end{cases} \] Che ha come soluzione $y_2(t)=\frac{t^2}{2}-t + \frac{1}{2}$.



Procedendo in questo modo si ottiene la soluzione su tutto $[0,T]$.

Osservazione 2.1.  Si prenda in considerazione l’esempio appena fatto per fare un confronto tra IVP e DDE. \[ \begin{cases} y'(t)=-y(t-\tau) \hspace{1cm} 0 \le t \le T \\ y(t)=1 \hspace{3.1cm} t \le 0 \end{cases} \] Se il ritardo $\tau$ è nullo allora il problema è un IVP e di cui si conosce già che la soluzione $x(t)=e^{-t}$.  Una proprietà importante di questa funzione è che tende a zero per $t\to\infty$. Introducendo un ritardo non trascurabile è lecito chiedersi se questa proprietà di decadimento si conserva. Riportiamo di seguito alcuni grafici.

Grafici al variare del ritardo

tau = 0, T=5	tau = 1, T=20
tau = 2, T = 50

Come è possibile vedere dai grafici in figura se il ritardo non è troppo grande la soluzione ancora tende a zero per $t\to\infty$ pur con qualche oscillazione. Ma se il ritardo diventa troppo grande (ad esempio $\tau=2$) la soluzione inizia ad oscillare divergendo.   Intuitivamente questo comportamento fa pensare all’esistenza di diagrammi di biforcazione.

2.2 – Ritardo variabile

Si consideri una DDE con ritardo non costante \[ \begin{array}{ll} y'(t) = f(t,y(t),y(t- \tau (t,y(t))))&\qquad t_0 \lt t \le T \\ y(t)=\phi(t) &\qquad t \le t_0 \end{array} \] Data la suddivisione \[ \Delta= \left \{ t_0, t_1, \dots , t_n , \dots , t_N = T \right \} \] risolvere la DDE equivale a risolvere $N$ problemi ai valori iniziali \[ \begin{array}{ll} y_{n+1}'(t)= f \left(t, y_{n+1}(t), x(t-\tau(t,y_{n+1}(t))) \right) &\qquad t_n \le t \le t_{n+1} \\ y_{n+1}(t_n)=y_n(t_n) \end{array} \] per $n=0, \dots, N-1$ dove si è posto \[ x(s):=\left\{ \begin{array}{ll} \phi(s) &\qquad s \le t_0 \\ y_{i}(s) &\qquad t_{i-1} \le s \le t_i \end{array}\right. \]

2.3 – Criterio per la scelta della suddivisione

Nell’algoritmo appena descritto, affinchè il problema sia ben posto è necessario scegliere nel modo opportuno $\Delta$.   Quello che si vuole è che la variabile $t-\tau$ nel terzo argomento della funzione $f$, coinvolga solo il passato e non il presente.   Ovvero si richiede che, data la suddivisione $\Delta$, valga $t-\tau \lt t_n$ per ogni valore $t \in [t_n,t_{n+1}]$. Tale condizione è assicurata se scegliamo la suddivisione $\Delta$ in modo che \[ t_{n+1} - t_n \lt \min{\tau}. \]      Se $\tau$ è una funzione positiva allora non si hanno problemi dato che questa assume minimo (positivo) in un intervallo compatto.   Se invece $\tau$ si annulla in qualche punto questo approccio non è più praticabile in questa forma che abbiamo descritto ed esistono metodi basati su tecniche più complesse per trattare questo caso più generale.

     D’ora in avanti si considererà $\tau$ positiva e $\Delta$ che rispetta la condizione $t_{n+1} – t_n < \min{\tau}$.

3 – Regolarità della soluzione delle DDE

Teorema 3.1.  Data una DDE \[ \begin{array}{ll} y'(t) = f(t,y(t),y(t- \tau (t,y(t)))&\qquad t_0 \lt t \le t_f \\ y(t)=\phi(t) &\qquad t \le t_0 \end{array} \] Se $f \in C^k[a,b]$ allora $y$ è una funzione $C^k$ a tratti.

        Si osserva che l’incollamento tra la storia e la soluzione può non essere $C^1$, basti considerare l’esempio 2.1, il grafico della soluzione è nella figura 5 ed è chiaro che questa funzione non è liscia in $0$.



Quello che si può provare (basta derivare la soluzione) è che la derivata seconda di $y$ è discontinua in $1$; in generale, provando a fare il calcolo, ci si accorge che $y^{(k+1)}$ è discontinua in $k$.

        Il teorema enunciato può esser in realtà rafforzato, si può infatti dire che: l’incollamento tra la storia e la soluzione non è liscio (nel senso che non è $C^k$) se e solo se la soluzione è $C^k$ a tratti, ma nella pratica tale incollamento non è mai liscio.   La dimostrazione di questo teorema segue dall’algoritmo che si utilizza per localizzare le discontinuità.

        Si vuole sottolineare inoltre che il teorema 3.1 evidenzia che la risoluzione numerica delle DDE è decisamente complessa.   Per approfondire questi aspetti computazionali si richiama il lavoro di tesi [6], dove questi problemi vengono esaminati in dettaglio.

4 – Conclusioni

Sono state introdotte le equazioni differenziali con ritardo e sono stati presentati degli esempi. È stato osservato che la presenza del ritardo può provocare un cambiamento qualitativo della soluzione, in particolare al variare del ritardo alcuni equilibri da stabili diventano instabili e viceversa.   Abbiamo infine dato la descrizione di un metodo per ridurre una DDE ad una successione di IVP. In questo modo la risoluzione numerica di una DDE viene ricondotta alla risoluzione numerica di una successione di IVP condotta, ad esempio, mediante i metodi di Runge-Kutta.

     Quello delle DDE è un campo molto vasto e molti aspetti della teoria sono ancora oggetto di ricerca.

Bibliografia

[1] Thomas Erneux, Applied Delay Differential Equations, Surveys and Tutorials in the Applied Mathematical Sciences, Springer 2009.   ISBN: 0387743715

[2] Alfredo Bellen e Marino Zennaro, Numerical methods for delay differential equations, Numerical Mathematics and Scientific Computation, Oxford University Press, USA 2003.

[3] Uri Ascher e Linda Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, SIAM 1997.

[4] Alan Feldstein e Kenneth Neves, High Order Solutions for state-dependent delay differential equations with nonsmooth solutions, Society for Industrial and Applied Mathematics 1984

[5] Paolo Acquistapace, Appunti di analisi matematica, Università di Pisa.

[6] Mele Giampaolo, Metodi numerici per equazioni differenziali ordinarie.   Tesi di laurea triennale in Matematica, Università di Pisa a.a. 2010-2011.

Tags: analisi numerica, applicazioni, equazioni con ritardo, equazioni differenziali, errori, modelli matematici