Rozważmy równanie postaci

Jest to tzw.równanie logistyczne z opóźnieniem. Dodatni współczynnik a jest bardzo istotnym parametrem. Okazuje się bowiem, że gdy a<π/2, to — dla dostatecznie długich czasów — rozwiązania (9.1) stabilizują się na poziomie a, natomiast dla a>π/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,π/2-0.1,π/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);

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.

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=π/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

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

Dlatego musimy napierw sprowadzić wyznaczone wartości do wspólnego zestawu chwil czasowych t0,…,tN. 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).