Niepokojący jest prawy koniec wykresu: wszak Σ⁢1 jest nieokreślone! Zatem doprecyzujmy nasze zadanie: w naszym konkretnym zastosowaniu, chcielibyśmy wiedzieć, ze względu na naturę badanego zjawiska, jak blisko 1 jest Σ-1⁢60 i narysować wykres Σ-1 na odcinku 0,60.

Napiszemy zatem w Octave funkcję, która będzie wyznaczać wartości funkcji odwrotnej do Σ dla dowolnego zadanego z∈0,∞. Idea jest następująca: z definicji funkcji odwrotnej, Σ-1⁢z=S, gdzie S spełnia Σ⁢S=z, a więc S jest miejscem zerowym funkcji

Aby dla zadanego z znaleźć S, wystarczy znaleźć miejsce zerowe skalarnej funkcji Fz — a to już powinno być łatwe, gdy skorzystamy z funkcji fzero. Powyższą ideę implementuje funkcja sigmainv z poniższego listingu.

function [S, fc] = sigmainv(z, M)
if z == 0
	S = 0;
	fc = 0;
else
	[S, FS, info,out] = fzero(@(s)Fz(s,z,M), [0,1-eps/2]);
	if(info ~= 1)
		warning(['Klopoty dla z=', num2str(z), '; info=', num2str(info)]);
	end
	fc = out.funcCount; % zliczamy liczbe wywolan Fz
end
end % sigmainv

function Z = Fz(S, z, M) % "z" jest parametrem
Z = sigmafun(S, M)-z;
end

Zanim przyjrzymy się meritum funkcji sigmainv, zwróćmy uwagę na kilka programistycznych szczegółów. Jak pamiętamy, fzero żąda funkcji jednego argumentu, dlatego, dla każdego zadanego z, konstruujemy taką funkcję ad hoc — na podstawie trzyargumentowej Fz(S,z,M) — korzystając z mechanizmu funkcji anonimowej: funkcja

f = @(s)Fz(s,z,M)

Jest oczywiste, że Σ⁢s-1 jest monotoniczna i może przyjmować wartości z przedziału 0,1. Dlatego wiedząc, że dla s=1 funkcja Σ⁢s jest nieokreślona, przedział lokalizujący miejsce zerowe ustaliliśmy na [0,1-eps/2], gdyż 1-eps/2 jest największą liczbą maszynową mniejszą od 1 (dlaczego?).

Pozornie ,,asekurancki” sposób zdefiniowania funkcji sigmainv, w której zawsze sprawdzamy wartość parametru info zwracaną przez fzero, jest w istocie bardzo ważną i sensowną decyzją, pozwalającą wstępnie zweryfikować wyniki. Wszak nie wiemy a priori, czy nasz solver sobie poradzi z Fz!

Chociaż wykres dostajemy bez trudu, stosując prostą metodę kontynuacji — biorąc za przybliżenie początkowe wartość z poprzedniego punktu wykresu (por. wykład z Matematyki Obliczeniowej II), to przychodzi nam nań czekać dość długo. Gdy ponownie spojrzymy na wykres Σ, sprawa staje się jasna: przecież dla s≈1, Σ ma piekielnie stromy wykres, a to znacznie utrudnia działanie solvera powodując, że praktycznie ogranicza się on do metody bisekcji. Aby uniknąć tej niedogodności i znacząco przyspieszyć rysowanie wykresu, skorzystamy ze starej, dobrej zasady prowadzenia obliczeń numerycznych (i nie tylko):

Tutaj po prostu zastąpimy ,,trudną” funkcję Fz inną, ,,łatwiejszą”, ale o identycznych miejscach zerowych. Prosty rachunek pokazuje nam, że s jest rozwiązaniem równania Σ⁢s-z=0 wtedy i tylko wtedy, gdy s jest miejscem zerowym funkcji

Różnica między Fz a Gz jest taka, że ta ostatnia jest prawie liniowa, więc rzeczywiście powinna być ,,łatwa” dla solvera równania nieliniowego. Przy okazji, Gz nie ma już osobliwości: jest określona nawet dla s≥1.

Wyznaczanie wartości Σ-1⁢z przez znalezienie miejsca zerowego funkcji Gz pokazuje szczegółowo kolejny listing.

function [S, fc] = sigmainv2(z,M)
if z == 0
	S = 0;
	fc = 0;
else
	[S, FS, info, out] = fzero(@(s)Gz(s,z,M), [0,1]);
	if(info ~= 1)
		warning(['Klopoty dla z=', num2str(z), '; info=', num2str(info)]);
	end
	fc = out.funcCount; % zliczamy liczbe wywolan Gz
end
end % sigmainv2

function Z = Gz(S, z, M) % "z" jest parametrem
p = 2/(1-M);
Z = 1 - S - ((p - S)./(p+z)).^(p/2);
end

Faktycznie, dzięki tej zmianie wykres Σ-1 (rysunek 3.2) dostajemy trzy razy szybciej aniżeli w przypadku poprzedniego podejścia, przy liczbie wywołań funkcji Gz sześciokrotnie mniejszej niż w przypadku Fz.

Na zakończenie warto być może zauważyć, że w dzisiejszych czasach obliczenia są na tyle tanie, że nie zawsze warto odwoływać się do aż tak wyrafinowanych metod. Gdyby naszym głównym zadaniem było narysowanie, wyłącznie w celach poglądowych, jednego wykresu na odcinku 0,60, to zapewne byłoby nam dość obojętne, czy wyznaczymy sto, sto tysięcy, czy nawet milion wartości sigmafun: nasz komputer najpewniej i tak jest wystarczająco szybki, by z tym sobie błyskawicznie poradzić. Komenda logspace(k,m,N) pozwala wygenerować zestaw N węzłów ,,równoodległych w skali logarytmicznej”, czyli postaci xi=10pi, gdzie pi tworzą zestaw N równoodległych punktów z przedziału k,m.

Generując więc zestaw punktów skupiających się wokół s=1, a następnie wybierając tylko te, dla których Σ⁢s≤60:

N = 10000;
s = 1-logspace(-13,0,N); % punkty zagęszczają się wykładniczo wokół 1
Z = sigmafun(s,M);
good = find(Z<=60); % odrzucamy te, które wykraczają poza zakres zmiennej "z"
plot(Z(good), s(good));

udałoby się nam w ten sposób (nawet jeszcze szybciej, niż poprzednimi metodami!) ,,dociągnąć” wykres Σ-1 do z=59.999, co jest chyba całkiem niezłym wynikiem, jak na metodę brutalnej siły…

(Ale i tutaj ważne jest wyczucie: gdyby nieopatrznie położyć s = 1-logspace(-300,0,N), dostalibyśmy znacznie gorszy wynik).