ONA 7 – Wartości i wektory własne

Dziś poznaliśmy numeryczne rozwiązania problemu znajdowania dominujacych wartości własnych metodą potęgową i wartości własnych bliskich zadanej wartości metodą odwrotną potęgową. Dużo więcej materiału można znaleźć na tej stronie

Jeśli chodzi o funkcje przydatne w pythonie, to przede wszystkim interesują nas funkcje:

scipy.linalg.eig(A)

oraz

scipy.linalg.eigvals(A)

Zadania na dziś:

1. Zaimplementuj metodę potęgową z wykładu i zastosuj ją dla rozważanej na wykładzie macierzy: M=array([[1,2.],[4,3]]). Czy Twoje rozwiązanie różni się od rozwiązania z funkcji eig? Jakie równanie spełniają wektory własne zwrócone przez tę funkcję?

2. Zaimplementuj metodę Rayleigha RQI (z wykładu ). zastosuje ją do znalezienia najmniejszej wartości własnej macierzy M z zad. 1.

3. Rozważ macierz A=array([[-1,2,2],[2,2,-1.],[2,-1,2]]). Spróbuj znaleźć jej wartości własne zaimplementowanymi przez siebie metodami. Skąd biorą się problemy? Czy podobne problemy spotkają nas dla macierzy rand(3,3)?

4. Rozważ macierz M=array([[a,1.],[0,b]]). Gdzie a i b są dodatnimi liczbami całkowitymi. Jak zmienią się wartości własne, gdy zamiast zera wstawimy do M[1,0] wartośći 1/k dla k rosnącego wykładniczo?

 

ONA 6 – Interpolacja funkcji

Dzisiaj zajmiemy się Interpolacją funkcji przy pomocy wielomianów i funkcji sklejanych.

Teoria z wykłądu jeśli chodzi o wielomiany znajduje się tu a jeśli chodzi o funkcje sklejane tu.

Większość interesujących nas dzisiaj funkcji znajdziemy w module scipy.interpolate, ale najprostsze funkcje polyfit  i poly1d znajdują się w module numpy.  W module interpolate interesują nas funkcje:

Na laboratorium będziemy rozwiązywać następujące problemy:

1. Spróbuj zinterpolować funkcję sinus(x) na przedziale [0,math.pi] korzystając z równoodległych N węzłów (np. wygenerowanych używając np.linspace(0,math.pi,N)). Jak zachowuje się błąd średniokwadratowy tej interpolacji dla punktów np.linspace(0.math.pi,1000), gdy N rośnie?

2. korzystając z przykładu pokazanego na wykładzie, gdzie wartości y_i=[1,1,1,2], czy potrafisz dobrać pozycje x_i=[1,x2,x3,4] tak aby uzyskać dowolnie dużą amplitudę interpolowanego (polyfit(x_i,y_i,s=0) wielomianu na przedziale [1,4]? Czy błąd aproksymacji (polyfit, k=2) wielomianem stopnia 2 jest tak samo duży?

3. Spróbuj interpolować te same dane przy pomocy krzywych sklejanych. Czym różnią się splajny dla różnych k (k=1,2,3)?

4. Wygeneruj zaburzone obserwacje wg funkcji y=log_2(x) + scipy.stats.normal() w wielu (~1000) punktach na przedziale 1..100. Czy lepiej będzie interpolować, czy aproksymować aby znaleźć kształt funkcji?

5. Spróbuj użyć splajnów do “wygładzenia swojego rozwiązania zadania zaliczeniowego nr 1. Czy może przyda się funkcja RectBivariateSpline

ONA 5 – zagadnienie najmniejszych kwadratów

Dzisiaj na wykładzie omówiliśmy zasadniczo tematy zawarte w wykładzie 12. z metod numerycznych. Jeśli ktoś chciałby doczytać to znajdzie materiały tutaj

Zadania na lab:

0. Rozważmy trzy punkty na płaszczyźnie: (0,6), (1,0) i (2,0). Jaka prosta przechodzi najbliżej nich? ułóż układ równań linioowych, który można rozwiązać metodą najmniejszych kwadratów. Wykorzystaj funkcje scipy.linalg.qr aby zobaczyć rozkład QR macierzy A. Użyj funkcji scipy. linalg.lstsq aby znaleźć rozwiązanie. jakie jest znaczenie wartośći zwróconych przez tę funkcje?

1. Rozważmy dane:

x f(x)
0.00 4.00000000000000e+00
1.25 3.28650479686019e+00
2.50 3.08208499862390e+00
3.75 3.02351774585601e+00
5.00 3.00673794699909e+00
6.25 3.00193045413623e+00
7.50 0.00055308437015e+00
8.75 3.00015846132512e+00
10.00 3.00004539992976e+00

Jak będzie wyglądało dopasowanie met. najlepszych kwadratów funkcji f(x)=a+b*exp(-x) do tych danych?

2. Rozważmy dane o obwodzie pnia, trees-stripped ( do wczytywania przyda się funkcja scipy.loadtext). Kolejne kolumny oznaczają tu:

  • obwód pnia
  • wysokość drzewa
  • objętość pozyskanego drewna.

Spróbuj dopasować (metodą najmniejszych kwadratów objętość drzewa jako funkcję:

  • kombinację liniową obwodu pnia i wysokości drzewa
  • iloczynu wysokości przez obwód
  • kombinację liniową powyższych

Gdzie uzyskujemy najmniejszy błąd przybliżenia?

3. Przedstaw na wykresie wynik zadania 1 (wykres punktowy obserwacji i dużo gęściej próbkowany wykres liniowy znalezionej funkcji). Jaki jest problem naszego rozwiązania? czy można jakoś pomóc sobie używając ważonego problemu średnich kwadratów? jak to zrobić dla wagi jednej z obserwacji=0? a jak dla wagi 0.1?

 

 

 

ONA – zadanie zaliczeniowe 1 – krajobrazy losowe

Naszym zadaniem będzie napisanie programu generującego losowe krajobrazy fraktalne.

Losowy krajobraz jest płaszczyzną dwuwymiarową w przestrzeni trójwymiarowej. Zwkle rozważa się płąszczyzny trójkątne lub kwadratowe, my zajmiemy się tym drugim przypadkiem. Płaszczyznę będziemy reprezentować jako macierz wymiaru 2^N+1*2^N+1, gdzie M[i,j] oznacza wysokość (współrzędną Y) punktu o x=i i z=j.

Aby utworzyć losowy krajobraz tego rodzaju, będziemy postępować wg następującego rekurencyjnego schematu:

Diamond_Square.svg

Zacznijmy od zainicjowania wartości narożnych macierzy (krok 0).

Następnie wykonujemy N zejść rekurencyjnych (k przebiega od N do 1):

  • (diamond step na schemacie) ustaw wartość każdego punktu środkowego (żółty) na wartość średnią pomiędzy czterema narożnikami (czarnymi)zaburzoną o wartość 2*k*sigma* norm()
  • (square step na schemacie) Ustaw wartości czterech wartości na środkach (żółtych) brzegów kwadratu na wartości średnie pomiędzy ich czterema sąsiadami (czarnymi). zaburzoną o ((2 *k)-1)*sigma*norm()

Naszym zadaniem jest nie tylko generowanie takiego krajobrazu, ale także jego wizualizacja w postaci mapy fizycznej (imshow) i płaszczyzny 3D.

Nasz program w wersji podstawowej (10pkt) powinien dać się uruchomić jako skrypt i pozwalać na ustawienie następujących parametrów:

  • N – rząd macierzy (2^N+1). Domyślnie 10, powinno działać co najmniej do N=3014
  • sigma – stopień “górzystości” terenu
  • map_file – nazwa pliku (pdf lub png) z mapą. Jeśli podane, to mapa jest zapisywana do pliku, jeśli nie, otwieramy okienko z wykresem
  • surf_file – nazwa pliku (pdf lub png), z wykresem płaszczyzny. Jeśli podane, to powierzchnia jest zapisywana do pliku, jeśli nie, otwieramy okienko z wykresem
  • colormap – mapa kolorów z matplotlib do użycia
  • matrix_file – zapis macierzy do pliku

Wersja nieco bardziej zaawansowana (dodatkowe 2 pkt) może obsługiwać także następujące opcje:

  • wczytywanie macierzy już częśćiowo wypełnionej i losowanie tylko od kroku K
  • Okresowe warunki brzegowe, tak, aby dało się “wykafelkować” większą powierzchnię

Rozwiązania wysyłamy do 3. maja 2016 (włącznie) w postaci e-mail’a do wykładowcy z [ONA-1-2016] w tytule.

Rozwiązania po terminie otrzymują połowę punktów, na które zasługiwałyby w terminie i można je przysyłać do 6 VI 2016).

ONA 4 – układy równań liniowych w postaci macierzowej

Dziś zajmujemy się ukłądami równań liniowych w reprezentacji macierzowej.

Wykład będzie prowadzony prz tablicy, więc slajdów nie ma, ale potrzebne materiały ( i dużo więcej niż nam potrzeba) są dostępne w materiałach do wykładu z Metod numerycznych . Nas w szczególności interesują wykłądy 5 (eliminacja Gaussa) i 7 (uwarunkowanie problemu). Można też zajrzeć do wykładu 8 (macierze rzadkie), ale jest on o dla nas zdecydowanie zbyt obszerny, a jest to dla nas tylko temat niejako poboczny.

Jeśli chodzi o operacje algebry liniowej w pythonie, przydadzą nam się następujące biblioteki:

numpy.linalg – podstawowe operacje na macierzach, m.in. cond (A), det(A)

scipy.linalg – nieco więcej operacji na macierzach, m.in. hilbert(N), solve(A,B)

scipy.sparse.linalg – operacje na macierach rzadkich m.in. spsolve(A,B)

scipy.io – wczytywanie różnych ciekawych formatów – np. macierzy rzadkich w formacie matlaba loadmat(f)

Zadania na laboratorium (ew. nierozwiązane jako praca domowa):

  1. Przypomnij sobie rozwiązywanie ukł. równań poprzez rozkład LU. Zaimplementuj ręcznię tę metodę (w wersji bez wyboru el. głównego) i spróbój rozwiązać tym sposobem układ postaci A= np.array([[1e-18,1.0],[1.0,2.0]]) B=np.array([1.0,4.0]). Porównaj wynik z wynikiem scipy.linalg.solve(A,B). Dlaczego w przypadku tej macierzy wybór elementu głównego ma takie znaczenie?
  2. Korzystając z metody scipy.rand(1000,1000) stwórz macierz losową. Użyj jej do przetestowania rozkładu macierzy  LU (scipy.linalg.lu) w wersji z permutacją i bez. Rozwiąż taki układ (dla b=scipy.rand(1000)) przy pomocy scipy.linalg.solve()
  3. Załóżmy, że mamy teraz jedną macierz A (my użyjemy losowej) i bardzo wiele różnych warunków brzegowych B (my wylosujemy 1000). Spróbuj rozwiązać wszystkie układy równań powstałe z przyrównania tej samej macierzy A do wielu różnych wektorów B. Czy można zamiast korzystać wielokrotnie z funkcji solve() coś przyspieszyć? Np. korzystając z funkcji lu()?
  4. Duże macierze równań liniowych często powstają w problemach inżynierskich. Pobierz jedną z macierzy z kolekcji układów powstałych przy projektowaniu elementów samolotów Boeing (np. nr 38 ) i spróbuj wczytać ją jako macierz rzadką do ipython’a. Następnie rozwiąż ją dla losowych warunków brzegowych b.

 

http://www.cise.ufl.edu/research/sparse/matrices/Boeing/bcsstk38.html