Laboratorium z Równań Różniczkowych zwyczajnych labem - do wykładu prof. Maksymiliana Dryji semestr letni 2005/06

wiosna 2006, czwartki godz. 16:15 (co drugi tydzień, pierwsze zajęcia: - termin mogą się zmienić w trakcie semestru np z powodu świąt itp).

Będziemy posługiwali się przede wszystkim octave'm , maplem (niestety starą wersją) i maximą. Maxima to darmowy klon (przynajmniej do użytku edukacyjnego) maple'a czyli pakiet obliczeń symbolicznych, octave to darmowy klon matlaba - pakietu do obliczeń numerycznych.

Wstępny plan labu (może się zmienić i to znacząco)

Lab 1 - sprawdzenie kont na tmoko niezbędnych do uruchomienai maple'a; rysowanie pól wektorowych i przybliżonych rozwiązań przy pomocy apletu java DFIELD na stronie www: http://math.rice.edu/~dfield/ , w szczeg. proszę o Zad1: przy pomocy apletu narysować pole wektorowe dla równania y'=|y|^{1/2} i narysować trajektorie przechodzące przez (0,0) i (0,1). Czy mamy jednoznaczność? Zad 2: dla równania y'=x + x/y narysować pole i przy pomocy różnych solverów narysować kilka trajektorii startując z punków z prawej górnej ćwiartki ekranu. Czy wszystkie działają poprawnie? To samo dla równania y'=-x/y. Zad 3: narysować pola wektorowe dla dotychczasowych przykładów z ćwiczeń i wykładu. 2 marca 2006.
Lab 2 - maxima i maple - podstawy - w szczególności jak rozwiązać proste równanie różniczkowe symbolicznie; tu plik tekstow z kilkoma przykladami jeśli chodzi o mnie nic więcej z maximy umieć nie trzeba: maxima-desolve.txt ; przykładowa sesja maple'a: sesja1.ms ( wersja tekstowa: sesja1.txt) - funkcja dsolve - rozwiązująca symbolicznie równania różniczkowe - znół wystarczy umieć się posługiwać tą funkcją, ewent. narysować wykres funkcji w maple'u., czy sprawdzic czy jakaś funckja jest rozwiązaniem. 16 marca 2006
Lab 3 - ode skalarne w maximie i maplu- rozwiązywanie problemów z ćwiczeń symbolicznie w maximie i maplu kontynuacja: rozwiązać równania z ćwiczeń; sprawdzenie czy dana funkcja jest rzeczywiście rozwiązaniem; czy spełnia warunek początkowy itp; sprawdzanie czy równanie jet różniczką zupełna, znajdowanie funkcji pierwotnej - sprawdzenie czy dana funkcja jest f. pierwotną dla danej różniczki! spraw.ms i spraw.txt - skrypty maple'a sprawdzające czy rozwiązanie jest OK i rozup.ms i rozup.txt - różniczka zupełna w maple'u 30 marca 2006
Lab 4- testowanie eksperymentalne rzędu oraz błędu np dla równania y'=ay y(0)=1 a=-10,-1,1,10., dla kilku prostych schematów (Eulera otwarty/zamknięty), midpoint x(n+2)=x(n)+2hf(n+1); x(0)=x0, x(1) - dobre przybliżenie rozwiązania w t0+h , Heuna: x(n+1)=x(n) +0.5h*(f(n)+f(t(n+1),K)) gdzie K= x(n)+2hf(n+1), otwarty Adamsa-Bashfortha rzędu 2: x(n+2)=x(n) +0.5h(3*f(n+1)-f(n)) itp programowanie schematy predyktor - korektor w Pascalu/C na bazie pary Euler otwarty zamknięty w 1 wymiarze (+np m. bisekcji czy m. iteracji Banachowskich). A jak się uda to może i 2 wymiarach. 20 kwietnia 2006
Lab 5- równania różniczkowe w w octavie, porównanie schematów otwartych i zamkniętych. sesja5.m - całkowanie równań różniczkowych zwyczajnych w octavie (2 proste przykłady). 4 maja 2006
Lab 6 - równania różniczkowe w w octavie, porównanie schematów otwartych i zamkniętych - tzn scałkować układ x'=998*x+1998*y;y'=-999*x-1999*y z war. pocz. x(0)=y(0)=1 an odcinku [0,10], [0,100] itp otwartymi i zamkniętymi schematami z lsode - tzn ustawiamy komendą lsode_option() rodzaj schematu np lsode_options("integration method","adams") ustawia otwarte schematy Adamsa, i wywołujemy tic;lsode;toc (tic ustawia timer - toc podaje czas który upłynąl od ostatniego tic). Implementacja zamkniętego schematu Eulera z wykorzystaniem funkcji fsolve w Octave'ie - w 1 i 2 wymiarach.; tzn zaimplementować schemat i przetestować na układzie r. sztywnych; porównać z otwartym Eulerem. skrypt octave'a (napisany w pospiechu) rozwiazujacy to ostanie zadania EZ.m 18 maja 2006
Lab 7 oddanie projektu zaliczeniowego. Portrety fazowe dla równań liniowych 2-wymiarowych w PPLANE

Projekt zaliczeniowy - termin ostani lab później obniżona ilośc punktów (2 propozycje do wyboru - moze będzie więcej) - istnieje możliwość zaliczenia innego projektu - możecie Pańśtwo przestawić propozycje do mojej ewent. akceptacji.

Projekt 1 Metoda strzałów. Mamy równanie y'=f(x,y) y(t0)=g(s) - znaleźć s i w szczeg. y(t,s) dla t w [t0,T] takie aby y -rozwiązanie spełniało G(y(T,s))=0 np y(T,s)=0. Rozwiązać numerycznie dla

Zadanie liniowe 2-go rzędu y''=4y+f(x), y(0)=0, y'(0)=s, z warunkiem y(1)=0 dla różnych f=0, 1+x*x, exp(3x), x/(1+x*x). Jako wynik podać s, błąd |y(1,s)| i narysować wykres y(t,s). (f. plot w octave'ie czy gnuplotem - albo jeszcze inaczej). Sprawdzić w maple'u czy maxima'ie (czy dfieldzie lub octavie'ie o ile rozwiązanie nie jest wyrażone przez funkcje elementarne) że rozwiązanie z tymi warunkami początkowymi rzeczywiście spełnia y(1,s)=0 (no y(1,s) bliskie 0).
Zadanie nieliniowe skalarne y'=(1-y*y)sin(y*y+x), y(0,s)=s z warunkiem y(1,s)=0. Jako wynik podać s, błąd |y(1,s)|, max i inf (dyskretne dla tn=h*n w [0,1]) rozwiązania y(t,s) i narysować wykres y(t,s).

Jako solver mozna zastosować samodzielnie zaimplementowany schemat rzędu co najmniej 4 (np Rungego) albo funkcję lsode z Octave'a (będzie na 5 labie - polecam) czy jakąś biblioteke np cvode.
W (1) mamy zadanie liniowe tzn y(1,s) jest funkcją liniową tzn y(1,s)=a+bs - wystarczy wyznaczyć a i b (dla 2 różnych s1 i s2 znaleźć y(1,s1)=a+bs1 i y(1,s2)=a+bs2 i rozwiązać układ równań liniowych) i potem wyliczyć s. W (2) funkcja y(1,s) jest nieliniowa - i trzeba albo zastosować jakąs metodę rozwiązywania r. nieliniowych (ale skalarnych) np bisekcji (tu warunki met. bisekcji są spełnione bo y(1,-1)=-1 i y(1,1)=1.... albo co oszczędzi sporo pracy funkcję octave'a fsolve rozwiązywania równań nieliniowych (ale trzeba ją poznać -> Lab 5 - o ile starczy czasu) czy znów jakąś bibliotekę (kinsol, minpack etc - ale to dla osób wprawnych).

Należałoby też uzasadnić teoretycznei dlaczego powyższe problemy mają rozwiązanie i czy ono jest wyznaczone jednoznacznie (w przypadku (2) już Państwo są w stanie to zrobić a w (1) po teorii r. liniowych.

Wykresy albo funkcja plot w octavie'ie albo przy użyciu gnuplota (polecam) albo samodzielnie (odradzam).

Można zauważyć że z wykorzystaniem octave'a ten projekt nie wymaga wielkiego wysiłku.

Projekt 2 Zaprogramować (C, C++, Pascal) schemat o zmiennym kroku całkowania (bazujący na schemacie Rungego rzędu 4) opisany w ksiażce A. Palczewskiego Równania Różniczkowe zwyczajne, WNT ,Warszawa 1999 rozdział 5.2, str. 151-155. W bibliotece jest też skrypt na którym oparta jest książka. Testować dla różnych równań porównując wyniki ze schematem Rungego 4 i z lsode (z Octave'a) tzn. podać o ile wyniki się różnią w końcach przedziału całkowania (biorąc dla poziomu dokłądności E - h=(E)^{1/3} w sch. Rungego) i podać ilość obliczeń pola wektorowego-prawej strony równania - dla obu schematów.
Przetestować schemat (dla różnych E - poziomów dokładności np 1e-3, 1e-4 itd ) na:
I. Rówaniu liniowym y''=-y z y(0)=0, y'(0)= -1 na odcinku [0,10] (to tak czy w ogóle działa)
II Równaniu van der Pola :
y''-a*(1-y*y)*y'+y=0 0<a - parametr.

Sprawdzić przypadki (dla różnych E)

a=1, y(0)=2, y'(0)=0 na odcinku [0,20]
a=1000, y(0)=2, y'(0)=0 na odcinku [0,3000]

Porównać ze schematem Rungego rzędu 4 i lsode z octave tzn porównać wyniki w końcach odcinków czyli w (I) dla t=10, (II.1) t=20 a w (II.2) t=3000 oraz ilość wywołań pola wektorowego (dla obu schematów - choć jak ktoś potrafi to można i dla lsode).

Narzędzia podstawowe

maxima - pakiet obliczeń symbolicznych dostępny za darmo dla wszytkich, manual maximy -- wersja on-line; również jako plik pdf

maple - pakiet do obliczeń symbolicznych (komercyjny mamy starą wersj e na wydziale zainstalowana na komputerze tomoko)

octave - pakiet obliczeń numerycznych (klon matlaba pakiety najczęściej używanego w praktycznych obliczeniach)

dokumentacja w htmlu

DFIELD - narzędzie w javie umożliwijące rysowanie pól wektorowych skalarnuch (nieautonomicznych) w przeglądarce
PPLANE- narzędzie w javie umożliwijące rysowanie pól wektorowych 2-wymiarowych autonomicznych w przeglądarce

Narzędzia dodatkowe

gnuplot - narzędzie Unixa pozwalające na wykonywanie prostych rysunków np wykresów funkcji

matlab - pakiet do obliczeń numerycznych

mupad - pakiet obliczeń symbolicznych; dostępny niegdyś w starszej wersji za darmo dla wszytkich, a do celów edukacyjnych bez ograniczenia pamięci - zainstalowany w labie (opcja dla chętnych - można porównać z maximą!!!)
dokumentacja (w j. angielskim)

Zaliczenie labu - projekt sprawdzający czy Państwo jesteście w stanie przy pomocy Maximy czy Octava znaleźć rozwiązanie prostego zagadnienia początkowego czy brzegowego szczegóły za jakiś czas - max. 30% punktów z ćwiczeń. REszta punktów za kolokwium (pełniące rółnież rolę egzaminu połówkowego) - 40% i zadania domowe - 30% - zaliczenie od ok 50%.