Weźmy pod uwagę prosty przykład. Nasz program wymaga wielokrotnych konwersji z układu współrzędnych kartezjańskich na biegunowe. Załóżmy, że konwersja jest dokonywana w pętli (na przykład dlatego, że każdy następny krok zależy od poprzedniego).
Matlab ma wbudowaną funkcję:
[THETA,RHO] = cart2pol(X,Y);
Do testowania wydajności przygotujemy sobie zbiór danych wejściowych:
v = randn(1000000,2)*1000;
Teraz dla każdej pary wartości (x, y) ze zmiennej v przeliczmy współrzędne:
for a = 1:length(v)
[t r] = cart2pol(v(a,1),v(a,2));
end
Wykonanie powyższej pętli zajmuje na moim komputerze 11.692266 sekund.
Przeanalizujmy co robi funkcja cart2pol:
function [th,r,z] = cart2pol(x,y,z)
th = atan2(y,x);
r = hypot(x,y);
end
O ile atan2 jest operacją matematyczną, to hypot jest kolejną funkcją, która zwraca pierwiastek sumy kwadratów dwóch argumentów.
Co się zatem stanie jeśli wykonamy pętle z przykładu wyżej, ale zamiast funkcji cart2pol wywołamy obliczenia bezpośrednio?
for a = 1:length(v)
t = atan2(v(a,2),v(a,1));
r = sqrt(v(a,1)^2+v(a,2)^2);
end
Efekt jest identyczny. Różnica polega na tym, że czas obliczenia został zredukowany do 0.329579 sekundy, co daje przyspieszenie ponad 36 razy…
Wniosek – kod w matlabie warto analizować i optymalizować!
Niedoświadczony użytkownik może szukać w Matlabie funkcji wycinającej kawałek ciągu znaków. Należy wtedy pamiętać, że ciąg znaków w Matlabie jest macierzą (wektorem) znaków. Aby wyciąć kawałek ciągu znaków odwołujemy się zwyczajnie do zakresu wektora jak w przykładzie poniżej:
string = 'To jest ciąg znaków.';
string2 = string(9:12);
Inexperienced user may by looking for the function cutting out a piece of string in MATLAB. You need to remember that a string in MATLAB is a matrix (vector) of characters. To cut a piece of string out you have to refer to the part of that vector as in the example below:
string = 'This is simple string.';
string2 = string(9:14);
O trzęsieniu ziemie na Haiti, które miało miejsce we wtorek 12 stycznia 2010 słyszeli chyba wszyscy. Był to silny i płytki wstrząs, więc został dobrze zarejestrowany przez stacje sejsmiczne na całym świecie. Tak się składa, że w Stanach Zjednoczonych, które leżą stosunkowo blisko od Haiti trwa największy pasywny eksperyment sejsmiczny w historii. Ponad 400 przenośnych stacji sejsmicznych jest ustawianych na planie gęstej siatki obejmującej “pas” przechodzący przez USA od południa do północy. Pas ten jest stopniowo przenoszony z zachodniego wybrzeża w kierunku wschodniego. W tej chwili stacje są ustawione w środkowej części kraju, co narysowałem na mapie poniżej:
Ponieważ stacje używane w tym eksperymencie są bardzo nowoczesne, szerokopasmowe i profesjonalnie ustawiane zapisy są wysokiej jakości – pozbawione szumów, artefaktów, trendów, stałych składowych itp.
Stosunkowo prostu udało mi się wygenerować wykres – hodograf tego trzęsienia:
Czas na osi pionowej jest zredukowany, widać pięknie wejścia fal P, PP, S itd.
Kolejnym, nieco bardziej skomplikowanym krokiem była animacja pokazująca przejście czoła fali P przez terytorium USA. Kropki na animacji mają kolor zależny od wychylenia składowej pionowej sejsmografu. Im bardziej zielona tym większe wychylenie do góry, im bardziej czerwona tym większe wychylenie do dołu. Na animacji widać przejście czoła fali P po którym widoczne są drgania wynikające ze skomplikowanej budowy ziemi. Gdyby animacja nie była widoczna poniżej, można ją zobaczyć tutaj:
Pięknie widać, prawda?
Ciekawy jestem, czy zagląda tu ktoś, kto docenia moje dzieła sejsmiczne :)?
Od ostatniego wpisu udało mi się usprawnić dwa aspekty modelowania. Dodałem tłumienie fal P i S zgodne z modelem PREM oraz przyspieszyłem obliczenia poprzez wykorzystanie wielowątkowości w środowisku matlab.
Wstępną (ze względu na brak czasu na razie jedyną) animację pokazującą rozchodzenie się fal z tłumieniem można zobaczyć poniżej (lub tutaj):
Szukając różnych materiałów dotyczących analizy danach sejsmicznych w środowisku matlab trafiłem na stronę internetową Mike’a Thorne’a – docenta w departamencie geologii i geofizyki uniwersytetu w Utah: http://web.utah.edu/thorne/. Można tam znaleźć między innymi różne animacje pokazujące propagację fal sejsmicznych. Jest między innymi animacja pokazująca model PREM. Skromnie jednak uważam, że moja jest… ciekawsza 🙂
Pracując ostatnio nad symulacją rozchodzenia się fal sejsmicznych natknąłem się na problem związany z wygenerowaniem animacji w środowisku Matlab. We wspominanym przed chwilą opisie można zobaczyć efekt jaki udało mi się uzyskać, czyli animację w wysokiej rozdzielczości.
Matlab oferuje dwie standardowe drogi wygenerowania animacji z serii figur (zawierających np. wykres). Jedna droga to stworzenie obiektu avi i dodawanie do niego kolejnych klatek animacji. Druga metoda to stworzenie zmiennej (struktury) z kolejnymi klatkami i zapisanie jej do pliku avi. Metody różnią się wykorzystaniem pamięci programu i zapewne wydajnością (coś za coś).
Niestety obie metody mają wadę, której nie dało mi się w nich wyeliminować: rozdzielczość i jakość filmu jest marna a na dodatek jego rozmiar ze względu na kiepską kompresję jest ogromny.
Mój pomysł na wygenerowanie animacji jest odmienny. Każdą klatkę filmu zapisujemy do pliku graficznego, nadając plikom kolejne nazwy. Wystarczy użycie funkcji print(). Dla przyspieszenia procesu można nie wyświetlać wykresów na ekranie i rysować w pamięci. W poniższym przykładzie rysujemy 50 wykresów z losowymi punktami i każdy wykres zapisujemy do pliku bez wyświetlania:
figure('Visible','Off') % Otwieramy "ukryte" okno
for i=1:50 % 50 powtórzeń
plot(rand(1,100),rand(1,100),'.'); % Rysujemy 100 punktów
print(gcf, '-r120', '-dpng', [num2str(i) '.png']); % Zapisujemy
% wykres do pliku .png
clf; % Czyścimi okno
end
close all; % Zamykamy okno
Szczegóły funkcji print() można znaleźć w dokumentacji.
Oczywiście to nie koniec przygody. Czas przerobić klatki w film. Bardzo prosto można to zrobić programem VirtualDub. Uruchamiamy program i przeciągamy pierwszą klatkę (ważne są kolejne numery klatek) do głównego okna programu.
W menu “video” możemy ustawić ilość klatek na sekundę “Frame Rate” i kompresję oraz jej parametry “Compression”. Następnie z menu file wybieramy “Save as AVI”. Otrzymujemy gotową animację. Możemy generować obrazu różnej wielkości (rozdzielczości) a co za tym idzie rozdzielczość animacji. Różnica w jakości wyjściowej i potrzebnej pamięci (zarówno w RAM i na dysku twardym) jest ogromna, więc polecam tą metodę.
Kroki w programie Virtual Dub opisałem pobieżnie, ale mam nadzieję, że każdy sobie poradzi.
Niektórych dziwi fakt, że na wydziale Fizyki Uniwersytetu Warszawskiego można się zajmować geofizyką, bo jest ona bardziej geo niż fizyką. Okazuje się, że są dwa rodzaje geofizyki: geofizyka i geofizyka. Postanowiłem skierować swoją uwagę w kierunku tej drugiej. Ponieważ z wielu tematów obejmowanych przez geofizykę najbardziej zainteresowała mnie sejsmologia to własnie tym tematem chcę się dalej zajmować.
Jestem fanem poglądu, że z samej teorii to człowiek się niczego nie nauczy postanowiłem wspierać swoją edukację elementami praktycznymi. Niestety, do wykonania doświadczenia związanego z sejsmologią potrzeba niemałych nakładów, więc pierwsze kroki skierowałem do modelowania komputerowego. Podobno można w tej dziedzinie jeszcze dużo zdziałać, więc mam przynajmniej kilka powodów dla których warto.
Ponieważ od około roku jestem miłośnikiem środowiska Matlab, wykorzystałem je do stworzenia swojego modelu.
Krótki wstęp teoretyczny dla (prawie) niefizyków
Model pokazuje kierunki rozchodzenia się fal sejsmicznych w ziemi i nie zajmuje się takimi zagadnieniami jak częstotliwości fal.
Fale sejsmiczne nie rozchodzą się po prostych, ponieważ w środku Ziemi prędkości fal nie są stałe i w modelu jednowymiarowym zależą wyłącznie od głębokości. Fala pokonuje drogę między dwoma punktami taką trasą, aby czas pokonania tej drogi był jak najkrótszy (dla zainteresowanych Prawo Snelliusa, Zasada Fermata).
W miejscach, gdzie prędkości fal zmieniają się wraz z głębokością płynnie sprawa jest dość prosta – promienie (czyli kierunki fal) po prostu sobie skręcają.
W miejscach gdzie gdzie istnieje ostra granica (czyli skokowa zmiana prędkości) dochodzi dodatkowo do zjawiska konwersji fal. W tym miejscu należy wspomnieć, że istnieją różne rodzaje fal sejsmicznych. We wnętrzu Ziemi mamy do czynienia z falami podłużnymi (oznaczane literką P) i poprzecznymi (oznaczane literką S) – dla zainteresowany: Fala sejsmiczna. Fale P i S rozchodzą się z różnymi prędkościami. W momencie gdy na granicę pada fala (P lub S) powstają w danym miejscu 4 fale: odbita P, odbita S, załamana P i załamana S. Z rożnych względów (całkowite wewnętrzne odbicie, zerowa prędkość fal S) czasami powstaje mniej niż 4 fale po konwersji. Dodatkowo energia i amplituda fal po konwersji jest rożna dla każdej z nich, co zależy od własności granicy.
Model
W dopiero co stworzonej pierwszej wersji modelu rozpatruję przekrój przez środek Ziemi czy badam tory promieni w dwóch wymiarach.
Budowę Ziemi (a dokładnie rozkład prędkości) przyjąłem jako jednowymiarową (czyli prędkości fal P i S oraz gęstość zależą wyłącznie od głębokości) według modelu Preliminary Reference Earth Model (PREM) (Dziewonski & Anderson, 1981).
Przyjąłem istnienie 6 granic, na których zachodzi konwersja fal (6368 km, 6346.6 km, 6151 km, 5701 km, 3480 km, 1221.5 km). Dla każdej granicy obliczyłem współczynniki konwersji fal (czyli amplitudy fal powstały w stosunku do amplitudy fali padającej w zależności od parametrów granicy i kąta padania).
Jeden z wyników
Ciekawych wyników z takiego modelu można uzyskać bardzo dużo. Tutaj chcę pokazać jeden z pierwszych (gdyż liczenie trwa bardzo długo) rezultatów – animację.
Realizując opisany model poprzez wygenerowanie w źródle (na głębokości 241 km) 360 promieni fal P i numeryczną symulację ich rozwoju (czyli torów, odbić i załamań) udało mi się stworzyć bardzo ładną animację pokazującą rozchodzenie się fal we wnętrzu Ziemi. Gdyby film poniżej nie działał to można go zobaczyć tutaj. Ponieważ przeciętnie udała mi się legenda wyjaśniam, że kolorem zielonym zaznaczone są fale P a czerwonym fale S.
Moim skromnym zdaniem efekt jest bardzo atrakcyjny :).
Tyle wstępu do problemu – jak dostanę kolejne ciekawe rezultaty na pewno się podzielę. Mam w planie podzielić się również spostrzeżeniami dotyczącymi Matlaba, gdyż nauczyłem się przez ostatnie kilkanaście dni wielu ciekawych rzeczy.
