Na tym etapie nie wydaje się, by powyższe równanie mogło sprawić nam (to znaczy: naszemu systemowi obliczeń numerycznych) jakikolwiek poważniejszy problem. Wystarczy skorzystać z gotowych narzędzi rozwiązywania równań różniczkowych zwyczajnych dostępnych w Octave (lub w MATLABie): na przykład, z solvera ode45 zaimplementowanego w obu wspomnianych systemach obliczeń naukowych.

W naszych obliczeniach przyjmiemy na początek następujące wartości liczbowe parametrów zadania:

Rozwiązanie problemu i wizualizacja składowej yN mogłaby być na przykład zrealizowana poniższym, całkiem standardowym kodem:

Jeśli chcesz pracować w serwisie internetowym WebOctave, najpierw musisz zdefiniować funkcje prawej strony rhs oraz v0.

function dY = rhs(t,Y)
% prawa strona rownania
N = length(Y);
q = [y0(t); Y(1:N-1)];
dY = -Y*(1+Y(N)) + q;
end

function v = y0(t)
% wymuszenie zewnetrzne; t musi byc skalarem!
global tau;
global Ti;
v = (min(abs(t - Ti)) < tau); % krótkie ($2\tau$) pulsy wokół $T_i$
end

Program działa w Octave i wydaje się, że produkuje poprawne wyniki (yN stabilizuje się), ale już dla T=153 obserwujemy ciekawe zjawisko: wyznaczone rozwiązanie wykazuje zachowanie oscylacyjne!

A więc mamy ,,numeryczny dowód”, że oczy skorupiaka (przynajmniej w ramach naszego modelu), poddane pulsacyjnemu działania światła zewnętrznego, reagują także w sposób periodyczny — co zasadniczo nie powinno nas dziwić… Warto więc byłoby zbadać okres tej reakcji, a następnie spróbować odnieść to zachowanie do wyników eksperymentów fizycznych. Jednak… fizyczny eksperyment nie może być przeprowadzony w tak długim czasie, pozostaje więc nam poruszać się na gruncie teorii.

Po rozwiązaniu zagadki oscylacji — nieistniejących w rzeczywistości — możemy teraz mieć ochotę na zabawę z wartością parametru τ, określającego długość trwania impulsu świetlnego. Z punktu widzenia modelowania, interesujące byłoby wyznaczenie wartości τ, przy której bezwładność oczu stwora jest na tyle duża, że przestają one reagować na (najwyraźniej zbyt krótkie) impulsy świetlne i ich reakcja wygasa.

W tym celu wystarczy stopniowo zmniejszać τ i obserwować, jak zmieni się odpowiedź układu; dla wystarczająco krótkich impulsów spodziewamy się, że bodźce będą zbyt słabe, by stabilnie wzbudzić niezerową wartość składowej yN: spodziewamy się yN gasnących do zera.

Przeprowadźmy zatem symulację dla τ=10-2 — to znaczy dwudziestokrotnie mniejszej niż dotychczas.

Jak można zorientować się z przeprowadzonej symulacji, początkowo yN rośnie (eksponencjalnie?) i nie wiadomo, czy potem stabilizuje się, jak poprzednio, czy też może zachowanie jest inne. Dlatego należy wydłużyć czas symulacji:

Okazuje się, że faktycznie, rozwiązania yN po chwilowym wzbudzeniu szybko gasną do zera. Aby upewnić się, czy nie mamy czasami do czynienia z numerycznym artefaktem, rozwiążemy to samo równanie, drastycznie zaostrzając kryteria tolerancji błędu, używając do tego funkcji odeset.

Jak widać, choć otrzymane krzywe nie pokrywają się idealnie, zgadzają się ze sobą co do charakteru: yN, po chwilowym wzroście, gaśnie do zera.

Gdybyśmy jednak wykazali mniej zaufania do wyników i rozwiązali zadanie… innym solverem, na przykład ode78 (dostępny w Octave, w MATLABie możesz wybrać na przykład ode115)?

Może okazać się, że otrzymamy zupełnie inny wykres! Co prawda pragmatycy staraliby się podkreślić, że ,,co do charakteru” dalej mamy zgodność, bo rozwiązania gasną (do zera?), ale większość z nas, widząc wykresy jak na rysunku 12.2, zaczęłaby powątpiewać w jakość prowadzonych obliczeń.

Może jednak problem leży w tym, że impulsy są na tyle krótkie, że nasz solver, stosując pewne strategie doboru kroku całkowania ,,opuszcza” część impulsu tak, że to z tego powodu rozwiązania gasną?

Gdyby bowiem zastosować jeszcze bardziej drastyczne parametry tolerancji błędu, na przykład

odeset('RelTol',1e-12,'AbsTol', 1e-12)

Po doświadczeniu z większymi τ bylibyśmy zapewne skłonni uwierzyć bez zastrzeżeń, że otrzymaliśmy wynik wykazujący, że dla τ=10-2 rozwiązania jednak nie zanikają w czasie.

Ale, jeśli przeprowadzilibyśmy symulację korzystając z dokładniejszej metody, ode78 i równie ostrych parametrach tolerancji błędu jak ostatnio, to zapewne zaskoczyłby nas przebieg tak wyznaczonego rozwiązania, przedstawiony na rysunku 12.4.

Bardzo możliwe, że jednak za bardzo przejęliśmy się narzędziami, a za mało samym problemem: przecież rozwiązania, które nas interesują, są bardzo mało regularne! A to z kolei znaczy, że stosowanie solvera wysokiego rzędu aproksymacji — takiego jak ode78 — nie bardzo ma sens; rzeczywiście, dla parametrów tolerancji błędu 10-14, wyniki symulacji (prawie) pokrywają się (por. rysunki 12.5 i 12.6: nie tylko dla obu dotychczas przez nas próbowanych solverów, ale także dla lsode, który wykorzystuje nie tylko zmiany długości kroku, ale także dostosowuje rząd metody.

A jak Ty uważasz, poniżej jakiej wartości τ rozwiązania yN naprawdę dążą do zera, gdy t→∞?

Zauważmy, że nasze zadanie jest mimo wszystko trudne do rozwiązania numerycznego! Wiąże się to z charakterem nieliniowości prawej strony równania, a dokładniej, z przebiegiem y0⁢t, która jest funkcją nieciągłą. Wszystkie używane przez nas metody (w postaci solverów typu ode45) zakładają tymczasem, że aproksymowane rozwiązania są funkcjami dostatecznie gładkimi, co oczywiście w naszym przypadku nie jest założeniem spełnionym. To typowa dla obliczeń naukowych sytuacja, gdy istniejące oprogramowanie (i zaimplementowane w nim metody numeryczne) zmuszamy do pracy w warunkach wykraczających poza ich założone zastosowanie.

Gdyby na nasz problem rzucić okiem praktyka, natychmiast zwróciłby nam uwagę na to, że w ,,trudnej” wersji zadania (dla małych τ) występują dwie skale czasowe: skala szybkozmienna, związana z zewnętrznymi impulsami i skala wolnozmienna, reakcji yN na wymuszenie. To zaś podpowiada nam, by w takiej sytuacji zastosować krok całkowania krótszy niż długość najkrótszej skali czasowej zadania.

Rozwiązaniem siłowym byłoby więc zmuszenie solvera do pracy z krokiem czasowym wyraźnie mniejszym niż τ. Jednak z drugiej strony oznacza to nieprawdopodobne zwiększenie kosztu obliczeniowego potrzebnego do wyznaczenia rozwiązania: dla T=10 i τ=10-3 musielibyśmy wykonać co najmniej 105 kroków czasowych (zapewne znacznie więcej).

I faktycznie, o ile uzyskanie (najprawdopodobniej) błędnego wykresu yN wymagało tylko kilku- kilkunastu sekund cierpliwości, to wygenerowanie (bliższego prawdy, miejmy nadzieję) wykresu yN kodem:

opts = odeset('RelTol', 1e-10, ...
	'AbsTol', 1e-10, ...
	'MaxStep', tau, ... % redukcja kroku całkowania do najmniejszej skali czasowej
	'InitialStep', tau/3);
[tcomp, Y] = ode45(@rhs, [0,T], Y0, opts);

ciągnie się niemalże w nieskończoność (sprawdź, ile herbat zdążysz sobie przyrządzić, gdy τ=10-4…).

Czy można byłoby więc jakoś temu zapobiec?

Musimy raz jeszcze spojrzeć na naturę naszych kłopotów. Wszystko zdaje się pochodzić stąd, że solvery, których używamy, dobrze działają, gdy rozwiązania są gładkie. Tymczasem nasze rozwiązania są niegładkie. No dobrze, ale przecież tym samym jest to nie tylko problem obliczeniowy, ale także problem czysto matematyczny: jak bowiem zdefiniować rozwiązania równania różniczkowego, w którym prawa strona jest… nieciągła?!

W naszym przypadku, ściślej rzecz biorąc, prawa strona jest kawałkami stała — ciągła z wyjątkiem punktów postaci Ti±τ — zatem, gdyby ograniczyć się do przedziałów postaci

tam moglibyśmy sensownie zdefiniować nasze równanie różniczkowe (za wartość początkową w danym przedziale biorąc wartość końcową z przedziału poprzedniego). Gdybyśmy więc postąpili identycznie w naszych obliczeniach, całkowicie uniknęlibyśmy w ten sposób kłopotów z nieciągłościami v0 i ich automatyczną lokalizacją! Ceną za to jest wbudowanie w program pętli po pododcinkach, jednak — jak sami za chwilę się przekonamy — będzie warto ją zapłacić!

opts = odeset('RelTol', RTOL, ...
	'AbsTol', ATOL, ...
	'InitialStep', tau/3);
tcomp = []; Y = []; Y0i = Y0;
for i = 1:length(Ti)-1
	opts = odeset(opts, 'MaxStep', (Ti(i+1)-Ti(i))*0.1);
	[tcompi, Yi] = ode45(@rhs, [Ti(i),Ti(i+1)-tau], Y0i, opts);
	tcomp = [tcomp; tcompi];
	Y = [Y; Yi];
	Y0i = Yi(end,:);

	opts = odeset(opts, 'MaxStep', tau*0.1);
	[tcompi, Yi] = ode45(@rhs, [Ti(i+1)-tau,Ti(i+1)], Y0i, opts);
	tcomp = [tcomp; tcompi];
	Y = [Y; Yi];
	Y0i = Yi(end,:);
end
I = Y(:,1); V = Y(:,N);
plot(tcomp,[I,V]);