4
prodotte: l’interpretazione di dati satellitari, in grado di osservare aree a livello globale, ha
prodotto una notevole letteratura a riguardo (si veda per esempio Okubo (1971, [84]), Klein
and Hua (1990, [58]), Lesieur and Sardouny (1981, [62]), Trathan et al. (1993, [119]),
Deschamps et al. (1981, [26]), Piontkovski et al. (1997, [93]), Gallagher et al. (1996, [36]),
Denman and Gargett (1995, [25]), Pasquero et al. (2005, [89])).
((a()(a)() ()(b)(t
()(c)(
Figura 0.1: immagini satellitari che riportano la concentrazione di clorofilla relative a: (a) costa
occidentale dell’Australia; (b) Mare Arabico; (c) costa della Florida (si noti in questo caso la
5
visualizzazione della concentrazione anche per il lago di Okeechobee, Everglades National
Park) [82].
Approcci alternativi nella modellizzazione della dinamica del plancton oceanico sono stati
quelli di Sandulescu et al. (2007, [100]), dove si elabora un modello che tiene conto anche dei
fenomeni di trasporto (analizzando l’insorgenza di vortici formatisi in presenza di un’isola ed
effettuando uno studio di stabilità mediante coefficienti di Lyapunov), e di Abraham (1998,
[1]), il quale ricostruisce la concentrazione del fitoplancton a partire dalle osservazioni spettrali
di grandezze fisiche superficiali e quella dello zooplancton mediante la formulazione di un
modello accoppiato “capacità portante-fitoplancton-zooplancton”. Altri modelli presenti in
letteratura: Mackas and Boyd (1979, [66]), Tsuda et al. (1971, [120]), Levin and Segel (1976,
[63]), Wroblewski and O’Brien (1976, [124]), Steele and Henderson (1992, [113]), Malchow
(1994, [67]), Kiorboe and Sabatini (1995, [57]), Gower et al. (1980, [42]), Bennett and
Denman (1985, [18]), Smith et al. (1996, [111]), Fasham (1978, [31]), Arístegui et al. (2004,
[11]), Bracco et al. (2000, [19]), Edwards and Beers (2001, [28]), Edwards and Brindley
(1996, [29]).
Figura 0.2: queste immagini relative al Golfo di Tehuantepec (Messico) mostrano la relazione
tra la temperatura superficiale e la concentrazione di clorofilla [82].
La relativa ricchezza di sviluppi di modelli di popolazione del plankton oceanico non è
compensata da un’adeguata presenza di studi sulla stima dei parametri nei modelli basati sulle
equazioni di reazione-diffusione.
Nella presente tesi si è analizzato un modello di reazione-diffusione descritto da un sistema
di equazioni differenziali di tipo parabolico non lineari. Il modello è stato inoltre stimato
applicando vari metodi e algoritmi. Ciò ha comportato l’insorgere di due problematiche
distinte:
• la necessità di metodi numerici appropriati per la risoluzione di equazioni differenziali
alle derivate parziali, spesso non lineari, e le relative implicazioni dell’approssimazione
delle soluzioni (come la presenza di fenomeni caotici);
• l’aumento del numero dei parametri da stimare, in quanto, rispetto ai modelli di tipo
concentrato, in questo caso giocano un ruolo importante i coefficienti dei termini di
diffusione.
6
L’integrazione di equazioni alle derivate parziali rispetto al tempo e alle coordinate spaziali
utilizza principalmente la tecnica delle differenze finite: il lavoro di Garvie (2007, [ ]), base di
partenza per gli sviluppi di questa tesi, costituisce uno degli esempi più importanti e dettagliati
di sviluppo di un modello discreto di reazione e diffusione della dinamica del plancton
oceanico presenti nella letteratura più recente.
Per risolvere questi problemi, nel lavoro di tesi sono stati inizialmente sviluppati degli
algoritmi di stima dei parametri del termine di reazione e di quello di diffusione partendo da un
modello discretizzato descritto da Garvie (2007, [38]) mediante le differenze finite: l’approccio
di stima si basa su uno schema innovativo, ovvero sulla disponibilità dell’informazione per un
determinato punto dell’ecosistema (situazione paragonabile al monitoraggio in situ), e su
eventuali successive iterazioni per raffinare i valori dei parametri e ridurre l’errore quadratico
medio.
La seconda fase della tesi ha riguardato la formulazione e la simulazione di un modello
discretizzato alternativo, basato sul metodo degli elementi finiti, tecnica numerica che presenta
caratteristiche migliori di convergenza rispetto alle differenze finite, e per questo largamente
impiegata non solo in vari settori ingegneristici (meccanica, scienza delle costruzioni,
progettazione aeronautica, elettromagnetismo applicato, ecc.), ma anche nella modellistica di
sistemi distribuiti aventi dominio spaziale di forma complessa. L’elaborazione di un modello
discretizzato con gli elementi finiti, riconducendo il sistema non lineare di equazioni alle
derivate parziali ad un sistema lineare di equazioni differenziali ordinarie, offre l’opportunità
di utilizzare la stima lineare ai minimi quadrati, per ottenere in forma chiusa, e non in modo
euristico come nell’approccio precedente, i valori dei parametri del termine di reazione e del
coefficiente di diffusione.
L’impiego di un metodo numerico alternativo mantiene la finalità di formulare tecniche di
stima idonee ai modelli di reazione-diffusione per la dinamica di popolazione, nell’intento di
fornire una buona base metodologica per sviluppi futuri.
Figura 0.3: esempio di applicazione del monitoraggio satellitare di dati oceanografici:
osservazione di fioritura tossica del cianobatterio Nodularia Spumigena nel Mar Baltico [82].
7
Capitolo 1
Il modello ecologico di reazione e
diffusione: descrizione, discretizzazione
e simulazione
Il modello ecologico trattato nel lavoro di tesi, introdotto da Murray (2003, [81]), è stato
studiato nella sua formulazione discreta da M. Garvie in una sua recente pubblicazione [38],
già citata nell’introduzione: si tratta di un sistema di equazioni alle derivate parziali del
secondo ordine di tipo parabolico, in cui si possono evidenziare, oltre alla derivata temporale,
un termine di reazione e uno di diffusione. Il modello descrive la dinamica spazio-temporale di
due specie in competizione tra loro: indicando con u(x,y,t) la densità di popolazione della preda
(fitoplancton) e con v(x,y,t) la densità di popolazione del predatore (zooplancton) il sistema
diventa
{ (1.1)
dove h(·)
2
C∈ è la risposta funzionale, che rappresenta il tasso di consumo della preda per
predatore come frazione del tasso di consumo massimo p e soddisfa le seguenti condizioni:
i) h(0)=0,
svkuqvhv
t
v
kupvh
w
u
ruu
t
u
−+Δ=
∂
∂
−
⎟
⎠
⎞
⎜
⎝
⎛
−+Δ=
∂
∂
)(
)(1
2
1
δ
δ
8
ii) 1)(lim =
∞→
xh
x
,
iii) h(·) è strettamente crescente in [0,+∞).
La costante k indica la velocità di saturazione del tasso di consumo in funzione della crescita
della preda, q e r rappresentano i tassi di nascita pro-capite massimi del predatore e della preda
rispettivamente, s è il tasso di morte pro-capite del predatore e w è la capacità portante della
preda. La crescita locale della preda è logistica, mentre la risposta funzionale del predatore è di
Holling del II tipo (Holling, 1959 [51]), generalmente la più utilizzata negli studi del settore
(Gentleman et al. (2003, [41]), Jeschke et al. (2002, [56]), Skalski and Gilliam (2001, [109])). I
parametri r, w, p, k, s e q costituiscono i coefficienti del termine di reazione, mentre δ
1
e δ
2
sono i coefficienti della parte di diffusione (tutti strettamente positivi).
Il modello tiene conto dell’invasione delle specie preda da parte dei predatori, ma non degli
effetti stocastici e delle influenze ambientali: le equazioni di reazione-diffusione mostrano
comunque un largo spettro di comportamenti ecologicamente rilevanti. Per una trattazione
completa ed esaustiva dei modelli di reazione-diffusione in modellistica ambientale si veda
Murray (2003, [81]).
1.1 IL MODELLO CONTINUO ADIMENSIONALE
Le equazioni del sistema vengono per convenienza riportate ad una forma non-dimensionale,
riscalando opportunamente sia le variabili che le coordinate spaziali e temporali come segue:
w
u
u =
~
,
⎟
⎠
⎞
⎜
⎝
⎛
=
rw
p
vv
~
, rtt =
~
,
2
1
1
~
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
=
δ
r
xx
ii
.
Lo stesso avviene per i parametri:
a = kw,
r
q
b = ,
r
s
c = ,
1
2
δ
δ
δ = .
Il sistema risultante è dunque
{ (1.2)
con parametri a, b, c e δ strettamente positivi. Tale sistema è definito su un dominio spaziale
limitato Ω ed ha condizioni iniziali appropriate e condizioni al contorno di flusso zero (per
soddisfare l’assunzione che le specie non possano uscire fuori dal dominio stesso).
Per quanto riguarda la risposta funzionale di Holling del II tipo, lo studio si focalizza su due
funzioni a parametri α, β e γ positivi
cvaubvhvvugv
t
v
auvhuuuvufu
t
u
−+Δ=+Δ=
∂
∂
−−+Δ=+Δ=
∂
∂
)(),(
)()1(),(
δδ
9
i) )(
1
)()(
1
auhh =
+
== η
η
η
ηη , con
α
1
=a , β=b , γ=c
ii) )(1)()(
2
auehh =−==
−
ηηη
η
, con γ=a , αβ=b , β=c
dovute rispettivamente a Holling (1965, [52]) e a Ivlev (1961, [54]). Si hanno così due
differenti tipologie di cinetica:
i) v
u
uv
vug
u
uv
uuvuf γ
α
β
α
−
+
=
+
−−= ),(,)1(),
ii) )1(),(),1()1(),(
uu
evvugevuuvuf
γγ
ααβ
−−
−−=−−−=
La regione biologicamente significativa si ha ovviamente per u ≥ 0 e v ≥ 0.
Considerando le nullclines di tale sistema, ovvero le curve per cui f = 0 e g = 0, si ricava,
dall’analisi della stabilità lineare, che si hanno punti di sella in (0,0) e in (1,0) per entrambi i
tipi di dinamica. I punti stazionari (u
*
, v
*
), corrispondenti alla coesistenza tra preda e predatore,
sono invece
i) ))(1(,
****
α
γβ
αγ
+−=
−
= uuvu , con
γ
γβ
αγβ
−
<> ,
ii)
*
1
)1(
,
1
ln
1
**
**
u
e
uu
vu
γ
α
α
γ
−
−
−
=
⎟
⎠
⎞
⎜
⎝
⎛
−
−= , con
⎟
⎠
⎞
⎜
⎝
⎛
−
−>>
α
α
γα
1
ln,1
Affinché il punto stazionario si trovi nel quadrante positivo, occorre che 0 < u
*
< 1. Per
determinate scelte dei parametri, si rileva un ciclo limite stabile intorno al punto stazionario
instabile (u
*
, v
*
), ovvero le densità di popolazione della preda e del predatore mostrano un
andamento ciclico (si veda a proposito Garvie and Trenchea, 2005 [39]). In figura 1.1 si ha una
visione del piano delle fasi per il modello in questione.
10
Figura 1.1: piano delle fasi per la cinetica di tipo (ii), con α = 1.5, β = 1.0, γ = 5.0 .
Le risposte funzionali di Holling e Ivlev sono state largamente usate nella modellistica
ecologica. Garvie e Trenchea [39] mostrano nel loro lavoro che i sistemi preda-predatore sono
ben posti, ovvero caratterizzati da un’unica soluzione per ogni istante di tempo, dipendente dai
dati in maniera continua. Inoltre, partendo da condizioni iniziali non negative, le soluzioni
rimangono non negative, vincolo necessario perché il modello sia biologicamente attendibile.
Nella letteratura scientifica vi sono numerosi esempi di sistemi di equazioni differenziali
ordinarie (ODE) spazialmente omogenei (Freedman, 1980 [33]; May, 1974 [69]; Murray, 2003
[81]), e in tal caso la cinetica di tipo (i) è conosciuta come modello di Rosenzweig-MacArthur
(Rosenzweig and MacArthur, 1963 [98]). Tuttavia, gli studi di sistemi di reazione-diffusione
che tengano conto della dinamica spazio-temporale sono più rari.
Un importante lavoro è quello di Medvinsky et al. (2002, [70]), dove si sviluppa un modello
di reazione-diffusione per descrivere la dinamica del plancton marino: dalle simulazioni del
sistema si nota la formazione di spirali per differenti condizioni iniziali, con sviluppo di forme
irregolari che preludono a comportamenti caotici.
Esempi simili sono: Malchow and Petrovskii (2002, [68]), Petrovskii and Malchow (2001,
[91]), Petrovskii and Malchow (1999, [90]), Petrovskii and Malchow (2002, [92]), Sherratt et
al. (1997, [106]), Sherratt et al. (2002, [107]), Sherratt et al. (1995, [108]), Alonso et al. (2002,
[8]), Gurney et al. (1998, [45]), Rai and Jayaraman (2003, [95]), Savill and Hogeweg (1999,
[102]).
Più numerosi sono invece gli studi sulla stazionarietà e sulle forme spaziali prodotte da
sistemi preda-predatore: instabilità guidata dalla diffusione (Turing, 1952 [121]), stima dei
coefficienti di diffusione (Segel and Jackson, 1972 [105]; Murray, 1993 [80]; Neurbert et al.,
2002 [83]), caos indotto dalla diffusione (Pascual, 1993 [88]; Rai and Jayaraman, 2003 [95]),
studio della turbolenza nella diffusività (Medvinsky et al., 2002[70]; Petrovskii and Malchow,
2002 [92]).