Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Obliczenia naukowe – 3. Funkcja odwrotna – MIM UW

Zagadnienia

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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.