Zagadnienia

9.1 Sformułowanie zadania i implementacja
9.2 Weryfikacja wyników

9. Równanie logistyczne z opóźnieniem

Klasyczne modele oparte na równaniach różniczkowych zwyczajnych zakładają, że w danej chwili reakcja systemu (zmiana wartości) jest natychmiastowa i zależy od jej wartości w tym momencie (czyli, że $y^{{\prime}}(t)=f(t,y(t))$ . Jednak bywają sytuacje, kiedy takie założenie jest nierealistyczne i należy uwzględnić fakt, iż — wskutek swego rodzaju bezwładności układu — wartość pochodnej rozwiązania może zależeć nie tylko od jego bieżącej kondycji, lecz również od wartości z chwil wcześniejszych. W pierwszym przybliżeniu taka koncepcja prowadzi to do tzw. równań różniczkowych ze stałym opóźnieniem $\tau\geq 0$ :

$y^{{\prime}}(t)=f(t,y(t),y(t-\tau)),\qquad t\in(t_{{\min}}t_{{\max}}].$

Pakiet OdePkg, rozszerzający możliwości Octave i dostępny w Octave-Forge [20], zawiera w sobie kilka solverów równań różniczkowych ze stałym opóźnieniem $\tau$ . Drobną ich wadą jest to, że dopuszczają jedynie stały warunek początkowy (w równaniach z opóźnieniem warunek początkowy musi być zadany nie w jednym punkcie, ale na odcinku, $[t_{{\min}}-\tau,t_{{\min}}]$ . Wszystkie solvery wykorzystują tzw. włożone jawne metody Rungego–Kutty, od metody niskiego rzędu ode23d, po metodę wysokiego rzędu ode78d. Warto pamiętać, że w przypadku równań z opóźnieniem rozwiązania mogą nie być zbyt gładkie i wówczas korzystniej może być stosować schemat niższego rzędu.

Sposób wykorzystania pakietu jest bardzo dokładnie opisany w [25], na którym, jak się wydaje, opiera się zestaw solverów z OdePkg. My ograniczymy się do krótkiego omówienia sposobu rozwiązania jednego z najprostszych równań.

9.1. Sformułowanie zadania i implementacja

Rozważmy równanie postaci

	$\displaystyle y^{{\prime}}(t)$	$\displaystyle=(a-y(t-1))\cdot y(t),\qquad\text{ dla }t\in(0,T],$		(9.1)
	$\displaystyle y(t)$	$\displaystyle=0.2,\qquad\text{ dla }t<0.$		(9.2)

Jest to tzw.równanie logistyczne z opóźnieniem. Dodatni współczynnik $a$ jest bardzo istotnym parametrem. Okazuje się bowiem, że gdy $a<\pi/2$ , to — dla dostatecznie długich czasów — rozwiązania (9.1) stabilizują się na poziomie $a$ , natomiast dla $a>\pi/2$ rozwiązania wpadają w niegasnące, regularne oscylacje. Więcej na temat logistycznego równania z opóźnieniem można dowiedzieć się z wykładu Modelowanie matematyczne w biologii i medycynie.

Zacznijmy od zapisania skryptu, rozwiązującego nasze równanie dla wartości parametru $a$ równych $0.1,\pi/2-0.1,\pi/2+.01$ , na przedziale czasowym długości $T=100$ (dostatecznie długim, byśmy zaobserwowali prawdopodobną tendencję rozwiązania).

global a;

lag = 1;
T = 100;

% parametry solvera
opts = odeset ('Stats', 'on', ...
	'RelTol', 1e-4, ... % tolerancja bledu wzglednego
	'AbsTol', 1e-4, ... % tolerancja bledu bezwzglednego
	'InitialStep', 1e-5, ... % poczatkowy krok schematu
	'MaxStep', 1e-2); % przy okazji steruje liczba punktow probkowania rozwiazania
	
filename = 'logdel';

for a = [0.1, pi/2 - 0.1, pi/2 + 0.1]
	
	sol = ode23d(@logdelf, [0,T], 0.2, lag, 0.2, opts);
	t = sol.x; Y = sol.y; sol.stats
	Ylag = interp1(t,Y,t-lag,'cubic'); % rozwiazanie opoznione
	Ydot = logdelf(t,Y,Ylag);  % pochodna rozwiazania
	
	plot(t,Y);
	xlabel('t'); ylabel('y(t)'); title('Rozwiazanie y(t)');
	print('-depsc', strcat(filename, '-a-', ...
		strrep(num2str(a,'%g'),'.','p'), '.eps') )
	
	pause(1);
	
	plot(Y,Ydot);
	xlabel('y(t)'); ylabel('y''(t)');  title('Wykres fazowy');
	print('-depsc', strcat(filename, '-a-', ...
		strrep(num2str(a,'%g'),'.','p'), '-phase.eps') )
	
	pause(1);
		
end

Kluczowa jest linijka

sol = ode23d(@logdelf, [0,T], 0.2, lag, 0.2, opts);

która w rzeczywistości rozwiązuje nasze równanie. Parametrami funkcji ode23d są kolejno:

funkcja określająca prawą stronę²³Funkcja musi być przekazana jako wskaźnik, a nie jako nazwa! (u nas jest to logdelf), do której parametr a przekazujemy z przyzwyczajenia za pomocą zmiennej globalnej²⁴W rzeczywistości funkcje z pakietu OdePkg umożliwiają wygodne przekazywanie listy dodatkowych parametrów do wnętrza funkcji prawej strony.),
przedział czasowy, na którym wyznaczamy rozwiązanie, tutaj: $[0,T]$ ,
warunek początkowy,
wartości opóźnień (u nas jest tylko jedno, bo równanie jest skalarne) lag = 1), no i ewentualnie
struktura, zmieniająca domyślne parametry pracy solvera.

function ydot = logdelf(t,y,ylag)
global a;
ydot = (a - ylag).*y;
end

W odróżnieniu od lsode, ode23d zwraca wartości rozwiązania w przez siebie wyznaczonych punktach przedziału $[0,T]$ . Aby więc określić jego wartość w zadanych przez nas punktach,

t = linspace(0,T,10*T);

musimy interpolować wartości; warto także wcześniej zadać w opcjach maksymalną długość kroku, która zagwarantuje nam wystarczającą rozdzielczość próbkowania; dla danych opóźnionych, a następnie wyznacz

Ylag = interp1(t,Y,t-lag,'cubic'); % rozwiązanie opóźnione
Ydot = logdelf(t,Y,Ylag);  % pochodna rozwiązania

która zwraca nam wartości obliczonego przybliżenia rozwiązania (i jego pochodnej!) w punktach $t$ . Potem wystarczy je tylko narysować, komendami

plot(t, Y);
plot(Y, Ydot);

— najpierw rozwiązanie w zależności od czasu, potem — wykres fazowy.

Wyniki, które, przynajmniej w sposób jakościowy, zgadzają się z przewidywaniami teoretycznymi, przedstawiają poniższe wykresy.

Rys. 9.1. Wykres rozwiązania równania logistycznego z opóźnieniem dla parametru $a$ równego $1.4708<\pi/2$ .

Krzywa fazowa jest przyciągana do jednego punktu

Rys. 9.2. Wykres fazowy rozwiązania równania logistycznego z opóźnieniem dla parametru $a$ równego $1.4708<\pi/2$ .

Rozwiązanie wpadają w niegasnące oscylacje wokół $a$

Rys. 9.3. Wykres rozwiązania równania logistycznego z opóźnieniem dla parametru $a$ równego $1.6708>\pi/2$ .

Krzywa fazowa nawija się na pewną orbitę.

Rys. 9.4. Wykres fazowy rozwiązania równania logistycznego z opóźnieniem dla parametru $a$ równego $1.6708>\pi/2$ .

9.2. Weryfikacja wyników

Aby nabrać więcej przekonania do ilościowych wyników, zobaczymy, czy stukrotne zmniejszenie parametrów tolerancji błędu spowoduje wyraźną zmianę uzyskanych wartości rozwiązania. Sprawdzenie przeprowadzimy dla parametru $a=\pi/2+0.1$ .

Najpierw wyznaczymy w standardowy sposób wartości rozwiązania w interesujących nas punktach:

a = pi/2 + 0.1;

opts = odeset ('Stats', 'on', ...
	'RelTol', 1e-4, ... % tolerancja błędu względnego
	'AbsTol', 1e-4, ... % tolerancja błędu bezwzględnego
	'InitialStep', 1e-5, ... % początkowy krok schematu
	'MaxStep', 1e-2); % przy okazji steruje liczbą punktów próbkowania

sol = ode23d(@logdelf, [0,T], 0.2, lag, 0.2, opts);
t = sol.x; Y = sol.y; sol.stats

a następnie powtórzymy obliczenia dla ostrzejszych kryteriów:

opts = odeset ('Stats', 'on', ...
	'RelTol', 1e-6, ... % tolerancja błędu względnego
	'AbsTol', 1e-6, ... % tolerancja błędu bezwzględnego
	'InitialStep', 1e-7, ... % początkowy krok schematu
	'MaxStep', 1e-2); % przy okazji steruje liczbą punktów próbkowania

sol = ode23d(@logdelf, [0,T], 0.2, lag, 0.2, opts);
tref = sol.x; Yref = sol.y; sol.stats

Jednak kłopot w tym, że — w przeciwieństwie do lsode — funkcja ode23d zwraca rozwiązanie nie we wskazanych przez nas wartościach $t$ , ale raczej zwraca wszystkie wyznaczone przez siebie wartości, co oznacza, że dla różnych parametrów wywołania solvera prawie na pewno dostaniemy inny zestaw sol.x.

Dlatego musimy napierw sprowadzić wyznaczone wartości do wspólnego zestawu chwil czasowych $t_{0},\ldots,t_{N}$ . Wydaje się rozsądne, by za ten bazowy zestaw wziąć po prostu t — wartości wyznaczone przy mniej restrykcyjnych parametrach pracy solvera. Następnie, musimy interpolować²⁵To może wprowadzić dodatkowe zaburzenia w wynikach, więc musimy być czujni! wartości (tref,Yref) na t, do czego idealnie nadaje się funkcja interp1.

Yreft = interp1(tref, Yref, t);

Następnie sprawdzamy, jak bardzo różnią się rozwiązania Y i Yref:

octave:> norm(Yreft-Y, Inf)
ans = 4.8362e-06
octave:> semilogy(abs(Yreft-Y));

A więc maksymalna różnica pomiędzy rozwiązaniami wynosi jedynie około $10^{{-6}}$ , co jest raczej uspokajającą wiadomością (choć należałoby jeszcze wykluczyć na przykład zjawisko aliasingu, które mocno dało się nam we znaki, gdy rozwiązywaliśmy zadanie o transformatorze, z rozdziału 6).

Ćwiczenie 9.1

Potwierdź eksperymentalnie, że właśnie dla $a=\pi/2$ następuje zmiana długoterminowego charakteru rozwiązań, z gasnących do zera, na oscylacyjne.

Ćwiczenie 9.2

Napisz funkcję, pozwalającą rozwiązać równanie logistyczne z opóźnieniem o zadanych parametrach i jednocześnie od razu zweryfikować w opisany powyżej sposób jakość rozwiązania. Zadbaj o komunikaty i wizualizacje dla niecierpliwego użytkownika. Następnie sprawdź pozostałe uzyskane powyżej wyniki.

Ćwiczenie 9.3

Napisz skrypt Octave lub MATLABa, rozwiązujący inną wersję równania logistycznego z opóźnieniem, pochodzącą z pewnego modelu wzrostu nowotworu u myszy [4]:

$N^{{\prime}}(t)=rN(t-\tau)\left(1-\dfrac{N(t-\tau)}{K}\right).$

Dodatnie współczynniki $r$ i $K$ są parametrami modelu. Zweryfikuj jakość uzyskanych aproksymacji rozwiązania i jego pochodnej w przypadku, gdy $r=4$ , $K=2$ , a opóźnienie wynosi $\tau=0.7$ . Jako funkcję początkową przyjmij stałą równą $N_{0}(t)=1.2$ .