Transcript mat_cw05 (ok. 547 kB)

Ćwiczenie V. „Panta rhei” - zastosowanie równań
różniczkowych do modelowania procesów biologicznych
•
Strona internetowa ćwiczeń:
http://www.home.umk.pl/~henroz/matm1112
Załóżmy, że mamy taksonomiczną grupę zwierząt lub roślin. Niech „e”,
będzie prawdopodobieństwem, że wymrą one w procesie ewolucji, a
„c” – niech będzie prawdopodobieństwem specjacji tej grupy w
ewolucji. Doświadczalnie, jest to do udowodnienia, że ogólne przystosowanie tej grupy do życia w danych warunkach („fitness”) spada
wraz z upływem czasu. Liczba zjawisk specjacji, spada w czasie, a
liczba zjawisk wymierania odpowiednio wzrasta. Załóżmy, że zależność pomiędzy obydwoma ww. zjawiskami jest liniowa, tzn., szybkości specjacji i wymierania są odpowiednio odwrotnie i wprost
proporcjonalne do czasu. Możemy dzięki temu założeniu opisać
zmianę liczby gatunków w czasie. Zmiana ta jest różnicą pomiędzy
ogółem zjawisk specjacji i wymierania. Proces taki może być
opisany przez równaniem różnicowe, lub – jeżeli weźmiemy pod
uwagę b. małe odstępy czasu – przez tzw. równanie różniczkowe.
Równanie różniczkowe wyznacza zależność pomiędzy nieznaną funkcją, a jej pochodnymi. Rozwiązanie równania różniczkowego polega
na znalezieniu funkcji f, której pochodne spełniają to równanie. W
praktyce rozwiązanie równania różniczkowego polega na sprowadzeniu go do standardowej formy, a następnie na scałkowaniu jednej
(prawej) lub obydwu stron równania. W rozpatrywanym przykładzie
zmiany liczby gatunków w obrębie danego taksonu/grupy w czasie
(w wyniku procesów ewolucji), D – oznacza zmianę, N – liczbę
gatunków, t – czas, e – prawdopodobieństwo wymierania, c –
prawdopodobieństwo specjacji: DN = Nt+1 – Nt = (ctNt – etNt)Dt; po
zamianie równania różnicowego na różniczkowe:
dN/dt = ctN – etN.
2
Czy możliwe jest obliczenie liczby gatunków w dowolnym czasie?
Tak, o ile rozwiążemy równanie względem N, tzn. gdy rozwiążemy
równanie różniczkowe. Dzięki prostemu przekształceniu (dzielenie
obu stron przez N i mnożenie przez dt) uzyskujemy: dN/N = (c – e)tdt.
Teraz całkujemy obydwie strony równania (podobnie, jak na ćw.
poprzednim) i uzyskujemy:
Jak w ćw.
c e 2
c e 2
poprzedt
t
ce 2
nim, N0 uln( N )  c1 
t  c2  N  Ce 2  N 0e 2
zyskano
2
ustawiając
t = 0. Dwa przykłady wykresów takiej funkcji zostaną pokazane na
następnym przeźroczu. Widać tu, że nawet mała różnica w szybkości
wymierania i specjacji może doprowadzić do znacznego wzrostu
liczby gatunków w obrębie taksonu lub do jego szybkiego wymarcia.
Dynamika wymierania i specjacji:
Wniosek: w przyrodzie obydwa te procesy są modyfikowane
przez dodatkowe mechanizmy,
które nie zostały uwzględnione w
ww. prostym modelu. Rozwiązane
powyżej równanie, jest liniowym
równaniem różniczkowym pierwszego rzędu. I-go rzędu, ponieważ
pochodna była I-go rzędu (I-sza) i
liniowe, ponieważ nie wystąpiły
wykładniki wyższe niż 1 przy N.
Równania różniczkowe I rzędu
są bardzo ważne w biologii i można je zapisać w formie ogólnej:
dy/dx + f(x)y = g(x), gdzie f i g są
funkcjami zależnymi od x. Równanie jest zwane homogenicznym
(jednorodnym), jeżeli g(x) = 0 dla wszystkich wartości x. Jeżeli rozpatrzymy tylko to jednorodne równanie, to po przekształceniu uzyskujemy: dy/y = –f(x)dx. Po scałkowaniu obu stron i ew. dodatkowych
przekształceniach uzyskujemy: y = Ce–f(x)dx. Stała C jest często oznaczana przez doprowadzenie x do odpowiedniej wartości (najczęściej
0) i rozwiązanie względem C.
Jeżeli równanie jest niejednorodne – sytuacja skomplikowana.
Stosujemy wtedy wzór na pochodną iloczynu: (uv)’ = u’v + uv’. Musimy
znaleźć taką funkcję a(x), przez którą moglibyśmy pomnożyć nasze
liniowe równanie różniczkowe I-go rzędu i w odniesieniu do czego
możnaby później zastosować wzór na pochodną iloczynu. Dlatego funkcja a(x), musi spełniać następujący warunek: a(x)dy/dx = a(x)f(x)y +
d(a(x)y)/dx. Równanie to ma taką samą strukturę, jak wzór na pochodną iloczynu, a ponadto: da/dx = a(x)f(x). Jest to jednorodne równanie
różniczkowe I-go rzędu, którego rozwiązanie podano powyżej. Ustalamy: a(x) = Cef(x)dx, a następnie podstawiamy to do naszego wzoru na
pochodną iloczynu i uzyskujemy:
Pamiętając, że: dy/dx + f(x)y = g(x), wnosimy, że g(x) powinna spełniać
następujący warunek:
Teraz możemy scałkować obie strony
tego równania, uzyskując:
Cef(x)dx y = Cef(x)dx g(x)dx + C1
Ostatecznie po przekształceniu tego równania przez podzielenie przez I-szą całkę:
y = e–f(x)dx ef(x)dx g(x)dx + ce–f(x)dx. Jest to
najbardziej ogólne rozwiązanie liniowego równania różniczkowego I-go
rzędu. Wygląda ono skomplikowanie, ale w większości przypadków jest
łatwe do wykorzystania, o ile funkcje f(x) i g(x) nie są zbyt skomplikowane.
Przykład praktyczny: stężenia niektórych hormonów i wielu leków w
krwioobiegu mogą być modelowane jako suma szybkości produkcji lub
indukcji i szybkości rozkładu. Załóżmy, że lek podawany jest z mniej
lub bardziej stałą szybkością „g”. Równocześnie jest on metabolizowany – proporcjonalnie do jego stężenia. Proces ten można opisać następującym równaniem: dF/dt = g – fF. Zmiana stężenia leku F w czasie t
jest zwykłą różnicą pomiędzy szybkością podawania (g) a szybkością
wykorzystywania/zużywania fF. Jest to niejednorodne (inhomogenous)
równanie różniczkowe I-go rzędu z g = f(t) = const. i c = g(t) = const.
Aby wyliczyć stężenie leku w dowolnym czasie t, wykorzystujemy rozwiązanie ogólne. Potrzebujemy: – f(t)dt = – fdt = – (ft + C1); po podsta-
wieniu: ef(t)dt g(t)dt = (eft+C1 )gdt = g(1/f)eft+C1 + C2. Podstawiamy te wyniki do rozwiązania ogólnego i uzyskujemy:
Stałe C, to wszystko stałe całkowania.
g Pozostałą stałą K można
 ft C1
 ( ft C1 )  g ft C1 
 ft
F  Ce
e
 e
  Ke  oznaczyć przyjmując t=0.
f
f

Dlatego K = F(0)–c/g.
Stąd ostateczny wynik możeg 
g
F (t )    F (0)   e ft
my zapisać następująco:
f

f 
Wynik ten jest ważny, ponieważ jest to ogólne rozwiązanie liniowego
równania różniczkowego I-rzędu, gdzie obie funkcje: f(t) i g(t) są stałymi.
Można również zastosować program matematyczny do rozwiązania naszego problemu (Mathematica, wxMaxima – cz. praktyczna). Lewa cz.
rysunku obok przedstawia zależność F(t) od czasu. Bez względu na stężenie początkowe, proces
zmierza do
stężenia stabilnego. Stężenie to jest
zwane punktem równowagi. Punkt ten wyliczamy, ustawiając dF/dt = 0. To daje nam bezpośrednio: Frównowagi = g/f, co jest widoczne na rys. Prawa cz. rys. przedstawia tzw. diagram fazowy procesu. dF/dt jest funkcją liniową i bez względu na stężenie początkowe, pierwiastkiem tej funkcji jest F(t) = g/f.
Inny przykład praktyczny: wykładnicze równanie wzrostu. Równanie to
jest realistyczne tylko dla b. niskich wielkości populacji. W rzeczywistości są pewne górne granice, zwane „zdolnościami przepustowymi”
(„carrying capacities”) i sprawiają one, że powyżej ich, rzeczywisty
wzrost populacji staje się wolniejszy dla wyższych wielkości populacji.
Rys poniżej przedstawia populację drożdży Saccharomyces cerevisiae,
hodowanych w chemostacie, które na początku doświadczenia rosły
szybko – zgodnie z wykładniczym modelem wzrostu. Lecz im wyższa
była wielkość populacji (objętość kolonii), tym wolniejszy był wzrost
drożdży. Jak modelować taki proces? Bierzemy wtedy wyjściową,
wykładniczą funkcję wzrostu i dodajemy
„drugie wyrażenie” („second term”), które
zmniejsza szybkość wzrostu populacji,
w miarę, jak N zbliża się do wartości granicznej. Jeżeli określimy tę górną granicę
wielkości populacji przez K, to cały proces
możemy modelować równaniem:
Gdy N staje się
dN
N
KN 

rN

rN

rN


coraz większe, K–N dąży do zera i drugi
dt
K 
K

czynnik [(K–N)/K] też dąży do zera i szybkość wzrostu populacji spada. Dla niskiego N, komponent (K–N)/K
jest bliski 1 i cały proces przypomina nasz wyjściowy wzrost wykładniczy. Powyższe równanie jest najprostszą formą modelowania tego
typu procesu i zwane jest logistycznym procesem wzrostu lub równaniem Verhulst’a-Pearl’a. Wykres na przeźroczu następnym pokazuje, że szybkość wzrostu jest b. mała dla małych wartości t i spada
do zera dla N=K. Maksymalna szybkość wzrostu występuje przy
N = K/2. Powyższe logistyczne równanie wzrostu nie daje nam możliwości wyliczenia N. W tym celu musimy rozwiązać odpowiednie rów-
dN
KN 
 rN 
 Aby je rozwiązać, aproksymudt
 N  jemy funkcję z jej rozwinięcia
w szereg Taylora. Zakładamy, że funkcja N(t)
daje się rozwinąć w szereg Taylora. Dlatego
dN/dt musi być funkcją algebraiczną w formie:
dN(t)/dt = a1 + a2N(t) + a3N(t)2 + a4N(t)3...
Równanie wzrostu populacji ma 2 pierwiastki:
dla N=0 i dla N=K. Najprostszym równaniem o
2 pierwiastkach jest funkcja kwadratowa. Dlatego możemy „obciąć” rozwinięcie funkcji
w szereg Taylora już po 3-cim składniku, uzyskując:
dN/dt = a1 + a2N + a3N2. Dla N = 0, uzyskujemy: a1 = 0 i dlatego możemy ten składnik opuścić: dN/dt = a2N + a3N2. Na tym etapie, dzielimy
całe równanie przez N(a2 + a3N). Dzięki temu staje się ono przekształcone:
Taka postać jest łatwiejsza do scałkowania:
nanie różniczkowe.
a3dN
dN

 a2 dt
N a2  a3 N
Stąd:
N
 Cea2t
a2  a3 N
Teraz można przekształcić to równanie, przyjmując a2/a3 = K i w ten sposób uzyskać:
W powyższym równaniu, t0 odpowiada populacji o wielkości N0 w czasię t = t0. K – „zdolność przepustowa” („the carrying capacity”), jest
górną granicą N i można ją łatwo obliczyć ustawiając dN/dt = 0. Przykład takiej funkcji jest przedstawiony na rys. poniżej. Dla przykładu z
drożdżami („2 rysunki wstecz”), udało
się
dopasować krzywą funkcji o następującym równaniu:
12, 74
N (t ) 
W genetyce lub w
1  9,32e0,26t epidemiologii, zjawiska mutacji lub infekcji
są zjawiskami masowymi. Oznacza to, że
efekt końcowy określany jest przez całkowitą liczbę czynników indukujących mutację lub infekcję. Dlatego
powinniśmy znać sumę wszystkich tych czynników. Zakładamy, że cały
badany proces daje się opisać za pomocą logistycznego modelu wzrostu. Stąd, całkowita liczba czynników może być równa całce, której odpowiada pole powierzchni pod krzywą funkcji logistycznej. Wyliczamy
to za pomocą programu matematycznego,
uzyskując: Ostatni przykład
K
K ln(1  e a (t t0 ) )
zmierza do ogólnego rozwiązania
dt


C
kwadratowego równania różnicz1  e a (t t0 )
a
kowego I-szego rzędu. Możemy
takie równanie zdefiniować następująco: dy/dx = ay + by2. Na początku
oznaczamy punkt równowagi i uzyskujemy dla y’=0: y(x) = –a/b.

Korzystamy znów z programu matematycznego (Mathematica, wxMaxima) i uzyskujemy rozwiązanie ogólne: Wynik – tzw. uogólnione logisrównanie wzrostu, a zarazem – rozwiązanie tzw.
ae ac e ax tyczne
autonomicznego równania różniczkowego. Można sprawy
ac ax
1  be e dzić, że dla wysokich wartości x, funkcja ta asymptotycznie zbliża się do wartości stanu równowagi: –b/a. Przy
użyciu Mathematica, można przekształcić to rozwiązanie i przedstawić
je w kategoriach funkcji sinh i cosh.
Inny przykład z genetyki i ekologii: załóżmy, że mamy populację gatunków zwierzęcych lub roślinnych, które zasiedlają „skrawki siedliska” („habitat patches”)(rys. – poniżej). Przyjmujemy, że wyjściowa częstość zasiedlonych fragmentów siedliska wynosi
p, a wyjściowa częstość niezasiedlonych fragmentów siedliska, to q = (1-p). Częstości są zawsze wybierane w taki sposób, aby ich całkowita
suma = 1 (lub 100%). Członkowie populacji migrują. Są oni zdolni do kolonizowania niezasiedlonych fragmentów, ew populacje na zasiedlonych fragmentach wymierają. Taka populacja,
która jest rozdzielona na szereg „skrawków/fragmentów” siedliska, powiązanych ze sobą zjawiskami migracji, nazywana jest metapopulacją. Problemem do rozwiązania jest modelowanie
zmian p (liczba zasiedlonych fragmentów) w zależności od czasu. Zmiana p w czasie, to: dp/dt. Definiujemy teraz m, jako tempo kolonizacji
(imigracji) i v – lokalne tempo wymierania (obie wartości wyrażane
jako
prawdopodobieństwo). Tak więc zmiany w p mogą być teraz modelowane jako suma 2 niezależnych procesów: zmian spowodowanych kolonizacją i zmian wywołanych lokalnym wymieraniem. Zakładamy, że
wzrost liczby zasiedlonych skrawków jest proporcjonalny do faktycznej
liczby zasiedlonych skrawków (p) i do liczby pustych skrawków q=1-p.
Stąd: dp/dt = mp(1–p). Z kolei spadek p, powinien być proporcjonalny do
lokalnej szybkości wymierania. Stąd: dp/dt = –vp. Teraz możemy modelować cały proces za pomocą prostego równania: dp/dt = mp(1–p) – vp =
–mp2 + p(m – v). Jest to bdb. znany model metapopulacji Richarda
Levinsa, zaproponowany w 1969 r. Aby oznaczyć liczbę zasiedlonych
skrawków dla stanu równowagi, ustawiamy: dp/dt = 0 i uzyskujemy
bezpośrednio: p = 1 – v/m. Metapopulacja utrzymuje się przy życiu tylko
wtedy, gdy: v < m. Załóżmy, że v jest odwrotnie proporcjonalne do powierzchni fragmentu siedliska (A). Oznacza to, że na większych skrawkach, prawdopodobieństwo wymierania spada, ze względu na większą
populację. Stąd v a 1/A. Załóżmy dalej, że szybkość kolonizacji (m) spada wraz ze wzrostem odległości pomiędzy skrawkami (D). Stąd: m a 1/D.
Uzyskujemy ostatecznie: dp/dt = (1/D)p(1–p) – (1/A)p. Wniosek I-szy:
liczba zasiedlonych siedlisk będzie wzrastała wraz ze wzrostem średniej powierzchni skrawka (z powodu spadku szybkości wymierania).
Wn. II-gi: liczba zasiedlonych siedlisk będzie spadała wraz ze wzrostem
odległości pomiędzy skrawkami (ze względu na ograniczone procesy
kolonizacji). Przedostatnie równanie jest wg przyjętej terminologii
kwadratowym równaniem różniczkowym I-go rzędu. Poniżej – podane
ogólne rozwiązanie równania ostatniego:
Brak jest jednak początkowej wartości C. Przy t = 0,
p = p0 – częstość początkowa. Biorąc to pod uwagę, uzyskujemy C:
Logarytm w liczniku musi być > 0.
Stąd właściwe rozwiązanie musi
spełniać warunek: p0 > 1 – v/m. Poniżej, przedstawiony jest wykres
zależności proporcji zasiedlonych skrawków [p(t)] w zależności od czasu: . Tylko 10 pokoleń wystarcza do osiągnięcia stanu równowagi przy:
p(t) = 1 – 0,3/0,5 = 0,4. Stąd na podstawie warunków początkowych można wnosić, że 40% fragmentów/skrawków siedliska zostanie zasiedlonych. Nie zawsze jednak jest to tak
proste i całki albo przyjmują bardzo
skomplikowaną postać lub nawet nie
można wyrazić całki w tzw. formie zamkniętej. Można spróbować metody uzyskiwania przybliżenia całek. Aby to zrobić, rozwijamy funkcję w szereg Taylora
i „obcinamy” szereg przy odpowiednim
wyrazie. Całkę można wtedy łatwo wyliczyć ręcznie, ponieważ będzie
potrzebna następująca całka:
Dla logistycznego równania wzrostu, program Mathematica wylicza
a
a( x  b) n dx 
( x  b) n 1  C automatycznie rozwinięcie w szereg Taylora. „Obcinamy” po 3-cim
n 1
wyrazie i uzyskujemy następujące
przybliżenie całki:
x x
x
x
 1  x dx   1  2  8 dx  x  8  56 |  0, 492

0,5
0,5
3
3
6
4
7
0,5
0
K
K aK
a3 K
K
aK
a3 K
3
2
4
dt


(
t

t
)

(
t

t
)
dt

(
t

t
)

(
t

t
)

(
t

t
)
C
0
0
0
0
0
 1  ea (t t0 )  4 4
48
4
8
192
0
0
Można porównać wyniki: dokładne rozwiązanie dało dla a = 1, K = 1
i t0 = 0 pomiędzy t = 0 a t = 3 wartość N = 2,355. Przybliżenie dało:
N = 1,45. Stąd nasz błąd wynosi: (2,355 – 1,45)/2,355 = 0,38 = 38% (..nie
najlepiej.....!).
Przy aproksymacji z wykorzystaniem rozwinięcia funkcji w szereg
Taylora, wynik znacznie zależy od tego jak szybko szereg osiąga zbieżność; szybka zbieżność oznacza wyniki bliskie dokładnego rozwiązania. Można to sprawdzić np. sporządzając wykres zależności sumy wyrazów szeregu od kolejnych numerów wyrazów.
W powyższym przykładzie aproksymacja do szeregu Taylora zawiodła. Dlaczego? Ponieważ szereg jest naprzemienny (wyrazy o znakach:
+, –, + –, +, –, etc.). „Obcięcie” za wyrazem o niewielkim numerze kolejnym (porządkowym) pociąga za sobą bardzo wysokie błędy. Praktyka
sugeruje, że powinniśmy wykorzystywać jako aproksymacje tylko te
rozwinięcia w szereg, które zawierają wyrazy wyłącznie dodatnie lub
ujemne, nadają się do wykorzystania jako aproksymacje funkcji. Rozwinięcia wielu znanych funkcji często dają szeregi, które „z trudem są
zbieżne” (wyjątek – prosta funkcja wykładnicza: y = ex). Inny przykład:
funkcja y = (1 – x3)0,5, w zakresie: 0 < x < 1. Mathematica wylicza b.
skomplikowaną całkę. Jednakże szereg Taylora jest bardzo prosty i
osiąga zbieżność bardzo szybko. Całka w granicach od 0 do 0,5 wynosi:
0,5

0
0,5
3
6
4
7
0,5
x
x
x
x
3
1  x dx   1   dx  x   |  0, 492
2 8
8 56 0
0
Rozwiązanie numeryczne przy użyciu Mathematica dało prawie identyczny wynik: 0,492041 + 3,33067  10–16.
Wskazówki do wykonania zadań praktycznych ćw. 5.
• Wskazówki do zadania 1:
Ogólna postać równania różniczkowego, to: dy/dx = f(x)y. Przystępując do jego rozwiązania, zazwyczaj dzielimy obie jego strony przez
y i mnożymy przez dx, wskutek czego uzyskujemy: dy/y = f(x)dx.
Całkujemy teraz obie strony, uzyskując: ln(y) = f(x)dx + C. Po przekształceniu: y = Cef(x)dx. C można wyznaczyć, przyjmując x0 = 0;
wtedy (o ile wykładnik przy e zostanie sprowadzony do 0): C = y0 i
ostatecznie: y = y0ef(x)dx. Dla I-szego równania:
dy/dx = ay
| :y *dx
dy/y = adx
dy/y = adx
ln(y) = ax + C
y = Ceax = y0eax.
dy/dx = y/x
dy/y = (1/x)dx
dy/y = (1/x)dx
ln(y) = ln(x) + C
y = Cx.
Dla II-go równania:
| :y *dx
Wskazówki do zadania 2:
Uruchamiamy program wxMaxima. Pierwsze równanie różniczkowe:
dy/dx = ex–y, wprowadzamy do pola „INPUT:” następująco:
‘diff(y,x)=%e^(x-y), co widać na zrzucie ekranu: Następnie naciskamy <Enter>, uzyskując:
Po tym, wracamy do pola przyciskami
i klikamy w przycisk „Solve ODE”,
który otwiera okno dialogowe rozwiązywania równań różniczkowych:
Klik
Wynik – okno dialogowe – na następnym przeźroczu
Okno dialogowe rozwiązywania równań różniczkowych:
Akceptujemy wszystkie opcje domyślne
i klikamy w przycisk OK. Znak „%” w
w polu „EquaKlik
ion” oznacza
wzięcie do obliczeń poprzedniego wyniku.
W następstwie tego, uzyskujemy rozwiązanie:
Nie jest to jednak rozwiązanie typu:
y = f(x), o jakie nam chodzi.
Aby uzyskać tego typu rozwiązanie,
należy uruchomić w wxMaxima
komendę: „Solve(y)”. Można ją wprowadzić bezpośrednio do pola:
„INPUT:”. Można też kliknąć w przycisk „Solve”, a po ukazaniu się okna
dialogowego rozwiązywania równań,
zamienić x na y w polu „for variable(s):” (następne przeźrocze).
Okno dialogowe rozwiązywania równań (nie różniczkowych!!!!!):
Zamieniamy: x  y, uzyskując:
Klik
Ostatecznie, uzyskujemy wynik:
Wynik ten, należy odczytać jako:
y = ln(ex + C).
Drugie równanie różniczkowe: dy/dx = (y+1)/x; wszystkie czynności
przeprowadzamy analogicznie, jak przy rozwiązywaniu równania
poprzedniego (poprawne wprowadzenie: ‘diff(y,x)=(y+1)/x) – z tym,
że pomijamy ostatni etap „solve(y)” (3 ostatnie zrzuty ekranu dla
równania poprzedniego), uzyskując od razu właściwe rozwiązanie:
Rozwiązanie to odczytujemy:
y = (C – 1/x)x.
Przy rozwiązywaniu trzeciego równania
różniczkowego: dy/dx = xy2 (r-nie kwadratowe; popr. wprowadzenie:
‘diff(y,x)=x*y^2), postępujemy analogicznie, jak w przypadku r-nia I-go:
Rozwiązanie to, należy odczytać:
2
y 2
x  2C
Czwarte równanie różniczkowe: dy/dx = –ayn; czynności rozwiązywania ogólnie wykonujemy tak, jak dla równań poprzednich (poprawne
wprowadzenie: ‘diff(y,x)=–a*y^n).
Jednakże po uruchomieniu „Solve ODE”, pojawia się pytanie: „Is
n – 1 zero or nonzero?”
Odpowiadamy w polu:
„INPUT:” „nonzero” I
wciskamy <Enter>. Pojawia się kolejne pytanie: „Is n an integer?”
(„Czy n jest liczbą całkowitą?). Odpowiadamy
: „yes” <Enter>.
Po tym pojawia
się ostateczne rozwiązanie, które odczytujemy:
y  (anx  ax  Can  Ca)
1
1n
Piąte równanie różniczkowe: dy/dx = axy – bxy; powtarzamy czynności rozwiązywania, jak w przypadku równań poprzednich. Poprawne
wprowadzenie: ‘diff(y,x) = a*x*y – b*x*y. Nie jest konieczne korzystanie z opcji „Solve(y)”.
Ostateczne rozwiązanie:
y  Ce
(a  b) 2
x
2
Biologiczne zastosowanie powyższego równania różniczkowego: modelowanie procesów ewolucji (specjacja,
wymieranie..) dowolnego taksonu organizmów żywych (zwierzęta
lub rośliny) – przedstawione na początku części teoretycznej. Znaczenie symboli (zmiana): y N (liczba gatunków w czasie t), C  N0
(stała C: liczba gatunków w czasie t0 = 0), a  c [prawdopodobieństwo (częstość) specjacji], b  e [prawdopodobieństwo (częstość)
wymierania], x  t (czas). Po zamianie:
 c  e  2  .
N  N0 exp 
t 

 2  
Wskazówki do zadania 3:
Model procesu ewolucji: N = N0exp[((c – e)/2)t2], został przedstawiony
(włącznie z podaniem opisów zmiennych i parametrów) na przeźroczu
poprzednim. Plik „ewolucja1.xls”, dostępny ze stron internetowych
ćwiczeń, umożliwia automatyczne zmiany podstawowych parametrów
równania funkcji N(t): N0, c i e.
Po pobraniu i zapisie pliku, dokończeniu obliczeń i sporządzeniu wykresu punktowego (XY), możliwe są łatwe zmiany parametrów równania i
obserwacje wynikłych stąd zmian krzywych N(t). Najważniejszą rolę
odgrywają tu zmiany wartości c i e – a w szczególności – wzajemna
proporcja c do e. Przykładowe zrzuty ekranów – patrz dalej.
Przykłady – funkcja rosnąca (1, 2 i 3) i stała (4)
Przykłady, cd. – funkcja malejąca:
Przy stałym N0, parametry równania mają następujący wpływ na
przebieg krzywej N(t):
c > e: funkcja rosnąca
c = e: funkcja stała (N = N0)
c < e: funkcja malejąca.
Wskazówki do zadania 4:
Logistyczne równanie wzrostu: N(t) = K/[1+C*exp(–a2t)], zostało przedstawione (włącznie z podaniem opisów zmiennych i parametrów) w
podstawowej instrukcji do ćwiczenia na WWW. Plik „logist1.xls”,
dostępny ze stron internetowych ćwiczeń, umożliwia automatyczne
zmiany podstawowych parametrów równania funkcji N(t): K, C i a2.
Po pobraniu i zapisie pliku, dokończeniu obliczeń i sporządzeniu wykresu punktowego (XY), możliwe są łatwe zmiany parametrów równania
i obserwacje wynikłych stąd zmian krzywych N(t). Najważniejszą rolę
odgrywają tu zmiany wartości C i a2 – a w szczególności – wzajemna
proporcja C do a2. Przykładowe zrzuty ekranów – patrz dalej.
Przykłady (1-4):
Przykład 5 (ostatni):
Wpływ parametrów równania logistycznego na kształt krzywej
(z wyjątkiem K):
C – z 1 strony decyduje o tym,
„która część krzywej” będzie
widoczna na wykresie (jeśli wysokie – to początkowa; jeśli niskie
– końcowa), a z II-giej – o wartości początkowej liczby organizmów (N0): im niższe – tym wyższa.
Wynika to z zależności:
C = K/N0–1.
Stała a2 decyduje o tym, jak szybko krzywa dąży do swej wartości
maksymalnej (asymptotycznie do K): im wyższa, tym szybciej (~współczynnik kierunkowy krzywej logistycznej).
Wskazówki do zadania 5:
Model opisujący podawania glukozy do krwioobiegu w zależności od
czasu G(t) = c/a + [G(0) – c/a]exp(–at), zostało przedstawione (włącznie
z podaniem opisów zmiennych i parametrów) w podstawowej instrukcji
do ćwiczenia na WWW. Plik „glukoza1.xls”, dostępny ze stron
internetowych ćwiczeń, umożliwia automatyczne zmiany
podstawowych parametrów równania funkcji G(t): G(0), c i a.
Po pobraniu i zapisie pliku, dokończeniu obliczeń i sporządzeniu wykresu punktowego (XY), możliwe są łatwe zmiany parametrów równania
i obserwacje wynikłych stąd zmian krzywych G(t). Najważniejszą rolę
odgrywają tu zmiany wartości c i a – a w szczególności – wzajemna
proporcja c do a. Przykładowe zrzuty ekranów – patrz dalej.
Przykłady – funkcja rosnąca (1, 2), stała (3) i malejąca (4)
Przykład 5 – funkcja malejąca:
Wpływ parametrów równania
na kształt krzywej:
parametry równania mają następujący wpływ na przebieg krzywej G(t):
G(0) > c: funkcja malejąca
G(0) = c: funkcja stała [G = G(0)]
G(0) < c: funkcja rosnąca
Bez względu na wartości G(0),
funkcja asymptotycznie dąży
do wartości równej:
G(t) = c/a (dla wysokich wartości
t)
Wskazówki do zadania 6:
Model metapopulacji, opisujący zmiany w częstości zasiedlonych
skrawków ekosystemuw zależności od czasu:
p(t)=(1–(v/m))/[(1–exp((m–v)C)*exp((v–m)t))], został przedstawiony
(włącznie z podaniem opisów zmiennych i parametrów) w podstawowej
instrukcji do ćwiczenia na WWW. Plik „metapop1.xls”, dostępny ze
stron internetowych ćwiczeń, umożliwia automatyczne zmiany
podstawowych parametrów równania funkcji p(t): C, m i v.
Po pobraniu i zapisie pliku, dokończeniu obliczeń i sporządzeniu wykresu punktowego (XY), możliwe są łatwe zmiany parametrów równania
i obserwacje wynikłych stąd zmian krzywych p(t). Najważniejszą rolę
odgrywają tu zmiany wartości m i v – a w szczególności – wzajemna
proporcja m do v. Przykładowe zrzuty ekranów – patrz dalej.
Przykład a) wysokie m i niskie v:
Przy wysokim m i niskim v – populacja nie wymiera całkowicie, lecz
dąży do stanu równowagi: p(t) = 1 – v/m. W konkretnym przypadku, stan
równowagi jest osiągany dla: p(t) = 1 – 0,1/0,9 = 1 – 0,1111 = 0,8889.
Przykład b) niskie m i wysokie v:
Przy niskim m i wysokim v – populacja wymiera całkowicie, po odpowiednio długim czasie (w niniejszym przykładzie, model osiąga wartości skrajnie niskie dla względnie wysokich wartości czasu).
Dziękuję
za uwagę ;-)