Numeryczne równania różniczkowe

semestr zimowy 2014-15

Egzamin ustny

1. równań różniczkowych zwyczajnych
2. równań różniczkowych eliptycznych
3. równań różniczkowych ewolucyjnych (paraboliczne i hiperboliczne pierwszego rzędu)

1. schematy dla równań zwyczajnych jednokrokowe i wielokrokowe
2. metodę różnic skończonych
3. metodę elementu skończonego

Skrypt

Literatura

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

W języku angielskim

Podstawowe podręczniki

1. Deuflhard, Peter, Bornemann, Folkmar, Scientific Computing with Ordinary Differential Equations, Texts in Applied Mathematics, Vol. 42, Springer-Verlag, New York, 2002. (teoria RRZ, schematy RRZ, zadania brzegowe w 1 wymiarze) Z serwerów MIMUW mozna sciagnac plik pdf (aktualne XII 2014): Springer link
2. David F. Griffiths, Desmond J. Higham, Numerical Methods for Ordinary Differential Equations, Springer-Verlag, 1st Edition, London, 2010. Schematy do RRZ - prosty podrecznik. Z serwerów MIMUW mozna sciagnac plik pdf (aktualne XII 2014): Springer link
3. Claes Johnson, Numerical solution of partial differential equations by the finite element method, Cambridge University Press, Cambridge, 1987.
4. Randall J. LeVeque, Finite difference methods for ordinary and partial differential equations, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007, Steady-state and time-dependent problems. (schematy dla RRZ, metoda różnic skończonych (MRS) dla równań eliptycznych i parabolicznych)
5. Alfio Quarteroni, Riccardo Sacco, and Fausto Saleri, Numerical mathematics, Texts in Applied Mathematics, vol. 37, Springer-Verlag, New York, 2000. (metody dla RRZ, teoria dla MRS dla RRCZ i elementy MESu). Z serwerów MIMUW mozna sciagnac plik pdf (aktualne XII 2014): Springer Link
6. John C. Strikwerda, Finite difference schemes and partial differential equations, second ed., Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2004. (MRS dla RRCZ - wszystkie typy równań)

Monografie i bardziej zaawansowane podręczniki

1. Dietrich Braess, Finite elements, third ed., Cambridge University Press, Cambridge, 2007, Theory, fast solvers, and applications in elasticity theory, Translated from the German by Larry L. Schumaker. (zaawasowany podręcznik - godny polecenia)
2. Susanne C. Brenner and L. Ridgway Scott, The mathematical theory of finite element methods, third ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008. (bardzo zaawansowany podręcznik - właściwie monografia)
3. J. C. Butcher, Numerical methods for ordinary differential equations, second ed., John Wiley and Sons Ltd., Chichester, 2008.
4. P. G. Ciarlet and J.-L. Lions (eds.), Handbook of numerical analysis. Vol. II, Handbook of Numerical Analysis, II, North-Holland, Amsterdam, 1991, Finite element methods. Part 1.
5. Philippe G. Ciarlet, The finite element method for elliptic problems, Classics in Applied Mathematics, vol. 40, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002, Reprint of the 1978 original [North-Holland, Amsterdam].
6. E. Hairer, S. P. Norsett, and G. Wanner, Solving ordinary differential equations. I, second ed., Springer Series in Computational Mathematics, vol. 8, Springer-Verlag, Berlin, 1993, Nonstiff problems.
7. E. Hairer and G. Wanner, Solving ordinary differential equations. II, second ed., Springer Series in Computational Mathematics, vol. 14, Springer-Verlag, Berlin, 1996, Stiff and differential-algebraic problems.
8. Alfio Quarteroni and Alberto Valli, Numerical approximation of partial differential equations, Springer Series in Computational Mathematics, vol. 23, Springer-Verlag, Berlin, 1994. Z serwerów MIMUW mozna sciagnac plik pdf (aktualne XII 2014): Springer Link

Lab

Program labów

• Lab 1 (21 X 2014) - wstępne zapoznanie się z octavem - octave jako kalkulator naukowy, operator dwukropek : tworzenie macierzy, 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. 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 @). Schematy Eulera i inne proste (jak zdazymy).
  • Zapoznaj się z kodem w nrrbasic.m - prostym skryptem z przykładami.
  • 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.
  • Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo - czyli bez użycia pętli).
  • 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 sa 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.
    Tu zamieszczam m-plik ze schematem otwartym Eulera i skrypt z jego testami:
    exEuler.m - schemat otwarty Eulera (dziala w wielu wymiarach)
    testexEul.m - proste testy schematu Eulera (skalarne i 2-wymiarowe)
• Lab 2 - Funkcja lsode() octave'a i proste schematy ciąg dalszy. 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 otwrtej tego samego rzędu np rzędu 2 np Taylora dla metody midpoint. Eksperymentalne badanie rzędu.
  • przeczytaj help dla lsode i narysuj wykres rozwiążania dla y'=-0.1y y(0)=1 na [0,20] za pomocą tej funkcji. Powtórz zadanie dla y'=(cos(x))^2 z y(0)=1 oraz y'=1+y*y z y(0)=0(jakie sa rozwiązania?) na [0,T] dla T=1,10,100.
  • 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 sa 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().
  • Przetestu 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 błąd lokalny schematu 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}
  • Przetestu eksperymentalnie rzędy 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 schematem midpoint, Taylor etc
    OrderExEuler.m - testy rzędu (lokalnego błędu) schematu otwartego Eulera
    midpoint.m - schemat midpoint <
    testmidpoint.m - testy schematu midpoint rząd zbieżności i niestabilnośc dla dużych T (dx/dt=-x z x(0)=1)
    taylor.m - schemat Taylor
    testtaylor.m - testy schematu Taylora (w 2 wymiarach)
• Lab 3 testowanie rzędu schematów i rzędu zbieżności. Start dla schematów wielokrokowych. Metoda strzałów.
  • Napisz m-plik ze schematem Adamsa-Bashforda 2giego 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śc trajektorii? Powtórz zadanie dla y'=-sin(y) z y(0)=0 y'(0)=1.
  • 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 jakos wymaga uzyskania s_k kolejnych iteracji fsolve().).
  
  AdamsB2.m - explicit Adams-Bashforth scheme (of order 2)
  testAB2start.m - testy wartosci startowej (x^h_1) dla explicit Adams-Bashforth scheme (of order 2)
  testAdamsBo2.m - testy dla schematu otwartego Adams-Bashforth (rzedu 2) jak dla midpoint - rzad zbieznosci etc
  testshoot.m -metoda strzalow dla zadania brzegowego: -y''+y=0; y(0)=1 y(b)=1 (b=1 - m. st. dziala OK, b=20 - nie dziala - dlaczego?)
  linshooting.m -m-plik z funkcja linshooting() rozwiazujaca: -y''+d(x)y+c(x)y=f(x); y(a)=ya y(b)=yb przy uzyciu metody strzalow
  shooting.m -m-plik z funkcja shooting() rozwiazujaca: -y''=F(t,y,y'); y(a)=ya y(b)=yb przy uzyciu metody strzalow (uzywajac lsode() i fsolve())
• Lab 4 - MRS dla -u''+cu=f z warunkami Dirichleta : u(a)=alpha u(b)=beta i mieszanymi: u'(a)=alpha, u(b)=beta
  • Utwórz 3diagonalną macierz z 2 na diagonali i -1 na pod i nad diagonalach z użyciem funkcji octave'a diag() i jako macierz rzadka bez tworzenia macierzy w formacie pełnym.
  • Rozwiąż zadanie dyskretne powstałe z dyskretyzacji MRS zadania: -u''=f , u(0)=u(\pi)=0 na [0,\pi] z rozwiązaniem u(x)=\sin(x) na siatce 10,20,40,80 pktow. Policz błąd dyskretny w dyskretnych normach max i L^2.
  • Policz rząd standardowego schematu MRS dla -u''=\sin(x) , u(0)=u(\pi)=0 w normach dyskretncyh L^2 i max.
  • Rozwiąż zadanie dyskretne powstałe z dyskretyzacji MRS zadania: -u''+u=0 , u(0)=u(b)=1 na [0,b] dla b=1,2,,4,8,16 etc na siatkach 20 i 100 pkt. Porównaj blad z rozwiazanie mdokladnym (ktore mozna policzyc - jest postaci u(x)=cexp(t)+dexp(-t) a wspolczynniki spelniaja uklad r. liniowych c+d=1;c+d*exp(-2b)=exp(-b)) Porównaj z metoda strzalów.
  • Policz rząd zbieżności MRS dla -u''=\sin(x) , u(0)=u(\pi)=0 w normach dyskretnych L^2 i max tzn. liczymy błędy e_h i e_{h/2} i potem e_h/e_{h/2} .
  • Powtórz poprzednie zadania ale dla dyskretyzacji MRS nowego zadania brzegowego: -u''=f na [0,1] , u'(0)=\alpha;u(1)=\beta bierzemy f , \alpha,\beta dla znanego rozwiązania u(x)=\sin(x+1). (unikamy takiego dla którego znikają niektóre pochodne w x=0 - mogłoby to sztucznie podwyższyć rząd schematu) Warunek Neumanna aproksymujemy różnicą wprzód (gdyby w b byl warunek typu Neumanna mana pochodna to bierzemy roznice wtyl). Policz eksperymentalnie rząd schematu, rząd zbieżności w normach dyskretnych itd Narysuj wykres błędu czy błąd wyraźnie wyższy blisko któregoś z końców? (Można się spodziewać, że blisko lewego końca błąd będzie wyższy). Przetetuje teraz dla rozwiazania takiego, ze u''(a)=0 np u=sin(x) z a=0.
  • Napisz funkcje rozwiazaujaca zadanie - eps*u''+u'=f z u(a)=\alpha,u(b)=\beta - eps stala >0 modyfikujac kod z poprzednich zadan. Pochodna u' w schemacie MRS aproksymujemy roznica centralna tzn schemat dla x_k=a+k*h z h=(b-a)/N i dla 0 < k< N to
    u_0=\alpha; u_N=\beta,
    eps(-u_{k-1}+2u_k-u_{k+1})/h^2 + (u_{k+1}-u_{k-1}/(2h)=f(x_k)
    Przetestuj dla a=0,b=1, f=10*\max(sin(5x-1),0)i eps=10^{-k} k=0,1,2,3,4 i ustalonego N=50 etc
  
  FDsolver.m funkcja rozwiązująca -u''+cu=f z warunkami Dirichleta : u(a)=alpha, u(b)=beta (MRS rzędu dwa)
  testFD.m testy MRS dla -u''+cu=f z warunkiem Dirichlea : u(a)=alpha u(b)=beta -zad 1 u''=u u(0)=u(b)=1 b=1,4,8,16 etc wykresy;zad 2 rzad standardowej MRS dla -u''=sin(x) u(0)=0 u(b)=sin(b); zas 3 - jak zad 2 ale z niedokladnym warunkiem brzegowym u_h(b)=sin(b+h/2) - order 1
  FDmixlft.m funkcja rozwiązująca -u''+cu=f z warunkami mieszanymi : u'(a)=alpha, u(b)=beta (MRS rzędu jeden)
• Lab 5 FDM for Poisson equation in 2D
  1. Macierz dyskretnego Laplacianu na kwadracie - siatko równomierna Utwórz macierz 5-diagonalną symetryczną - dyskretny Laplacian w 2 wymiarach na kwadracie na siatce jednorodnej. Wsk: taka macierz to T_N\tensorpoduct I + I \tensorproduct T_N dla T_N macierzy 1-wym Laplacianu; funkcja octave'a korn(A,B) tworzy produkt tensorowy 2 macierzy
  2. Klasyczny problem - blad Rozwiąż MRS $ -Laplacian u=f$ w $[0,1]^2$, u=0 na brzegu ze znanym rozwiązaniem $u(x)=\sin(\pi*x)\sin(\pi*y)$ na siatce 10,20,40,80 pktów w każdym kierunku. Policz błędy dyskretne w normach L^2 i max. Narysuj błąd i rozwiązanie MRS (funkcja octave'a mesh()).
  3. Rzad schematu. Policz eksperymentalnie rzad zbieznosci dla dyskretyzacji MRS $-\triangle u=f$ w $[0,1]^2$, $u=0$ na brzegu dla znanego rozwiazania $u(x)=\sin(\pi*x)\sin(\pi*y)$ w normach dyskretnych $L^2$ i maksimum (mozna narysowac wykres bledu punktowego L_hu-f_h)
  4. Klasyczny problem - rzad bledu Policz eksperymentalnie rząd zbieżności schematu MRS dla $-Laplacian u=f$ in $[0,1]^2$, $u=0$ na brzegu ze znanym rozwiązaniem $u(x)=\sin(\pi*x)\sin(\pi*y)$ w normach dyskretnych $L^2$ i max metodą połowionego kroku. I.e liczymy $e_h$ i $e_{h/2}$ a potem: $e_h/e_{h/2}$ .
  5. Wplyw wspolczynnika dyfuzji (1/c) - prawa strona aproksymacja delty Diraca Rozwiąż MRS $ - Laplacian u + cu =f$ w $[-1,1]^2$, u=0 na brzegu z $f$ równego N^3 w ustalonym punkcie siatki bliskim srodka kwadratu a zero w pozostalych punktach (przyblizenie delty Diraca w (0,0)). na siatce N punktów dla N=40,80,160 pktów w każdym kierunku dla róznych $c=0,1,10,100,1000$. Narysuj rozwiązanie MRS. Czy widac wieksze 'rozmycie' rozwiazania dla c=0?
  6. Rzad bledu dla malo regularnego rozwiazania Policz eksperymentalnie rząd zbieżności schematu MRS dla $-Laplacian u=f$ w $[-1,1]^2$, $u=0$ na brzegu ze znanym rozwiązaniem u(x) które ma niską regularność np. $C^2$ lub $C^3$ - weź np. $u(x,y)=\sin(\pi*x)*f(y)$ dla f(x) kawałkami wielomianu kubicznego $[-1,0]$ na $[0,1]$ takiego, że $f(-1)=f(1)=0$ i $f,f',f''$ ciągłe w $x=0$. (interpolacja Hermite'a). Weź siatki zawierające punkty na prostej y=0 i takie że na tej prostej nie ma punktów siatki.
  7. Niedokladny warunek brzegowy Dirichleta - punkty brzeg. siatki NIE na brzegu obszaru Policz eksperymentalnie rząd zbieżności schematu MRS $-Laplacian u=f$ in $[-1,1]^2$, $u=0$ na brzegu ze znanym rozwiązaniem $u(x)=\sin(\pi*x)\sin(\pi*y)$ w dyskretnych normach L^2 i max dla siatek które NIE MAJĄ PUNKTÓW na brzegu oryginalnego obszaru np. i wprowadź warunki brzegowe zerowe dla punktów z k,l = 0,N+1 (czyli najbliższych brzegowi obszaru wyjściowego).
  8. Stabilnosc . Policz normy indukowaneprzez dyskretne L^1_h, L^2_h i max - macierzy A+cI i jej odwrotnej dla roznych c (A macierz dyskretyzacji MRS 2wymiarowego Laplacianu - na siatce rownomiernej na kwadracie) dla N=10,20,40,80,160.- UWAGA! jak to zrobic z wykorzystaniem funkcji octave'a norm()?
  
  FDsolver2d.m -funkcja rozwiazujaca zadanie -Laplacian u = f w (a,b)^2 u=g na brzegu metoda roznic skonczonych
  testFDM2d.m - testy rzędu zbieżności dla -Laplacian u=f na [0,1]^2 z u=g na brzegu; - 5-punktowy stensyl na siatce jednorodnej dla rozwiązania u = sin(pi*x_1)*sin(pi*x_2), (wtedy f=2*(pi)^2*u; u=0 na brzegu)
• Lab 6 MES dla $-u'' +cu=f$ on $[a,b]$ warunki Dirichleta. Wszedzie uzywamy standardowa przestrzen MESu funkcji ciaglych kawalkami liniowych
  • Siatka jednorodna;war. brzegowe Dirichleta jednorodne - zerowe Policz rozwiazanie MES $-u'' +cu=f$ na $[a,b]$dla $u=\sin(\pi*x)$ i $a=0;b=\pi$ uzywajac kwadratury trapezow dla wektora prawej strony dla $N=10,20,40,80,160$. Policz normy $L^2$, dyskretna $\infty$ i $H^1$ bledu $u_h-I_h u$. przetestuj rzad zbieznosci.
  • Siatka niejednorodna;war. brzegowe Dirichleta jednorodne - zerowe Policz rozwiazanie MES $-u'' +cu=f$ na $[a,b]$ dla $u=\sin(\pi*x)$ i $a=0;b=\pi$ uzywajac siatki z krokiem $2h$ na [a,(b+a)/2] i krokiem $h$ na [(a+b)/2,b] oraz kwadratury trapezow dla wektora prawej strony dla $N=10,20,40,80,160$. Policz normy $L^2$, dyskretna $\infty$ i $H^1$ bledu $u_h-I_h u$. przetestuj rzad zbieznosci.
  • Siatka niejednorodna;war. brzegowe Dirichleta jednorodne - zerowe Policz rozwiazanie MES $-u'' +cu=f$ na $[a,b]$ dla $u=\sin(\pi*x)$ i $a=0;b=\pi$ uzywajac nieregularnej siatki postaci: x_0=a;$\{x_{n+1}\}\cup\{b\}$ dla $x_{n+1}=x_n+h_n\leq b$ z $h_n=(b-a)/(sqrt(n)*N)$ oraz kwadratury trapezow dla wektora prawej strony dla $N=10,20,40,80,160$. Policz normy $L^2$, dyskretna $\infty$ i $H^1$ bledu $u_h-I_h u$. przetestuj rzad zbieznosci podwajajac $N$
  • Siatka jednorodna;war. brzegowe Dirichleta niejednorodne - niezerowe Policz rozwiazanie MES $-u'' +cu=f$ na $[a,b]$ dla $u=\sin(\pi*x)$ i $a=0;b=1$ uzywajac kwadratury trapezow dla wektora prawej strony dla $N=10,20,40,80,160$. Policz normy $L^2$, dyskretna $\infty$ i $H^1$ bledu $u_h-I_h u$. przetestuj rzad zbieznosci.
  • Siatka jednorodna;war. brzegowe mieszane jednorodne - zerowe Policz rozwiazanie MES $-u'' +cu=f$ na $[a,b]$ z $u(a)=0;u'(b)=0$ dla $u=\sin(x)$ i $a=0;b=1$ uzywajac kwadratury trapezow dla wektora prawej strony dla $N=10,20,40,80,160$. Policz normy $L^2$, dyskretna $\infty$ i $H^1$ bledu $u_h-I_h u$. przetestuj rzad zbieznosci.
  • Ogólny elliptic operator eliptyczny; siatka jednorodna; jednorodne warunki brzegowe D. Policz rozwiazanie MES $-(a(t)u')'+c(t)u=f(t)$ $u(le)=u(re)=0$ na siatce jednorodnej. Macierz sztywnoci i wektor prawej strony policz uzywajac kwadratury trapezow. rule. Przetestuj zbieznosc dla $u=\sin(t)$, $a=1+t^2$, $c=0$ na $[0,\pi]$. Powtorz obliczenia dla $c(t)=1+t$.
  • Siatka jednorodna;war. brzegowe Dirichleta jednorodne - zerowe; dokladniejsze obliczanie prawej strony Policz rozwiazanie MES $-u'' +cu=f$ na $[a,b]$dla $u=\sin(\pi*x)$ i $a=0;b=\pi$ uzywajac kwadratury SIMPSONA dla wektora prawej strony dla $N=10,20,40,80,160$. Policz normy $L^2$, dyskretna $\infty$ i $H^1$ bledu $u_h-I_h u$. przetestuj rzad zbieznosci. Porównaj z przypadkiem gdy prawa strone obliczalismy kwadratura trapezow.
  
  FEM1Dsolver.m - solver MES w 1 wymiarze dla -u''+cu=f z war brzegowymi Dirichleta - na siatce dowolnej
  testMESDbc.m funckja w m-pliku z testami rzedu MESu liniowego na roznych siatkach dla znanych rozwiazan
  testMESor.m skrypt wywolujacy funkcje octave'a testMESDbc() (z pliku testMESDbc.m - trzeba go sciagnac) - testy MESu liniowego na roznych siatkach dla znanych rozwiazan (malo skomentowane trzeba zajrzec do kodu by zrozumiec o co chodzi - ogolnie macierze ratio/rratio zawieraja stosunki bledu dla siatki x podzielone przez blad dla siatki powstalej z zageszczenia x poprzez podzial kazdego odcinka na dwa pododcinki - ratio zawiera stosunki odpowiednich bledow dyskretnych - tylko na siatce na ktorej szukamy rozwiazan - rratio odpowiednie stosunki bledow na duzo gestszej siatce)

Skrypty m-pliki octave'a