Modelowanie matematyczne - kilka oczywistych faktów,
które warto sobie uświadomić
Co to jest Model:
" reprezentacja badanego obiektu w postaci innej, niż ta, w której
występuje on w rzeczywistości.
" W nauce model jest rozumiany jako uproszczona, celowo,
reprezentacja rzeczywistości.
Typy modeli:
" modele o podobieństwie geometrycznym (mapy, makiety)
" modele o podobieństwie kinematycznym
" modele o podobieństwie dynamicznym (makiety stosowane w
tunelach aerodynamicznych)
" modele tworzone przez analogie (hydrauliczno elektryczny)
" modele matematyczne
Model matematyczny to skończony zbiór symboli i relacji
matematycznych oraz ścisłych zasad operowania nimi, przy czym
zawarte w modelu symbole i relacje majÄ… interpretacjÄ™ odnoszÄ…cÄ… siÄ™ do
konkretnych elementów modelowanego wycinka rzeczywistości.
Zbiór symboli i relacji matematycznych - to twór abstrakcyjny;
czynnikiem przekształcającym go w model matematyczny jest
fizyczna interpretacja.
1
Jak tworzymy model?
Problem
b d
Cele modelowania
Kategoria modelu Dane
Hipotezy
Struktura modelu
Teorie
Identyfikacja
Prawa
Wiedza
Empiry
Algorytmy
czna
Obliczenia
Prognozy
Weryfikacja
Zweryfikowa
ny
model
2
Przykład: wzrost populacji bakterii
Niech:
N(t) - ilość bakterii w chwili t
r - prędkość reprodukcji
W takim razie, pomijając wymieranie, możemy przypuszczać, że:
N(t + "t) H" N(t) + rN(t)"t
N(t + "t)- N(t)
= rN(t)
"t
Zakładamy teraz, że rozmiar populacji N(t) jest wielkością ciągłą.
Ma to sens przy następujących założeniach:
1. Liczebność populacji jest duża tak, że dodanie lub odjęcie kilku
osobników nic nie zmienia.
2. Rozmnażanie się pojedynczych osobników nie ma wpływu na
rozmnażanie się innych osobników.(nie ma skorelowanych
gwałtownych zmian liczebności)
Wtedy przechodząc w granicy "t 0 powyższe równanie różnicowe
można zastąpić równaniem różniczkowym:
dN
= rN
dt
Równanie to jest znane jako prawo Malthusa (1798). Zastosował je on
do opisu populacji ludzi na Ziemi i jego wnioski były nieco zatrważające
...
Równanie to możemy scałkować przez rozdzielenie zmiennych
dN
= rdt
N
tt
dN
= rds
+"+"
00
N
t
lnN = rt
0
Czyli mamy wzrost exponencjalny !
ln N(t)- ln N(0) = rt
ëÅ‚ öÅ‚
N(t)÷Å‚
ìÅ‚
lnìÅ‚ ÷Å‚ = rt
N0
íÅ‚ Å‚Å‚
N(t) = N0ert
3
Dyskretne modele jednej populacji
Opis dyskretny jest dobry dla populacji dla których nie ma zachodzenia
pokoleń na siebie.
Będziemy analizować modele opisywane równaniami różnicowymi typu:
Nt+1 = Nt F(Nt ) = f (Nt )
Sztuka modelowania polega na umiejętnym dobraniu f (N) tak aby
równanie powyższe dobrze oddawało obserwowane fakty.
Uwaga: nie ma prostego związku między równaniami różnicowymi i
różniczkowymi:
" Najprostszy przykład:
Nt+1 = rNt Ò! Nt = rt N0
czyli mamy wzrost wykładniczy (a jaki był wzrost dla analogicznego
równania różniczkowego?)
" Jak uwzględnić wpływ pojemności środowiska?
Np.: Nt+1 = rNs gdzie Ns = Nt1-b b takie, że Ns < Nt
Interpretacja: Ns - frakcja populacji przeżywająca do rozrodu
" A co się stanie jeśli wez
miemy dyskretnÄ… kalkÄ™ modelu
logistycznego?
Nt
Nt+1 = rNt ëÅ‚1- öÅ‚ r,K > 0
ìÅ‚ ÷Å‚
K
íÅ‚ Å‚Å‚
Np. dla Nt > K Nt+1 < 0
Można to poprawić tak:
Nt
îÅ‚ Å‚Å‚
öÅ‚
Nt+1 = Nt expïÅ‚rëÅ‚1- ÷łśł
r, K > 0
ìÅ‚
K
íÅ‚ Å‚Å‚
ðÅ‚ ûÅ‚
rNt
öÅ‚
expëÅ‚- ÷Å‚ - czynnik Å›miertelnoÅ›ci zwiÄ…zany z pojemnoÅ›ciÄ… Å›rodowiska
ìÅ‚
K
íÅ‚ Å‚Å‚
4
Pajęczynki - graficzna metoda badania równań
różnicowych
Stany stacjonarne
*
N = 0
* * * *
N = f (N )= N F(N ) Ò!
*
F(N )= 1
Å‚atwo znalez
ć graficznie:
1
0.9
0.8
Nt+1=Nt
0.7
0.6
Nt+1 0.5
Nt+1=f(Nt)
0.4
0.3
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1
Nt
Zachowanie w pobliżu punktu stacjonarnego zależy od sposobu
przecięcia się f (N)i Nt+1 = Nt
Ponieważ prosta Nt+1 = Nt ma stały kąt nachylenia w układzie wsp.
Wystarczy rozważyć jedynie lokalne nachylenie krzywej Nt+1 = f (Nt )
*
(czyli pochodnÄ… f '(N ))
*
f '(N )
1
-1 0
5
*
f '(N )
1
-1 0
Nie stabilny Stabilny  Stabilny  Nie stabilny
zbiega zbiega
oscylujÄ…c mono-
tonicznie
Bifurkacja  modele są na ogół zależne od pewnych parametrów.
Stabilność stanów stacjonarnych może zależeć od wartości parametrów.
Jeśli przy małej zmianie jakiegoś parametru r w sąsiedztwie wartości r0
następuje jakościowa zmiana zachowania modelu (stan stacjonarny
zyskuje lub traci stabilność, pojawiają się nowe stany stacjonarne itp.) to
mówimy że wartość r0 jest wartością bifurkacji.
Zadanie 1:
Dla modelu:
Nt
Å‚Å‚
öÅ‚
Nt +1 = Nt îÅ‚1+ rëÅ‚1-
ìÅ‚ ÷łśł
ïÅ‚
K
íÅ‚ Å‚Å‚
ðÅ‚ ûÅ‚
" Przejść do zmiennych bezwymiarowych
" określić nieujemne stany stacjonarne
" przedyskutować ich liniową stabilność
" znalez
ć minimalną i maksymalną liczebność populacji
" znalez
ć wartość pierwszej bifurkacji
Zadanie 2:
" Skonstruować pajęczynkę dla modelu:
(1+ r)Nt
Nt +1 =
1+ rNt
" Przedyskutować jakościowo globalne zachowanie rozwiązań.
" Podać minimalną i maksymalną wartość Nt .
6
Równanie logistyczne  chaos deterministyczny
Przeanalizujmy zachowanie siÄ™ odwzorowania
xt +1 = rxt(1- xt ) (*)
możemy interpretować je jako bezwymiarową wersję równania na
liczebność populacji z uwzględnieniem pojemności środowiska.
Równanie to ( jak pewnie wiecie ) ma ciekawe zachowanie zależne od
wartości parametru r
1. Znajdz
my jego stany stacjonarne i zbadajmy jak zmienia siÄ™ ich
stabilność w zależności od wartości parametru r
2. Narysujmy pajęczynkę dla tego odwzorowania dla r = 2i dla r = 3.1
3. Znajdz
my stany stacjonarne dla odwzorowania
xt +2 = rxt+1(1 - xt +1) = r(rxt (1 - xt ))(1 - rxt (1 - xt )) (**)
Musimy w tym celu rozwiązać równanie
x = r(rx(1- x))(1 - rx(1- x))
zauważmy, że stany stacjonarne poprzedniego odwzorowania (*) są też
stanami stacjonarnymi odwzorowania (**), a także wyższych iteracji
równania (*). Zatem wielomian:
p(x) = r(rxt (1- xt ))(1 - rxt (1 - xt ))- x możemy podzielić przez
1
ëÅ‚1 öÅ‚
x - -
ìÅ‚ ÷Å‚
r
íÅ‚ Å‚Å‚
otrzymujemy sfaktoryzowany wielomian:
1 1 1 1
öÅ‚ëÅ‚ öÅ‚
ëÅ‚1 öÅ‚÷Å‚ìÅ‚ x2-xëÅ‚1 + öÅ‚ ëÅ‚ öÅ‚÷Å‚
p(x) = -r3x ëÅ‚ x + - + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
2
r r r r
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
íÅ‚ Å‚Å‚íÅ‚ Å‚Å‚
Ostatecznie punktami stacjonarnymi naszego odwzorowania (**) sÄ…
0
Å„Å‚
1
ôÅ‚
1-
ôÅ‚
r
ôÅ‚
1
x =
òÅ‚
(r +1+ (r - 3)(r +1))
ôÅ‚2r
ôÅ‚
1
(r (r
ôÅ‚2r +1- - 3)(r +1))
ół
dwa ostatnie pierwiastki sÄ… rzeczywiste tylko dla r < -1 (" r > 3
Tak więc dla dodatnich r te dwa stany stacjonarne istnieją gdy r > 3,
1
wtedy też stan stacjonarny 1 - traci stabilność. Można pokazać
r
7
1
rachunkiem wprost, że punkty (r +1 ą (r - 3)(r +1)) są stabilne dla
2r
3 < r < r2 H" 3.3.
W analogiczny sposób można znalez
ć i zbadać stabilność punktów
stacjonarnych kolejnych iteracji równania (*) -> orbity o okresach 2n .
Nowy jakościowo efekt pojawia się w momencie gdy r H" 3.83. Pojawia
się wówczas orbita o okresie 3. Udowodnione jest twierdzenie, że jeśli
dla pewnego rc istnieje rozwiązanie dla okresów nieparzystych (co
najmniej 3) to powyżej rc istnieją rozwiązania aperiodyczne.
4. Zbadajmy graficznie co dzieje siÄ™ z kolejnymi iteracjami (*)
8
Układy dwóch równań różnicowych
Rozszerzymy teraz naszą metodologię na układy dwóch równań
różnicowych. Układ dany jest równaniem :
xn+1 = f (xn, yn)
f, g  funkcje nieliniowe
yn+1 = g(xn, yn)
W punktach stacjonarnych mamy:
x* = f (x*, y*)
y* = g(x*, y*)
Aby wnioskować o stabilności stanu stacjonarnego zbadamy co dzieje się
z małymi zaburzeniami (x* + x', y* + y').
Rozwijamy funkcje f i g w szereg Taylora wokół punktu
stacjonarnego:
"f "f
f (x* + x', y* + y')= f (x*, y*)+ x' + y' + K
"x "y
x*, y*
x*, y*
"g "g
g(x* + x', y* + y')= g(x*, y*)+ x' + y' + K
"x "y
x*, y*
x*, y*
i zostawiamy tylko wyrazy liniowe
Oznaczmy:
"f "f
a11 = , a12 =
"x "y
x* , y*
x* , y*
"g "g
a21 = , a22 =
"x "y
x* , y*
x* , y*
wówczas nasze równania w przybliżeniu liniowym wyglądają
następująco:
x'n+1 = a11x'n +a12 y'n
y'n+1 = a21x'n +a22 y'n
w notacji wektorowej
a11 a12 x'n
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
x'n+1 = Ax'n gdzie A = ìÅ‚
ìÅ‚a a22 ÷Å‚ , x'n = ìÅ‚ y'n ÷Å‚
÷Å‚ ìÅ‚ ÷Å‚
íÅ‚ 21 Å‚Å‚ íÅ‚ Å‚Å‚
Teraz musimy znalez
ć pierwiastki równania charakterystycznego:
det(A - I)= 0
2 - tr A + det A = 0
9
Aby punkt stacjonarny był stabilny potrzeba aby oba pierwiastki co do
modułu były mniejsze niż 1.
Jakie wynikają stąd warunki na ślad i wyznacznik?
tr A Ä… tr A2 - 4det A
Pierwiastki tego równania to: 1,2 =
2
Dla rzeczywistych pierwiastków ( tr A2 > 4 det A ):
Wierzchołek paraboli musi być pomiędzy -1 a 1
tr A
-1 < <1
2
2 > tr A
zaś odległość pomiędzy wierzchołkiem paraboli a miejscem zerowym
musi być mniejsza niż odległość wierzchołka od ą1
tr A tr A2 - 4det A
1 - >
2 2
tr A2 tr A2 - 4det A
1 - tr A + >
4 4
1 + det A > tr A
Å‚Ä…czÄ…c oba warunki mamy
2 > 1 + det A > tr A
0.8
0.6
0.4
0.2
1/2 (trA2-4detA)
"
<>
0
trA/2
-0.2
-0.4
-1.5 -1 -0.5 0 0.5 1 1.5
10
Warunek stabilności dla układów wyższego rzędu
Dla układów wyższych rzędów metoda postępowania teoretycznie jest
taka sama - linearyzacja układu w pobliżu punktu stacjonarnego, a
następnie zbadanie wartości własnych jakobianu. Problem praktyczny
polega na tym, że dla wyższych stopni równania charakterystycznego
trudno jest znalez
ć analityczną postać pierwiastków. Pomocne jest tu
kryterium Jury1 (1971) .
Rozważmy wielomian:
P() = n + a1n-1 + a2n-2 +K+ an-1 + an
Potrzebne nam będą kombinacje parametrów konstruowane w
następujący sposób:
2 2 2 2
bn = 1- an cn = bn - b12 dn = cn - c2
bn-1 = a1 - anan-1 cn-1 = bnbn-1 - b1b2 dn-1 = cncn-1 - c2c3
M M M
bn-k = ak - anan-k cn-k = bnbn-k - b1bk +1 dn-k = cncn-k - c2ck +2
M M M
M M d3 = cnc3 - c2cn-1
M c2 = bnb2 - b1bn-1
b1 = an-1 - ana1
aż do momentu gdy zostaną nam trzy wielkości np.:
2 2
qn = pn - pn-3
qn-1 = pn pn-1 - pn-3 pn-2
qn-2 = pn pn-2 - pn-3 pn-1
1
Jury, E.I. (1971) The inners approach to some problems of. IEEE Trans. Automatic Contr. AC-
12,233-240 system theory
11
Na tak zdefiniowanych parametrach test Jury wygląda następująco:
Warunki konieczne i wystarczające na to aby wielomian P() miał
pierwiastki  < 1 są następujące:
P(1)=1+ a1 + a2 +K+ an > 0
n n n n-1
(-1) P(-1)= (-1) [(-1) + a1(-1) +L + an-1(-1)+ an]> 0
an < 1
bn < b1
cn < c2
M
qn < qn-2
Przykład zastosowania testu Jury:
Mamy wielomian:
P(x) = x4 + x3 + x2 + x +1
n = 4
1. P(1) = 4 > 0
2. (-1)4 P(-1) = 1(1-1+1-1+1) = 1 > 0
3. an = 1 nie spelniony warunek
Ponieważ trzeci warunek nie jest spełniony więc ten wielomian ma
przynajmniej jeden pierwiastek nie mniejszy niż 1.
12
Przykład
Prostym przykładem modelu z układem dwóch równań różnicowych jest
model opisujący pasożyty i ich nosicieli w świecie owadów.
Wspólne założenia dla tego typu modeli są następujące:
1. Nosiciele zakażeni pasożytami wydadzą następne pokolenie
pasożytów
2. Nosiciele nie zakażeni wydadzą kolejne pokolenie własnego
gatunku
3. Frakcja nosicieli zakażonych zależy od częstości spotkań obu
gatunków, najczęściej zależy od gęstości obu gatunków.
Do opisu modelu wprowadzimy następujące zmienne:
Nt - gęstość nosicieli w pokoleniu t
Pt - gęstość pasożytów w pokoleniu t
f = f (Nt, Pt )- frakcja nosicieli, która nie została zakażona
 - prędkość reprodukcji nosicieli
c - średnia ilość jaj pasożytów złożona na pojedynczym nosicielu
Założenia 1-3 prowadzą do następującej postaci równań:
Nt +1 = Nt f (Nt, Pt )
Pt +1 = cNt[1- f (Nt , Pt )]
Jest to ogólna postać równań dla układu nosiciel  pasożyt
13
Dla skonkretyzowania dalszych rozważań przeanalizujemy konkretny
model opisany przez Nicholsona (biolog) i Baileya(fizyk) (1935).
Panowie ci dołożyli do poprzednich założeń dwa następne:
4. Spotkania pasożytów i nosicieli są przypadkowe. Liczba spotkań
jest zatem proporcjonalna do iloczynu gęstości obu populacji:
ne = aNtPt
5. Tylko pierwsze spotkanie dwóch osobników jest znaczące. (Przy
pierwszym spotkaniu pasożyt składa jaja w nosicielu, kolejne
spotkania zakażonego nosiciela z pasożytami nie zmieniają ilości
potomnych pasożytów (c), które wylęgną się z tego zakażonego
nosiciela.
Ponieważ spotkania następują przypadkowo to trzeba je opisywać
jakąś funkcją rozkładu prawdopodobieństwa. Zajście pewnej średniej
liczby zdarzeń w jednostce
Dla przypomnienia:
czasu dobrze opisuje rozkład
Prawdopodobieństwo zajścia w
Poissona.
jednostce czasu k zdarzeń których
Prawdopodobieństwo tego, że
Å›rednio zachodzi µ dane jest
nosiciel w czasie swojego
życia nie zostanie zarażony e-µµk
p(k)=
wynosi:
k!
t
f (Nt , Pt )= p(0)= e´-aP
Zatem otrzymujmy model:
t
Nt +1 = Nte-aP
t
Pt +1 = cNt(1- e-aP )
Zadanie
Znalez
ć stany stacjonarne. Zbadać ich stabilność.
Zadanie kolejne
Zbadać powyższy model przy dodatkowym założeniu, że pod
nieobecność pasożytów populacja nosicieli ma wzrost limitowany
pojemnością środowiska:
(1-Nt K)
(Nt )= er
14