Strona domowa ćwiczeń z Metody Numerycznych II A
w roku 2000/2001

Metody Numeryczne II A

WYNIKI z Metod Numerycznych II A (w formacie HTML): wyniki testu, punkty za obecność i programy (zadania) oraz proponowane oceny. Wyniki zostały uaktualnione (patrz przedostatnia kolumna). Skończył się już termin oddawania programów w sesji letniej!

Wydruk wyników wisi na tablicy ogłoszeń przy pracowni komputerowej, oraz przy pokoju 125 (te wyniki nie są aktualne). Można na nich zobaczyć, ile punktów otrzymało się za test, zadanie programistyczne na teście, ile się ma obecności i ile punktów za programy oddane przed początkiem sesji.

Wyniki wraz z propozycjami ocen wywieszone są na drzwiach pokoju wykładowcy, dr hab. Tomasza Wernera (pokój 115). Wszystkie uaktualnienia zostały wysłane do dr Wernera za pomocą e-mail; w razie gdy wywieszone oceny nie zostały uaktualnione można prosić o sprawdzenie poczty.
Patrz także progi ocen.

Zasady zaliczania:

• 50% - zadania na ćwiczeniach (programy)
• 35% - test z wiadomości z wykładu
• 15% - obecność na ćwiczeniach

Punktacja testu

• Jedno pytanie testowe - 2,5 pkt
• Program do napisania - do 10 pkt

Osoby które nie były na teście lub uzyskały z niego słabe wyniki, i nie wystarczy im do zaliczenia oddanie zaległych programów, mogą zdawać test u dr Wernera. Nie będzie ogólnego testu poprawkowego.

Punktacja zadań w grupach 2 i 3:

• Oddane zadanie - 1 pkt.
• Zadanie oddane w wersji rozszerzonej - 1 1/3 pkt.
• Zadanie dodatkowe - +1 pkt.
• Dodatkowe rozszerzenie - od +1/3 pkt.
• Szczególnie dobre rozwiązanie - od +0,2 do +1 pkt.

Z obowiązkowych zadań można zdobyć 11 1/3 pkt., które będą przeliczone na 50 pkt (co daje około 4,41 pkt. za program).

Ostatecznie progi ocen oraz zaliczenia  
 

Programy można wysyłać pocztą (mój adres e-mail to jnareb@fuw.edu.pl), ale preferowane jest oddawanie programów osobiście. Przynajmniej połowa programów powinna być oddana osobiście (czy to na ćwiczeniach, czy też umówiwszy się ze mną na konkretny termin). Oceniane jest działanie programu oraz sposób jego zapisania (t.j. ocena bazuje także na źródłach programu).
Po sprawdzeniu programów za pomocą e-mail wysyłam komentarze na adres zwrotny, tak więc proszę sprawdzać pocztę!. W szczególności wywieszone listy (oprócz listy wywieszonej na drzwiach pokoju prof. Wernera) nie są aktualizowane po sprawdzeniu każdego programu.

W poradach znajduje się opis preferowanego sposobu (stylu) kodowania (zapisywania programu). Nie jest on obowiązkowy. Służy on głównie jako wskazówka, jak należy pisać programy. Za przestrzeganie tego można będzie dostać dodatkowe punkty.  
 

Zagadnienia przedstawiane na ćwiczeniach

1. Wartości własne i wektory własne macierzy symetrycznych

1.1 Opis teoretyczny

• W eigen.ps opisana jest metoda potęgowa znajdowania wartości własnej o największym module i odpowiadającego jej wektora własnego, odwrotna metoda potęgowa (metoda iteracji odwrotnej) znajdowania wartości własnej najbliższej zadanej liczbie \lambda^* i odpowiadającego jej wektora własnego, oraz metoda Gaussa rozkładu macierzy na iloczyn macierzy trójkątnej dolnej i trójkątnej górnej LU (potrzebny do rozwiązywania układu równań w odwrotnej metodzie potęgowej). Przykładowa implementacja metody Gaussa (z częściowym wyborem elementu podstawowego) oraz metoda rozwiązywania układu równań z macierzą o zadanym rozkładzie LU znajduje się w pliku gauss.c. Uwaga: wersja ta używa tablic alokowanych dynamicznie (o zmiennej wielkości); tablice te należy zaalokować a następnie zwolnić.
• W eigen1.ps opisane są metody znajdowania wszystkich wartości własnych (oraz wektorów własnych) za pomocą transformacji podobieństwa. Przedstawiona jest metoda Householdera redukcji macierzy symetrycznej do macierzy trójdiagonalnej (przykładowa nieprzetestowana implementacja znajduje się w pliku tridiag_reduce.c) oraz metoda wyznaczania wartości własnych macierzy trójdiagonalnej (redukcji macierzy trójdiagonalnej do diagonalnej) za pomocą rozkładu LR.

1.2 Programy

• gauss.c - redukcja macierzy do postaci LU metodą Gaussa z częściowym wyborem elementu podstawowego; rozwiązywanie układu równań z macierzą zadaną przy pomocy rozkładu LU. Dostosowany do dynamicznej alokacji pamięci; przy stosowaniu macierzy o stałej wielkości wymaga drobnych zmian. Powinien działać poprawnie, ale tego nie gwarantuję!
• gauss.h - plik nagłówkowy do powyższego modułu.
• tridiag_reduce.c - redukcja macierzy symetrycznej do macierzy trójdiagonalnej za pomocą transformacji Householdera. Nie przetestowane! Wkrótce powinna pojawić się nowa wersja, bazująca na bibliotece lapack (Linear Algebra PACKage).

1.3 Zadania

1. Napisać program znajdujący największą wartość własną macierzy i odpowiadający jej wektor własny za pomocą metody potęgowej
2. Napisać program znajdujący wartość własną macierzy najbliższą zadanej wielkości i odpowiadający jej wektor własny za pomocą metody iteracji odwrotnej (odwrotnej metody potęgowej). Można stosować metodę eliminacji Gaussa bez wyboru elementu podstawowego.

2. Równania różniczkowe zwyczajne

2.1 Opis teoretyczny

• W pliku rozn.ps znajduje się krótkie wprowadzenie do równań różniczkowych zwyczajnych i wstęp do metod różnicowych (numerycznych metod rozwiązywania równań różniczkowych w określonych punktach). Przedstawiona (i opisana) jest metoda Eulera i jej modyfikacje: ulepszona metoda Eulera (metoda Heuna) oraz zmodyfikowana metoda Eulera. Opisana jest ogólna postać metody Runge-Kutty i podana jest klasyczna metoda Runge-Kutty. Opisane są także liniowe metody wielokrokowe oraz metody typu predyktor-korektor.
• W pliku ode.ps.gz przedstawiona została klasyczna metoda Runge-Kutty czwartego rzędu oraz ulepszona metoda Eulera (metoda Runge-Kutty drugiego rzędu). We wstępie opisana została notacja; znajduje się tam także ogólna informacja o sposobie konstrukcji metod różnicowych (na podstawie kwadratur). Sposób "wyprowadzenia" metody Runge--Kutty podaję za wykładem Ryszarda Kutnera "Symulacje w Materii Skondensowanej".
  

2.2 Zadania

1. Ruch w polu siły oporu lepkiego: v'(t) = - k v(t); należy obliczyć prędkość cząstki po zadanym czasie. Uzupełnienie: zadanie należy rozwiązać za pomocą metody Eulera, ulepszonej metody Eulera (metody Heuna), zmodyfikowanej metody Eulera (metody Runge-Kutty 2. rzędu) oraz metody Runge-Kutty 4. rzędu. Należy zrobić wykres rozwiązania od czasu i błędu (różnicy rozwiązania numerycznego i analitycznego); pod Linuksem zapisując wyniki do pliku i korzystając z programu gnuplot, pod DOS-em tworząc własną procedurę do robienia wykresu lub zapisując wyniki do pliku i korzystając z jakiegoś programu.
2. Wahadło matematyczne; siła dana jest wzorem F(x) = - k sin(x); należy rozwiązać zagadnienie metodą Runge-Kutty dla kilku prędkości początkowych. Wyniki powinny być przedstawiane za pomocą wykresu fazowego (prędkość od położenia); można to zrobić za pomocą zapisywania wyników do pliku i użycia programu gnuplot. Należy normalizować położenie (kąt) do przedziału [-pi,pi].
   Uzupełnienie: należy także rozwiązać za pomocą tego samego (bądź zmodyfikowanego) programu równanie oscylatora harmonicznego. Użyć metody Eulera oraz metody Runge-Kutty 4. rzędu. Porównać obie metody, korzystając z analitycznego (ścisłego) rozwiązania. W tym przypadku nie jest potrzebna normalizacja położenia do przedziału [-pi,pi] (i nie należy tego robić). Do robienia wykresu fazowego "biegnącego" można się posłużyć skryptem w bashu wykresik.sh. Można się z niego nauczyć posługiwania gnuplotem.
3. Metody wielokrokowe jawne: zbadać zbieżność i koszt (ilość obliczeń funkcji prawej strony równania) dla zagadnienia o znanym rozwiązaniu, np. oscylatora harmonicznego bądź równania u' = a u, gdzie a > 0. Uzupełnienie proszę także rozwiązać zagadnienie:
   u' = ((2 log u + 8)/t -5) u z warunkiem początkowym u(1) = 1 dla przedziału t = [1, 6]. Rozwiązanie ścisłe to u(t) = exp(-t² + 5 t - 4). Należy zrobić wykres rozwiązania od czasu (czy to bezpośrednio w programie, czy za pomocą programu gnuplot). Porównać dokładność i ilość obliczeń funkcji z klasyczną metodą Runge-Kutty. Porównać błąd dla różnych wielkości kroku i ewentualnie wyliczyć rząd metody (podpowiedź: log(xⁿ) = n log(x)). Uwaga: metoda Milne'a nie jest stabilna i jej rząd wyliczany na podstawie równania oscylatora harmonicznego może być niższy niż 4, a nawet może dawać błędne wyniki.
4. Dla chętnych: metody typu predyktor-korektor. Porównać właściwości z metodami wielokrokowymi, stosując je do tych samych zagadnień

3. Dyskretna transformata Fouriera i jej zastosowania

3.1 Opis teoretyczny

• W pliku four.ps znajduje się wstęp do transformaty Fouriera i opis dyskretnej transformaty Fouriera. Opisane są szeregi Fouriera, przejście do szeregów w zmiennych zespolonych, a następnie do transformaty Fouriera. Następnie opisana została dyskretna transformata Fouriera, którą potrafimy obliczać na komputerach (dla dowolnych funkcji/danych).
• W pliku fft.ps.gz (tex) znajduje się opis algorytmu szybkiej transformaty Fouriera. Przedstawione są tam również własności transformaty Fouriera i jej dyskretnej wersji. Podane są także najprostsze zastosowania transformaty Fouriera.
• W pliku fourrsc.ps znajduje się opis obliczania rzeczywistej transformaty Fouriera, transformaty sinusowej i cosinusowej.
• FOURIER TRANSFORM - opis z Numerical Recipies (kopia)
  1. Introduction "INTFTNR.ps"
  2. Fourier Transform of Discretely Sampled Data "FTNR.ps"
  3. Fast Fourier Transform (FFT) "FFTNR.ps"

3.2 Programy

• Zestaw funkcji z Numerical Recipies do obliczania szybkiej transformaty Fouriera jako archiwum tar: fft_NRC.tar. Patrz także "Gotowe programy i skrypty" na tej samej stronie. Opis stosowania tych metod znajdziesz tamże oraz w fft.ps.gz (tex) (funkcji four1) i w fourrsc.ps (pozostałych, w szczególności realft).

3.3 Zadania

1. Napisać program obliczający zespoloną dyskretną transformatę Fouriera z jej definicji:
   
   H_n = sum(k=0..N-1) h_k e^{2 piikn/N}
   gdzie i oznacza jednostkę urojoną (t.j. ze wzoru nr 15 w four.ps). Sprawdzić, czy odwrotna transformata Fouriera, zadana wzorem:
   
   h_k = 1/N sum(n=0..N-1) H_n e^{-2 piikn/N}
   przeprowadza otrzymaną tablicę H z powrotem w oryginalną funkcję. Można przykładowo wziąć funkcję zadaną wzorem (wypełnić tablicę h w następujący sposób):
   
   • h_k = sin(2 pi k/N) + i0
   • h_k = exp(i8*2 pi k/N)
   • h_k = (k < N/2 ? 1 : i)
   Zrobić wykres oryginalnej funkcji, jej transformaty Fouriera oraz odwrotnej transformaty Fouriera transformaty Fouriera. Narysować część rzeczywistą i urojoną lub moduł (wartość bezwględną). Pod Linuksem można to zrobić za pomocą programu gnuplot poleceniem (zakładając, że do pliku h.dat zostały zapisane kolejne elementy tablicy h, w każdym wierszu część rzeczywista i część urojona (oddzielone np. spacją)):
    plot 'h.dat' using 1 title "re(h)" with linespoints
    replot 'h.dat' using 2 title "im(h)" with linespoints
    replot 'h.dat' using (abs({$1,$2})) title "abs(h)" with linespoints
   Polecenie replot powoduje wykreślenie dodatkowych danych na bieżącym wykresie. Można pominąć fragment with linespoints: służy on do określenia sposobu wykreślania danych (tutaj: punkty połączone liniami). {a,b} oznacza w gnuplocie liczbę zespoloną a + i b. Argument using z następującym po nim wyrażeniem określa jakie dane będą wykreślone na wykresie.
2. Napisać program sprawdzający własności transformaty Fouriera. Sprawdzić jakie własności ma transformata Fouriera funkcji
   1. rzeczywistej, parzystej
   2. rzeczywistej, nieparzystej
   3. urojonej, parzystej
   4. urojonej, nieparzystej
   gdzie nieparzystość oznacza, że h_k = h_{N - k}.
   Uzupełnienie: sprawdzić własność tzw. aliasingu, t.j. wykreślić funkcję o częstości większej od krytycznej częstości Nyquista 2/N, jej wartości w punktach próbkowania, jej dyskretną transformatę Fouriera oraz transformatę odwrotną. Przykładową funkcją może być h(t) = sin(2*N*2 pi t), gdzie t_k = k/N, zaś h_k = h(t_k).
3. Zrobić prosty filtr przy użyciu dyskretnej transformaty Fouriera, tzn. obliczyć transformatę Fouriera, wyciąć odpowiednie częstości (pamietając o tym, że dla zwykłej, zespolonej dyskretnej transformaty Fouriera po indeksie N/2 następują częstości ujemne). Uzupełnienie: zrobić to także przy użyciu rzeczywistej transformaty Fouriera, wylicznej za pomocą własnej procedury (na podstawie fourrsc.ps) lub procedury realft z NRC (zawartej w pliku fft_NRC.tar). Uwaga: procedura realft w tym archiwum została poprawiona. Poprzednia wersja wymagała poprawienia deklaracji procedury four1 i podawania połowy wielkości tablicy jako parametru.
4. Dla chętnych: mnożenie wielomianów przy użyciu [szybkiej] transformaty Fouriera. Za pomocą transformaty Fouriera przechodzimy od reprezentacji wielomianu P(t) za pomocą współczynników a₀ + a₁ x + ... do reprezentacji za pomocą wartości w deg(P) + 1 punktach P(z_i); w naszym przypadku (stosowania TF) są to zespolone pierwiastki z jedności. Wymnożyć wielomiany punkt po punkcie. Obliczyć współczynniki iloczynu z wartości w punktach za pomocą odwrotnej transformaty Fouriera. Uwaga: należy pamiętać, że do odtworzenia wielomianu potrzeba wartości w przynajmniej deg + 1 punktach. Należy wypisać oba czynniki oraz iloczyn (wynik) w czytelnej formie; bardzo małe współczynniki można przyjąć równe zeru przy wypisywaniu.

4. Splatanie i rozplatanie przy użyciu transformaty Fouriera

4.1 Opis teoretyczny

• W pliku fft-splot.ps.gz (tex) znajduje się opis splatanie i rozplatania za pomocą transformaty Fouriera i iteracyjne rozplatanie, oraz filtr optymalny Wienera. Podana jest definicja splotu dla funkcji ciągłych oraz dyskretnego, podane są twierdzenia wiążące transformatę Fouriera splotu (zarówno w przypadku ciągłym jak i dyskretnym) z transformatami czynników. Opisany jest sposób eliminacji efektów brzegowych przy liczeniu splotu (i odplataniu) za pomocą transformaty Fouriera przez dopełnienie sygnału i funkcji odpowiedzi jednostkowej zerami. Opisana jest operacja rozplatania przy użyciu transformaty Fouriera. Podany jest sposób odplatania funkcji odpowiedzi przy obecności szumu: filtr optymalny Wienera oraz iteracyjna metoda odplatania.
• Splatanie przy pomocy FFT i filtr Wienera - opis z Numerical Recipies (kopia)
  1. Convolution and Deconvolution Using the FFT "splotnr.ps"
  2. Optimal (Wiener) Filtering with the FFT "wiennr.ps"

4.2 Zadania

1. Przetestować twierdzenie o splocie dyskretnym. Obliczyć splot sygnału (np. funkcji okresowej piłokształtnej) z~funkcją odpowiedzi impulsowej o~małej szerokości (np. funkcji Gaussa lub prostokątnej bariery) z definicji splotu:
   (r * s)_j = sum(k = -M/2+1..M/2) r_k s_j-k, oraz za pomocą transformaty Fouriera z twierdzenia o splocie dyskretnym:
   (R * S)_j = R_j S_j. Należy pamiętać o dopełnieniu sygnału zerami (jeśli sygnał traktujemy jako zwykłą funkcję) lub 'zawijaniu' indeksów (jeśli sygnał jest funkcją okresową) i ułożeniu funkcji odpowiedzi w 'wrap-around order' (r_-k = r_N-k). Odpleść funkcję odpowiedzi ze splotu obliczonego z definicji i zobaczyć, jak dokładnie odtwarzany jest sygnał oryginalny. Odpleść sygnał za pomocą wzoru iteracyjnego:
   q -> q + (c - r * q), gdzie c jest wynikiem splotu, q wynikiem w kolejnej iteracji, zaś r funkcją odpowiedzi jednostkowej. Dla chętnych: sprawdzić wpływ szumu na dokładność odtwarzania.
2. Dla chętnych: przetestować filtr Wienera dany wzorem
   F(f) = (|H(f)|²)/(|H(f)|² + |N(f)|²), (gdzie h = r*s, n oznacza szum) i porównać wynik jego stosowania:
   S'(f) = (C(f) F(f))/R(f), z iteracyjnym odplataniem (dla sygnału splatanego z funkcją odpowiedzi w obecności szumu c = (r * s) + n).

5. Rozwiązywanie równań różniczkowych cząstkowych metodą różnic skończonych (MRS)

5.1 Opis teoretyczny

• W pliku czastk.ps przedstawiony jest problem numerycznego rozwiązywania równań różniczkowych cząstkowych. We wstępie opisane jest przykładowe równanie: równanie dyfuzji. Następnie przedstawione są założenia metody różnic skończonych i różne metody przybliżania pochodnej za pomocą ilorazu różnicowego. Przedstawione zostały różne schematy rozwiązywania równań różniczkowych cząstkowych: FCTS (w tym sposób wyboru stosunku kroku czasowego do przestrzennego), BCTS, DuFort-Frankel, Crank-Nicholson oraz metoda Richardsona poprawiania dokładności schematu.

5.2 Zadania

1. Rozwiązać numerycznie równanie dyfuzji: D_t T(x,t) = D_xx T(x,t) (dla współczynnika dyfuzji D = 1), dla 0 <= x <= 1, t >= 0, z warunkami początkowo-brzegowymi
   • T(x,0) = sin(pi x/2) + 1/2 sin(2 pi x)
   • T(0,t) = 0
   • T(1,t) = exp(-pi^2/4 t)
   Równanie to ma rozwiązanie T(x,t) = exp(-pi^2/4 t) sin(pi x/2) + 1/2 exp(-4 pi^2 t) sin(2 pi x) Rozwiązać to zagadnienie przynajmniej dwoma metodami jawnymi. Dla chętnych: przetestować inne podane na wykładzie metody różnicowe.

Zadania dodatkowe (na ocenę celującą)

Tutaj będą podawane propozycje zadań (programów) zaliczeniowych na ocenę celującą (5+). Można zgłaszać także własne propozycje, po uzgodnieniu z prowadzącym ćwiczenia.

• zadanie-planety.ps.gz (tex) - symulacja ruchu planet (rozwiązywanie równań ruchu w polu grawitacyjnym). Zagadnienie: równania różniczkowe zwyczajne.
• Rozwiązywanie równania Schroedingera zależnego od czasu kilkoma metodami. Zagadnienia: równania rózniczkowe zwyczajne, równania rózniczkowe cząstkowe, [szybka] transformata Fouriera,...

Pojawiło się już pierwsze rozwiązanie z zadań dodatkowych. Program został napisany przez Władysława Saja (grupa 1.), na system operacyjny Linux. Program rozwiązuje zagadnienie ruchu satelity w polu grawitacyjnym Ziemi i Księżyca (zagadnienie trzech ciał). Wykorzystano metodę Runge-Kutty ze stałym i ze zmiennym krokiem. Program dostępny jest w katalogu ~wmsaj/grav3/ (czyli fs98/wmsaj/grav3/). Na stronie Pawła Jankowskiego znajduje się krótki opis tego programu (w PostScripcie).
Uwaga: zagadnienie trzech ciał nie jest już przyjmowane jako rozwiązanie zadania dodatkowego!

Gotowe programy i skrypty

• gauss.c - redukcja macierzy do postaci LU metodą Gaussa z częściowym wyborem elementu podstawowego; rozwiązywanie układu rownań z macierzą zadaną przy pomocy rozkładu LU. Dostosowany do dynamicznej alokacji pamięci; przy stosowaniu macierzy o stałej wielkości wymaga drobnych zmian. Zobacz także informacje w poradach dotyczących programowania rozdział dotyczący przekazywania tablic do funkcji.
• gauss.h - plik nagłówkowy do gauss.c.
• tridiag_reduce.c - redukcja macierzy symetrycznej do macierzy trójdiagonalnej za pomoca transformacji Householdera. Nie przetestowane!
  
  
• dmetodyn.c - różne metody różnicowe rozwiązywania układu równań różniczkowych zwyczajnych.
• dmetodyn.h, dtypyn.h, real.h - pliki nagłówkowe do dmetodyn.c
• rk4.c - klasyczna metoda Runge-Kutty 4. rzędu z "Numerical Receipies in C"
  
  
• Zestaw funkcji z Numerical Recipies do obliczania szybkiej transformaty Fouriera: jako archiwum tar: fft_NRC.tar (które pod Linuksem rozpakowywuje się za pomocą polecenia tar xvf fft_NRC.tar, pod Windows np. WinZip czy WinRAR, a także Windows Commander). Pliki w archiwum są dołączane w swoim formacie (w tym przypadku jako tekst), więc w razie potrzeby mozna otworzyć to archiwum w edytorze i wyciągnąć z niego zawartości odpowiednich plików. Pliki przykładowe x*.c (oprócz xfour1.c) wymagają zamiany #include "nr.h" na #include "nrcfft.h". Znajdujący się w archiwum plik Makefile (korzystający z rozszerzeń GNU make) służy do kompilacji programu przykładowego xfour1.c, testującego właściwości dyskretnej [zespolonej] transformaty Fouriera. Pod Linuksem możemy skompilować ten program za pomocą polecenia make -k (gdzie opcja -k oznacza, że program make ma się nie zatrzymywac na błędach kompilacji), zaś pod Emacsem za pomocą polecenia M-x compile.
  1. four1.c - FFT w jednym wymiarze
  2. twofft.c - FFT dwu funkcji rzeczywistych na raz
  3. realft.c - FFT funkcji rzeczywistej
  4. sinft.c - transformata sinusowa
  5. cosft.c - transformata cosinusowa
  Uwaga: funkcje z NRC wymagają tablic o elementach od 1 do N (mozna to obejść podając jako parametr tablica-1), a rozmiar tablicy musi być potęgą dwójki. Procedura four1 jako parametr dostaje ilość liczb zespolonych, a nie rozmiar tablicy liczony w liczbach rzeczywistych.

Lista obecności i lista oddanych zadań

Stan na dzień 01.06.2001.

Grupa 2. (czwartkowa): lista2 (lista oddanych zadań w grupie 2. w formacie HTML).
Grupa 3. (piatkowa): lista3 (lista oddanych zadań w grupie 3. w formacie HTML).

Patrz także: Strona Pawła Jankowskiego do Metod Numerycznych II A.
Zadania oraz notatki do ćwiczeń z Metod Numerycznych II A znajdują się także w bibliotece IFD.

Powrót na górę strony