Program labu do wykladu z RRzw 2014/15

semestr letni 2014-15

• Bieżący lab
• Projekt zaliczeniowy

Literatura (schematy)

1. Janina Jankowska, Michał Jankowski, Metody numeryczne, tom 1, Wydawnictwo Naukowo Techniczne (WNT), Warszawa, 1981.
2. Andrzej Krupowicz, Metody numeryczne zagadnień początkowych równań różniczkowych zwyczajnych., Państwowe Wydawnictwo Naukowe (PWN), Warszawa, 1986.
3. Krzysztof Moszyński, Rozwiązywanie równań różniczkowych zwyczajnych na maszynach cyfrowych., Wydawnictwo Naukowo-Techniczne (WNT), Warszawa, 1971.
4. Andrzej Palczewski, Równania różniczkowe zwyczajne. Teoria i metody numeryczne z wykorzystaniem komputerowego systemu obliczeń symbolicznych., Wydawnictwo Naukowo-Techniczne (WNT), Warszawa, 1999.

Lab

Narzędzia

Program labów

• Lab 1 (3/03/2015) - Podstawy linuxa tzn krótkie omówienie podstaw - tzn podstawowych narzędzi wywoływanych w terminalu typu: np ls - listowanie zawartości katalogu itd, prostych edytorów tekstowych, pakowaczy itp
  Applety javadfield.jar i pplane.jar - czyli interaktywne rysowanie pól wektorowych i trajektorii w 2 wymiarach.
  Uwaga W razie gdyby aplety nie zadziałały to proszę w w przeglądarce uruchomić któryś z tych: pierwszy lub drugi - niestety generują tylko obrazy pól wektorowych bez trajektorii
  Pakiety symboliczne maxima (wolne oprogramowanie) i mathematica (wolfram)
  • Utwórz katalog. W nim utwórz plik edytorem vi o nazwie test.sh z jedna linia "ls -l" - wylistuj zawartość katalogu. Spraw by skrypt był wykonywalny - wywołaj go usuń skrypt. Wylistuj zawartość katalogu. Usuń katalog.
  • Utwórz 2 pliki - połącz je w jeden o nazwie new2 - spakuj do archiwum tar (zgzipowany) -usuń pliki - rozpakuj katalog. usuń katalog - zmień nazwę pliku - przenieś do nadkatalogu. Znajdź plik o nazwie new2 z katalogu domowego.
  • Ściągnij aplety ze strony: strona z apletami javy lub bezpośrednio:
    dfield.jar i pplane.jar
  • Uruchom dfield: java -jar dfield.jar - narysuj kilka trajektorii, zmień parametry etc
    (w razie klopotów z apletami prosze spojrzeć powyżej i załadować pierwszy lub drugi - niestety generują tylko obrazy pól wektorowych bez trajektorii)
    Używając apletów (najlepiej dfield lub pplane) narysuj następujące pola wektorowe i kilka trajektorii (o ile to możliwe)
    1. y'=ay dla a=-1,1,10 etc (a<0 rozpad promieniotwórczy; a>0 np wzrost populacji
    2. y'=f(t) dla f=sin(t), t , -t, 1 etc
    3. x'=-kx, y' = kx - cy dla różnych dodatnich c,k (można to zrobić oboma apletami...)
    4. równanie wahadła y''=-ksin(x) k>0 np dla k=1 i zlinearyzowane r. wahadła y''=-ky
    5. Inne przykłady z ćwiczeń.
  • Uruchomić pakiet symboliczny maxima i z niego wyjść za pomocą polecenia exit();
  • Przeczytać tutorial do maximy (a przynajmniej przejrzeć) w szczególności rozdział dotyczący RRzw w maxima'ie.
  • zdefiniować równanie różniczkowe y'=y i je symbolicznie rozwiązać w maximie tzn. znaleźć rozwiązanie ogólne i spełniające y(0)=-10
  • sprawdzić czy maxima rozwiąże równanie y''=-y oraz y''=-sin(y)
  • sprawdzić czy maxima rozwiąże równanie y'=(sin(y)+y)*(exp(-t*t)) czy ogólnie postaci y'=f(y)g(t) dla f,g funkcji skomplikowanych (np których nie potrafimy scałkować)
• Lab 2 i 3 (17 i 31 III 2015) - Wstępne zapoznanie się z octavem - octave jako kalkulator naukowy, operator dwukropek : tworzenie macierzy wprost i funkcjami, wektorów, zapisywanie/czytanie do/z plików w obu formatach - tekstowym i binarnym, tworzenie macierzy z podmacierzy, wycinanie podmacierzy itp, podstawowe operacje na macierzach - mnożenie, dodawanie, transponowanie, funkcje od macierzy, normy wektorów/macierzy. Wykresy funkcji (1 wymiarowe)
  Podstawy octave'a cd. Skrypty i funkcje (m-pliki) w octavie. Instrukcje warunkowe (if else endif; switch case endswitch), pętle (while( ) do .. endwhile; do .. until( );for .. endfor). Wskaźniki do funkcji (handle - operator @). Rozwiązywanie równań nieliniowych, numeryczne odwracanie funkcji - uzywając fsolve() i zeros().
  • Użyj octave jako kalkulatora.
  • Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x8 C której pierwsze 3 kolumny to A a kolejne to B. Teraz z tej macierzy 'wytnij' podmacierz D składającą się z 1 głównego minora tzn 3x3 od C(1,1) do C(3,3). Zamień kolejność kolumn D. Wstaw D z powrotem do C jako główny minor. Policz sin od D. Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a jako DD. Policz normy (różne) macierzy DD.
  • Powtórz powyższe dla macierzy losowej. Czy zapisywanie/odczytywanie tekstowe zmniejsza dokładność?
  • narysuj wykres sin(x) na [-1,4] i na [-1000,1000].
  • narysuj wykres 2 funkcji rożnymi kolorami na 1 wykresie np sin na czerwono i cos na niebiesko na [-1,3]. Podpisz wykresy. Zaznacz punkty przecięcia (znamy je...dla x=pi/4,) jako +.
  • narysuj wykres funkcji odwrotnej do (x+sin(x))^p na [0,1] dla p=2,3,4,..,10 BEZ użycia solvera nieliniowego tzn BEZ rozwiązywania równan
  • Skrypty - czyli pliki tekstowe z komendami octave'a - wywołuje się linii komend octave'a: np rrzbasic wywola skrypt rrzbasic.m .
  • Zapoznaj się z kodem w rrzbasic.m - prostym skryptem z przykładami.
  • Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo - czyli bez użycia pętli).
  • Zsumuj liczby od 1 do n z pętlą while i for. Zrób to samo be z pętli.
  • Utwórz skrypt w octave'ie rozwiązujący LZNK dla danych (x_k,y_k) znajdujący a,b najlepiej je przybliżające w sensie arg min \sum_k |ax_k+b-y_k|^2 z użyciem operatora backslasha.
  • napisz swój skrypt (jakikolwiek np. dodający 2+2). Wywołaj go z linii komend.
  • korzystając z funkcji vander() (tworzy macierz Vandermonde'a) oraz polyval() - oblicza wartości wielomianu - rozwiąż zadanie regresji liniowej tzn dla danych x_k,y_k znaleźć a,b takie \sum_k |ax_k+b -y_k|^2 osiąga minimum. Napisać swój skrypt z rozwiązaniem tego zadania i wywołaj go z linii komend octave'a. Zbadaj uwarunkowanie macierzy Vandermonde'a. Porównaj z wynikami funkcji polyfit().
  • napisz funkcje w m-pliku która dla 2 argumentów - wektorów tego samego wymiaru zwróci wektor współczynników. rozwiazujacy zadanie liniowej regresji (modyfikując skrypt z porpzedniego zadania..)
  • przeczytaj help do fsolve() i fzero() - rozwiąż równanie x^2-2=0 z ich pomocą.....
  • rozwiąż równanie cos(x)=0 na [0,2] z użyciem solverów octave'a.
  • Napisz skrypt obliczający funkcje odwrotna do cos na [-1,1] rozwiązując równanie g(y)=y-cos(x)=0 ze względu na y na siatce równomiernej na tym odcinku. narysuj wykres. porównaj wynik z funkcja octave'a arccos(x)
  • narysuj wykres funkcji odwrotnej do (x+sin(x))^p na [0,1] dla p=2,3,4,..,10 na siatce równomiernej rozwiązując odpowiednie równanie. Porównaj z wynikiem z zadania labu.
  • Dla równania y'=(1+x)/(1-0.5sin(y))z y(0)=s narysuj kilka rozwiązania dla s=-1,0,1 na [0,1] rozwikłując odpowiednia funkcje numerycznie.
    
    nonlin.m - kilka rozwiazan zadan z uzyciem nieliniowych solverow - w tym wykres funkcji odwrotnej na siatce rownomiernej.
• Lab 4 (13 IV 2015) - Rozwiązywanie RRzw w octave'ie czyli funkcja lsode. - przypadki skalarne i wielowymiarowe; rysowanie pola wektorowego quiver() i krzywej w 3 wymiarach plot3()
  • przeczytaj help dla lsode()
  • Rozwiązać kilka prostych równań różniczkowych używając lsode() np.
    1. y'=-a*y;y(0)=-y na [0,1] albo [-10,0].
    2. Równanie eksplodujące: y'=y*y; y(0)=1 na odcinku [0,1/2] czy [0,3].
    3. y'=(cos(x))^2 z y(0)=1
    4. y'=1+y*y z y(0)=0(jakie są rozwiązania?) na [0,T] dla T=1,10,100.
    5. y'=1+x/(1-0.5sin(y))z y(0)=s narysuj kilka rozwiązania dla s=-1,0,1 na [0,1]
    6. Równanie 2-wymiarowe y''=-y z y(0)=1;y'(0)=0 na [0,2*pi](rozwiązanie cos(x)). narysuj wykresy (t,y(t) i (t,y'(t) oraz trajektorie (y,y').
    7. równanie wahadła y''=-sin(y) z y(0)=1;y'(0)=0 - narysuj wykresy (t,y(t) i (t,y'(t) oraz trajektorie (y,y')
    8. Równanie wahadła z oporem powietrza y''=-a*sin(y)-b*y' z y(0)=1;y'(0)=0 - a,b>0 - narysuj wykresy y, y' i trajektorii y vs y' - narysuj wykres a->y(1,a) i y->y(1,b) czyli wykres zalenoznosci od konkretnych parametrów Powtórz dla równania zlinearyzowanego (y''=-ay'-by')
    9. Dodaj do równania wahadła sile zewnetrzna wiatru np f(t)=d/(1+c(t-1)^2) d, c stale
  • Narysować wykres funkcji F(s)=y(1,s) zdefiniowanej jako F(s)=y(1,s) gdzie y(t,s) rozwiązanie równania y'=f(y,t) z y(0)=s dla s w [0,1], [-1,1] dla f=y,10y,-10y 0.1*y*y itp
  • Przy pomocy lsode i fsolve znaleźć s0 takie że rozwiązanie równania y'=y*(1-y)*sin(y*y+x) z war pocz. y(0)=s0 spełnia y(1,s0)=0.5. Narysować wykres F(s)=y(1,s) - potok fazowy dla t=1 i y(t,s0) dla t w [0,1]. Skąd wiemy że istnieje rozwiązanie?
  • Narysuj pole wektorowe wahadła i zlinearyzowanego wahadła tzn y''=-ay i y'=-asin(y) dla a=0.1,1,10 na [-1.5,1.5]
  • Narysuj pionowa spirale tzn krzywa: (sin(a*t),cos(a*t),t) dla t w [0,20] dla a=0.1,1,10
  
  lablsode.m - kilka rozwiazan zadan z uzyciem lsode() - wykresy zaleznosci od warunku poczatkowego s dla t=1 i a=0,1dla y'=y(1-ay) y(0)=s
• Lab 5 - schematy Eulera, Taylora i midpoint - najprostsze schemat rzędu 2 W schematach wielokrokowych za x_1,x_2 etc weź dokładne rozwiązanie jak je znamy jak nie weź za odp. pierwsze x_k k=1,..,p- p-1 kroków metody otwartej tego samego rzędu np rzędu 2 np Taylora dla metody midpoint. Eksperymentalne badanie rzędu.
  • Zaimplementuj metody Eulera w skrypcie - dla rozwiązania równania y'=ay y(0)=1 dla a=1,10,100,-1,-10,-100 na [0,T] (rozwiązania znamy...)- sprawdź jak działają tzn narysuj wykresy dla T=1,10,100, policz błąd absolutny dla t=T ae(T;h)=|y(T)-y_h(T)| i względny re(T;h)=e(T;h)/y(T) oraz dla t=h tj: ae(h;h)=|y(h)-y_h(h)| i re(h,h)=e(h;h)/y(h). Narysuj wykresy błędów ae(t;h) i re(t;h) na [0,T].
  • Powtórz zadania poprzednie dla równania drugiego rzędu y''=-y z y(0)=0;y'(0)=1 - narysuj dodatkowo wykresy trajektorii (y,y') czy obliczone trajektorie są zawarte w okręgach? (W schemacie zamkniętym Eulera musimy odwracać macierz 2x2 ze stałymi współczynnikami. Proszę nie liczyć macierzy odwrotnej - choć to w tym przypadku prostym to ma akurat sens ale w ogólnym nie)
  • Zastosuj schematy Eulera dla równania y'=1+y' z y(0)=0 (rozwiązanie znane) na odcinku [0,T) dla T=0.5,1,1.5. Policz błąd dla t=0.5,1.5 dla różnych h=0.5,1/4,1/8 etc.
  • Zaimplementuj metodę midpoint i schemat Taylora rzędu 2 w skryptach - dla rozwiązania równania y'=ay y(0)=1 dla a=1,10,100,-1,-10,-100 na [0,T] (rozwiązania znamy...)- sprawdź jak działają tzn narysuj wykresy dla T=1,10,100, policz błąd e(T;h)=|y(T)-y_h(T)| i ew(T;h)=e(T;h)/y(T) oraz e(h;h)=|y(h)-y_h(h)| i ew(h,h)=e(h;h)/y(h). Narysuj wykresy błędów ae(t;h) i cew(t;h) na [0,T]. Porównaj wykresy i błędy z wynikami dla schematów Eulerem otwartym/zamkniętym.
  • Powtórz zadanie poprzednie dla równania drugiego rzędu y''=-y z y(0)=0;y'(0)=1 - narysuj dodatkowo wykresy trajektorii (y,y') czy obliczone trajektorie są zawarte w okręgach? Powtórz zadanie dla równania wahadła y''=-sin(y) czy obliczone rozwiązania są okresowe? porównaj z rozwiązaniami z lsode().
  • Przetestuj schemat midpoint dla y'=-y z y(0)=1 na [0,T] dla coraz większych T dla h=0.1,1e-2,1e-3 etc narysuj wykresy.
  • Przetestuj rząd błedu lokalnego schematu Eulera, Taylora i midpoint dla równania y'=ay y(0)=1 (dla różnych a=0.1,1,10,-0.1,-1,-10 etc) dla t=1,10,100 tzn. Dokładniej wstaw wartości rozwiązania x(t+kx) w miejsce x_{n+k} itd w schemat x_{n+k}-\Phi(t+nh,h,x_n,...,x_{n+k-1}=0) tzn policz e(t;h)=|x(t+kh)-\Phi(t,h,x(t),...,x(t+(k-1)h))| (np dla otwartego Eulera policz e(t;h)=|x(t+h)-x(t)-hf(t,x(t))|) dla ustalonych t, powtórz to samo dla h/2 tj. e(t,h/2) i policz iloraz tzn. e(t;h)/e(t;h/2). W przypadku schematu rzędu p powinniśmy dostać ok. 2^{p+1}
  • Przetestuj eksperymentalnie rzędy błedu schematu Taylora lub midpoint dla różnych równań ze znanymi rozwiązaniami : y'=ay,y(0)=1; a=-100,-10,-1,1,10,100 dla t=1,10,100 (o ile zakres arytmetyki pozwoli) to samo dla y''=y y(0)=0 y'(0)=1; y'=cos^2(y) y(0)=0 dla t=1 i 100 itd Można dla równania wahadła biorąc rozwiązania lsode() za dokładne.
  
  Tu zamieszczam m-plik ze schematami otwartym Eulera, midpoint, taylor (rz. 2) i ich testami - rzad lokalengo bledu etc:
  exEuler.m - schemat otwarty Eulera (działa w wielu wymiarach)
  taylor.m - schemat Taylor
  midpoint.m - schemat midpoint
  testexEul.m - proste testy schematu Eulera (skalarne i 2-wymiarowe)
  OrderExEuler.m - testy rzędu (lokalnego błędu) schematu otwartego Eulera
  OrderExEuler2d.m - testy rzędu (lokalnego błędu) schematu otwartego Eulera dla zadania 2 rzedu
  testmidpoint.m - testy schematu midpoint rząd zbieżności i niestabilność dla dużych T (dx/dt=-x z x(0)=1)
  Ordermidpoint.m - testy rzędu (lokalnego błędu) schematu midpoint
  OrderTaylor.m - testy rzędu (lokalnego błędu) schematu Taylora
  testtaylor.m - testy schematu Taylora (w 2 wymiarach)
• Lab 6 Modelowanie - budynek w czasie trzesienia ziemi. Zalozmy ze sztywne pietra w budynku sa polaczone elastycznie - czyli jesli pietro i+1 sze o masie m_{i+1} wysunie sie wzgledem pietra l o x_{i+1}-x_{i} to na oba dziala sila sciagajaca o wartosci proporcjonalnej do przesuniecia tzn k_i(x_{i+1}-x_i) tzn na i-te dzialaja 2 sily k_i(x_{i+1}-x_i) i -k_{i-1}(x_{i}-x_{i-1}), na pierwsze (parter) dziala analogiczna sila sciagajaca miedzy ziemia a tym pietre, tzn -k_0x_1 i oczywiscie sila sciagajaca miedzy parterem a kolejna kondygnacja. W czasie trzesienia ziemi dodatkowo dziala sila f(t)=Gsin(gamma t) na parter dla gamma czestotliwosci trzesienia np. rownej trzy.
  • Napisz rownanie modelujace przesuniecia pieter w postaci MX''=SX+ CX' + F(t). (F =(f(t),0..,0)^T, C - macierz oporu np powietrza przyjmijmy ze diagonalna. Napisz procedury liczace macierze A=M^{-1}S, M^{-1}C i funkcje M^{-1}F
  • Dla budynku 2 pietrowego o m_k=m=10000kg z k_l=k=5000kg/s^2 policz wartosci wlasne macierzy M^{-1}S : \lambda_l - a nastepnie czestotliwosci tego budynku zdefiniowane jako \omega_l=\sqrt{-\lambda_l} i okresy budynku tzn T_l=2\pi/\omega_l. Czy okresy sa w zaresie [2,3] (typowy zakres okresów trzesienai ziemi)
  • zasymuluj budynek 2 pietrowy przyjmujac losowe warunki poczatkowe (male) bez sily zewnetrzenej. Narysuj wykresy przesuniec pieter.
  • Powtorz zadanie dla budynku 10 pietrowego (o pietrach o masie 10000kg) i k_l=k=5000kg. Okresl czy okresy budynku sa poza [2,3]. Zasymuluj rozwiazania na odcinku [0,200]. Powtorz dodajac tlumienie tzn C=alfa*Id dla alfa=0.1,1,10
  • Policz czestotliwosci budynku w ksztalcie piramidy tzn takiego ze pierwsze pietro wazy 10000kg a kazde kolejne o 20% mniej (k stale = 5000kg/s^2) Porownaj ze zwyklym budynkiem 10 pietrowym o m_i=m=10000kg i k_i=k=5000gh/s^2 - policz czestotliwosci - Zasymuluj trzesienie ziemi (war poczatkowe zerowe) gamma=3 a G=10,20,40,80cm. narysuj wykresy ruchu najwyzszego pietra w obu przypadkach - policz max przesuniecia Dodaj tlumienie alfa=0.1,1,10. Zmniejsz sile zewnetrzna tzn niech G = G(t) np G_{max}exp(-t^2) czy funkcja stala przez 10 sekund a potem gwaltownie schodzaca do zera (sami dobierzcie Panstwo odpowiednia)
  
  earthquake.m - skrypt z rozwiazaniami modelu trzesienia ziemi >
• Lab 7 testowanie rzędu schematów i rzędu zbieżności. Start dla schematów wielokrokowych i ocena projektu
  • Napisz m-plik ze schematem Adamsa-Bashforda drugiego rzędu - rozwiązujący dane zagadnienie początkowe na odcinku [t_0,T] na podziale na N+1 punktów tzn wyznaczający wektor x=(x_k)_{k=0}^n z x_k \approx x(t_k) dla t_k=t_0+k*h z h=(T-t_0)/n. Na wejściu wskaźnik do pola wektorowego, t_0,x_0,T,n, i x_1. Na wyjściu wektor x. W praktyce x_1 musi być obliczony jakimś schematem 1-krokowym odpowiedniego rzędu.
  • Przetestuj rząd lokalnego błędu (rząd schematu) otwartego schematu Adamsa-Bashforda rzędu dwa dla y'=ay z y(0)=1 z x_1=exp(ah) dla a=1,-1.
  • Przetestuj rząd zbieżności metodą połowienia kroków otwartego schematu Adamsa-Bashforda rzędu dwa dla y'=ay z y(0)=1 dla a=-1,1 dla t=0.1,1,10,100. Biorąc różne x_1: x_1=exp(ah), x_1 obliczone schematem 1-krokowym rzędu dwa np Heuna i otwartym schematem Eulera x_1=x_0+h*f(t_0,x_0). Czy widać różnicę w wynikach? (Liczymy błąd e(t,h)=|x_n-sol(T)| dla h=(T-t_0)/n dla połowionych h; następnie liczymy stosunek e(T;h)/e(T;h/2)).
  • Przetestuj rząd schematu i rząd zbieżności schematów Heuna i zmodyfikowanego Eulera dla y'=ay z y(0)=1 z x_1=exp(ah) dla a=1,-1 dla t=0.1,1,10,100.
  • Przetestuj rząd schematu i rząd zbieżności schematów Heuna lub zmodyfikowanego Eulera i Adamsa otwartego rzędu dwa dla y'=1+y^2 z y(0)=0 (dla Adamsa x_1 dostateczne dokładne) dla t=0.1,1,1.5, 1.54. Czy widać eksplozję wyniku blisko \pi/2?
  • Przetestuj rząd schematu i rząd zbieżności schematów Heuna lub zmodyfikowanego Eulera i Adamsa otwartego rzędu dwa dla y''=-y'' z y(0)=0 y'(0)=1 (dla Adamsa x_1 dostateczne dokładne) dla t=0.1,1,10,100,1000. Czy widać okresowość trajektorii? Powtórz zadanie dla y'=-sin(y) z y(0)=0 y'(0)=1.
  • Policz dla y''=y*y +a y(0)=1 pochodna dy/da(x) tzn pochodna względem parametru a dla a=0 i x=0.1 raz używając lsode i ilorazu różnicowego a raz używając równań w wariacjach (i tez lsode())
  
  adams2s.m - schemat jawny Adamsa-Bashfortha 2-krokowy
  Orderadams2s.m - testy rzedy schematu A_B 2krokowego
  testadams2s.m - testy zbieznosci schematu Adamsa-B 2krokowego (w 2 wymiarach)
• Lab 8 Portrety fazowe w 2D i 3D i metoda strzałów.
  • Dla najprostszej wersji równania ofiara-drapieżnik : x'=(a-alfa y)x, y'=(-c+gamma x) y np. \quad a=4,alfa=2,c=gamma=3 narysuj portret fazowy dla (x,y) w kwadracie [0,c/gamma] X [0,a/alfa] tzn. narysuj kilka trajektorii np. w aplecie pplane (w octave'ie: lsode() oraz dla pola wektorowego - quiver(), funkcja maximy plotdf(), istnieją inne narzędzia on-line). Następnie narysuj porterety fazowe zlinearyzowanych równań w pkt stacjonarnych tj. (0,0) i (c/gamma, a/alfa) . Policz wartości własne macierzy pochodnych w pkt stacjonarnych. Określ ich stabilność na bazie portretu fazowego ale również wartości własnych.
  • narysuj pola wektorowe y'=Ay w 2 wymiarach dla rożnych macierzy - kanoniczne tzn dla diagonalnych z diagonala [a1,a2], [a1,-a2], a1,a2, dodatnie np a1=1 a2=0.1,1,10 z klatka Jordana [a,1;0,a] dla a=0.1,1,10,100, z para zespolonych w.własnych [a,b;-b,a] dla a,b nieujemnych np weź a=0 i b=1 lub a=1 i b=0.1,1,10,100. W każdym przypadku znajdź kilka rozwiązań i naszkicuj trajektorie. (Mozna uzyc roznych narzedzi....
  • naszkicuj kilka trajektorii dla y'=Ay z A=CJC^{-1} dla danej macierzy C nieosobliwej np C=[1,1;1,2] i J w postaci kanonicznej - wystartuj z y0 na prostych o kierunkach kolumn C tz np y0=[1,;1] lub [-2,;4]
  • Naryuj kilka trajektorii w 3D dla rownania y'=Ay z y(t0) losowym w postaci kanonicznej tzn albo diagonalnej albo z 1 klatka Jordana wymiaru 2 lub 3 - za warunki startowe wezmy losowe wektory z [-1,1]^3 np
    1. A=[-c,0,0;0,-a,-b;0,b,-a] dla a,b,c dodatnich np a=1,b=10,c=1 -zobacz co sie stanie jak zwiekszymy b i zmniejszymy c
    2. A=[-c,1,0; 0-c,1;0,0,-c] klatka Jordana 3x3 dla c=0.1,1,10
    3. A=[-c,0,0;0,-a,0;0,0,b]- 2 ujemne wartosci wlasne 1 dodatnia
    4. A=[-c,0,0;0,-a,eps;0,0,-a] - 1 klatka Jordana 2x2 za eps=0.1,1,10
  • narysuj kilka trajektorii dla y'=Ax -b dla A ujemnie okreslonej losowej i danego wektora stalego b - biorac za b=A[1;1]. czy z portretu fazowego jestesmy w stanie okreslic pkt stacjonarny ([1;1])?
  • Znajdź rozwiązanie równania y'=Ay y(0)=y0 używając expm() i lsode() porównaj wyniki z dokładnymi. Za A weź macierze [1,1;1,-1], [0,1;-1,1], [1,3,3,1] itp [-1,10;1,-1]. naszkicuj pola wektorowe. narysuj trajektorie czy wszystko się zgadza z teoria?
  • zaimplementuj metodę strzałów dla liniowego zadania y''-d(x)y'- c(x)y=f(x) z y(a)=ya; y(b)=yb, a, b, ya,yb dane liczby, c,f dane funkcje na [a,b] (2 strzały z pochodną s_1=0 i s_2=1 tzn rozwiązujemy równanie z warunkami początkowymi y(a)=ya; y'(a)=s_k otrzymujemy rozwiązanie dla t=b tzn y(b;s_k) k=1,2 i stąd wyliczamy s takie że rozwiązanie z y(a)=ya; y'(a)=s_k spełnia y(b)=tb o ile ono istnieje - zazwyczaj tak).
  • rozwiąż metodą strzałów równanie -y''+y=0 z y(0)=1; y(b)=1 dla b=1,4,10,15,20. narysuj wykresy.Policz błąd |y(b;s)-yb|. Raz użyj funkcji z poprzedniego zadania a raz wylicz na papierze s (daje się!) i policz lsode() rozwiązanie dla tego s. Następnie sprawdź czy to oba policzone rozwiązania spełniają y(b)=1.
  • rozwiąż metodą strzałów zagadnienia brzegowe:
    1. -y''+y=f z y(0)=0; y(b)=0 dla f=2*sin(x) i b=\pi (znamy rozwiązanie - zadanie - można też zgadnąć).
    2. -y''+y=-2+(1-x*x) z y(-1)=0; y(1)=0. (znamy rozwiązanie - zadanie - można też zgadnąć)
    3. -y''+exp(-x) y =cos(x) z y(-1)=0; y(1)=0.
    4. -y''+ay'+ y =cos(x) z y(-1)=0; y(1)=0 dla a=0,1,10,100,-1,-10,-100.
    Narysuj wykresy. Policz błąd |y(b;s)-yb|.
  • Zaimplementuj metodę strzałów dla równania y''=F(t,y,y') z z y(a)=ya; y(b)=yb. Korzystając z funkcji lsode() i fsolve() (lub swojego solvera dla równań nieliniowych np. metody Newtona - jak liczyć pochodną F(s)=y(b;s) po s?). Przetestuj dla F(t,y,y')=sin(y) z y(0)=1 y(1)=2 rysując rozwiązania dla kolejnych iteracji s_k (to jakoś wymaga uzyskania s_k kolejnych iteracji fsolve().).

  ---

  Skrypty octave'a

• nrrbasic.m -prosty skrypt octave'a z kilkoma podstawowymi operacjami
• nonlin.m - kilka rozwiazan zadan z uzyciem nieliniowych solverow - w tym wykres funkcji odwrotnej na siatce rownomiernej.
• lablsode.m - kilka rozwiazan zadan z uzyciem lsode() - wykresy zaleznosci od warunku poczatkowego s dla t=1 i a=0,1dla y'=y(1-ay) y(0)=s
• exEuler.m - schemat otwarty Eulera (działa w wielu wymiarach)
• taylor.m - schemat Taylor
• midpoint.m - schemat midpoint
• testexEul.m - proste testy schematu Eulera (skalarne i 2-wymiarowe)
• OrderExEuler.m - testy rzędu (lokalnego błędu) schematu otwartego Eulera
• OrderExEuler2d.m - testy rzędu (lokalnego błędu) schematu otwartego Eulera dla zadania 2 rzedu
• testmidpoint.m - testy schematu midpoint rząd zbieżności i niestabilność dla dużych T (dx/dt=-x z x(0)=1)
• Ordermidpoint.m - testy rzędu (lokalnego błędu) schematu midpoint
• OrderTaylor.m - testy rzędu (lokalnego błędu) schematu Taylora
• testtaylor.m - testy schematu Taylora (w 2 wymiarach)
• earthquake.m - modelowanie trzesienia ziemi
• adams2s.m - schemat jawny Adamsa-Bashfortha 2-krokowy
• Orderadams2s.m - testy rzedy schematu A_B 2krokowego
• testadams2s.m - testy zbieznosci schematu Adamsa-B 2krokowego (w 2 wymiarach)
• testAB2start.m - testy wartości startowej (x^h_1) dla explicit Adams-Bradford scheme (of order 2)

Projekt zaliczeniowy

-wykresy zazleznosci od k wartosci x dla t=200 i analogiczny wykres pochodnej wzgledem k