Zagadnienia

3.1 Doprecyzowanie zadania
3.2 Transformacja zadania do wygodniejszej postaci
3.3 A może wystarczy po prostu… wziąć większy młotek?

3. Funkcja odwrotna

W analizie pewnego modelu dotyczącego transportu morfogenów w tkankach [16] istotną rolę odgrywa funkcja

$\Sigma(s)=\dfrac{2}{M-1}\,(1-(1-s)^{{M-1}})-s\,(1-s)^{{M-1}},\qquad s\in[0,1),$

a ściślej: funkcja odwrotna do niej, $\Sigma^{{-1}}$ . Naszym zadaniem jest naszkicowanie wykresu $\Sigma^{{-1}}$ dla różnych wartości $0\leq M<1$ i zbadanie, jak bardzo $\Sigma^{{-1}}(z)$ jest bliskie 1, gdy $z$ jest rzędu kilkudziesięciu.

Jasne jest, że wszelkie nasze komputerowe działania musimy zacząć od zdefiniowania funkcji $\Sigma$ . Zrobimy to w Octave:

function v = sigmafun(s,M)
v = (2/(M-1))*(1-(1-s).^(M-1)) - s.*(1-s).^(M-1);
end

(Dalsze rozważania będziemy prowadzić dla $M=1/2$ ). Jeśli chcieć tylko z grubsza naszkicować wykres odwrotnej, to zadanie robi się banalne: wykres $\Sigma^{{-1}}$ to zbiór punktów postaci

$(z,\Sigma^{{-1}}(z)),\quad z\in[0,\infty),$

a to jest przecież ,,to samo”, co zbiór

$(\Sigma(s),s),\quad s\in[0,1).$

Wobec tego wystarczy skorzystać z funkcji plot, ale ,,na opak”:

M = 0.5;
s = linspace(0,1);
plot(sigmafun(s,M), s, 'r-', 'linewidth', 2);

Rys. 3.1. Wykres funkcji odwrotnej do $\Sigma$ , uzyskany metodą ,,inżynierską” dla $M=1/2$ . Ze względu na bardzo szybki wzrost $\Sigma$ (i nieokreśloność $\Sigma(1)$ ), wykres dostajemy jedynie dla $z<26$ .

Z wykresu na rysunku 3.1 widzimy, że $\Sigma^{{-1}}$ bardzo szybko przybliża się do wartości granicznej równej 1, co sugerowałoby wybór innej skali dla wizualizacji: tym, co mogłoby nas naprawdę zainteresować, mogłaby być szybkość przybliżania się $\Sigma^{{-1}}$ do 1: wybralibyśmy więc wówczas, zamiast zwykłego plot, znacznie bardziej wyraziste semilogy(sigmafun(s,M), 1-s). Jednak dla ustalenia uwagi, w dalszym ciągu pozostaniemy przy wykresie w skali liniowej.

3.1. Doprecyzowanie zadania

Niepokojący jest prawy koniec wykresu: wszak $\Sigma(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 $\Sigma^{{-1}}(60)$ i narysować wykres $\Sigma^{{-1}}$ na odcinku $[0,60]$ .

Napiszemy zatem w Octave funkcję, która będzie wyznaczać wartości funkcji odwrotnej do $\Sigma$ dla dowolnego zadanego $z\in[0,\infty)$ . Idea jest następująca: z definicji funkcji odwrotnej, $\Sigma^{{-1}}(z)=S$ , gdzie $S$ spełnia $\Sigma(S)=z$ , a więc $S$ jest miejscem zerowym funkcji

$F_{z}(S)=\Sigma(S)-z.$

Aby dla zadanego $z$ znaleźć $S$ , wystarczy znaleźć miejsce zerowe skalarnej funkcji $F_{z}$ — 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 funkcją jednej zmiennej, z wartościami $z$ i $M$ takimi, jakie były w chwili jej definiowania. Samą funkcję Fz określiliśmy zaś jako funkcję lokalną dla sigmainv. Wreszcie, aby uniknąć konsternacji fzero dla $z=0$ (czym spowodowanej?) zwracamy zawsze $\Sigma^{{-1}}(0)=0$ , bez niepotrzebnych obliczeń.

Jest oczywiste, że $\Sigma(s)^{{-1}}$ jest monotoniczna i może przyjmować wartości z przedziału $[0,1)$ . Dlatego wiedząc, że dla $s=1$ funkcja $\Sigma(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 $F_{z}$ !

3.2. Transformacja zadania do wygodniejszej postaci

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 $\Sigma$ , sprawa staje się jasna: przecież dla $s\approx 1$ , $\Sigma$ 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):

Gdy zadanie jest za trudne, należy zmienić… zadanie!

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

$G_{z}(s)\equiv 1-s-\left(\dfrac{p-s}{p+z}\right)^{{p/2}},\text{ gdzie }p=2/(1-M).$

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

Wyznaczanie wartości $\Sigma^{{-1}}(z)$ przez znalezienie miejsca zerowego funkcji $G_{z}$ 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 $\Sigma^{{-1}}$ (rysunek 3.2) dostajemy trzy razy szybciej aniżeli w przypadku poprzedniego podejścia, przy liczbie wywołań funkcji $G_{z}$ sześciokrotnie mniejszej niż w przypadku $F_{z}$ .

Rys. 3.2. Wykres funkcji odwrotnej do $\Sigma$ , uzyskany metodą przez znalezienie miejsca zerowego funkcji $G_{z}$ .

3.3. A może wystarczy po prostu… wziąć większy młotek?

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 $x_{i}=10^{{p_{i}}}$ , gdzie $p_{i}$ 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 $\Sigma(s)\leq 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 $\Sigma^{{-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).

Ćwiczenie 3.1

Czasami funkcja może być na tyle złośliwa, że zamiast zmieniać jest postać — jak to zrobiliśmy powyżej — możemy spróbować skorzystać z mniej wyrafinowanej, a za to ogólniejszej metody rozwiązywania równania nieliniowego $F_{z}(S)=0$ . Taką bardzo prostą metodą, o minimalnych wymaganiach wobec $F_{z}$ , jest metoda bisekcji (zadowoli się funkcją ciągłą, zmieniającą znak). Zaimplementuj metodę bisekcji i zobacz, jak sprawdzi się w warunkach naszego problemu dla oryginalnej funkcji $F_{z}$ .

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Obliczenia naukowe