Numeryczne równania różniczkowe

• 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 midpoin i jego testami etc
  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)

• 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 Taylora, 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 schemat predyktor korektor - predykotem jawnym Eulerem, korektorem niejawnym Eulerem - rozwiazujac rownania nieliniowe metoda iteracji prostych tzn x^{k+1}=x_n+h*f(x^k,t^k+h) (jesli x^{k+1}=x^k to x^k=x_{n+1} n+1 krk niejawnego Eulera) za x^0=x_n+hf_n (predyktor jawny Euler) - napisz ten schemat biorac ilosc nieliniowych iteracji M jako opcje - testuj dla M=1,2,3 - mozna napisac tez wersje z warunkiem stopu zaleznym od bledu np zatrzymujemy sie jak |x^{k+1}=x_n+h*f(x^k,t^k+h)|\leq h lub \leq x_n*h
• Powtórz testy rzedu lokalnego bledu i zbieznosci dla jawnego 3 krokowego schematu Adamsa: x(n+1)=x(n)+(h/12)*[23f(n)-16f(n-1)+5f(n-2)] z f(n)=f(t(n),x(n)) Za x(1),x(2) mozna brac dokladne wartosci rozwiazania lub obliczyc je schematem Heuna (rz 2) - czy wyniki w obu przypadkach beda rozne jakosciowo?
  AdamsB2step.m - jawny Adams-Bashforth 2 krokowy
  Taylor.m - jawny Taylor (rzad 2)
  modEuler.m - jawny zmodyfikowany Euler (Runge-Kutta rzad 2)
  testAB2start.m - testy wartosci startowe (x^h_1) dla 2krokowego Adams-Bashforth
  testAdamsB2st.m -testy wrzedu etc dla 2krokowego Adams-Bashforth
  testTaylor.m - testy - Taylor (rzad 2)
  testmodEul.m - testy - jawny zmodyfikowany Euler (Runge-Kutta rzad 2)
  LTEtests.m - testy rzedu lokalnego bledu kilku schematów np schematów Eulera
• Lab 4 (22 i 29 listopada 2015) - MRS dla -u''+cu=f z warunkami Dirichleta : u(a)=alpha u(b)=beta i mieszanymi: u'(a)=alpha, u(b)=beta
  • Porownaj czas obliczania X(T) dla X rozw. dX/dt=AX z X(0)=[1;2] i A=[-1,-10;-1,-11] i T=10,100,1000,1000 etc raz z uzyciem lsode z opcja stiff i raz z non-stiff (lsode_options("integration method","non-stiff")). Czas obliczen uzyskamy za pomoca: (tic;lsode(..);toc)
  • 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().).
  • 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} .
  
  linshoot16.m -m-plik z funkcja linshoot15() rozwiazujaca: -y''+p(x)y+q(x)y=f(x); y(a)=alpha y(b)=beta przy uzyciu metody strzalow
  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?)
  shootexample16.m -skrypt : rozwiazujacy: y''=y^2; y(0)=1 y(1)=-4; y(a)=ya y(b)=yb przy uzyciu metody strzalow (uzywajac lsode() i fsolve()) - mamy 2 rózne rozwiazania problemu
  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())
  fdmdir16.m funkcja rozwiazujaca MRS zadanie -u''+c*u=f(t) z : u(a)=alpha, u(b)=beta; c -nieujemna stala
  fdmdirtest16.m funkcja testujaca rzad zbieznosci MRS w dyskretnych normach dla zadanie -u''+c*u=f(t) z : u(a)=alpha, u(b)=beta (dla znanych rozwiazan - gladkich i niekoniecznie gladkich); c -nieujemna stala
• Lab 5 MRS w 1wym cd i w 2 wymiarach
  1. Policz rzad zbieznosci 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 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.
  2. Napisz funkcje rozwiazujaca 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
  3. 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 kron(A,B) tworzy produkt tensorowy 2 macierzy
  4. 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()).
  5. Rzad schematu. Policz eksperymentalnie rzad lokalnego bledu 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)
  6. Klasyczny problem - rzad bledu zbieznosci 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. tzn liczymy $e_h$ i $e_{h/2}$ a potem: $e_h/e_{h/2}$ .
  7. 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?
  8. 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.
  9. 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).
  10. 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()?
  
  fdmix16.m - funkcja dla -u''+c*u=f in (a,b) u'(a)=al u(b)=be z uzyciem MRS rzedu jeden (wb Neumanna - aproks. roznica wprzod w a) lub schematem podwyzszonego rzedu (z uzyciem rownania)
  fdmixtest16.m - funkcja z testami dla -u''+c*u=f in (a,b) u'(a)=al u(b)=be z uzyciem solvera MRS rzedu jeden (wb Neumanna - aproks. roznica wprzod w a) lub schematem podwyzszonego rzedu (z uzyciem rownania)
  fddirng16.m - funkcja dla -u''+c*u=f in (a,b) u(a)=al u(b)=be z uzyciem MRS rzedu jeden z ostatnim prawym pkt siatki x_N < b (b=x_N+0.5*h) - wb Dirichleta przyblizone bezposrednio u_N=be lub przez aproksyamcje Collatza i.e. ostatnie rownanie MRS l(b)=be z l(t) liniowym t.z. l(x_k)=u_k k=N,N-1
  fdmdirngtest16.m - - funkcja z testami dla -u''+c*u=f in (a,b) u'(a)=al u(b)=be z uzyciem solvera MRS - solver jak w fdmdirng16.m
  fdmdir2D16.m - funckja dla -Laplacian u+c*u=f in (a,b)x(a,b) - zerowe wb Dirichelta - MRS na siatce rownomiernej - krok ten sam w obu kierunkach- 5 pkt stensyl
• Lab 7 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.
  
  FEM1Dmats.m - macierze std liniowego MES w 1 wymiarze dla -u''+u=f - na siatce dowolnej tzn. macierz sztwynosci A=(\int_a^b \phi_k'\phi_l')_{k,l} i masy B=(\int_a^b \phi_k\phi_l)_{k,l} i \phi_k nodalna f bazowa zwiazana z x_k
  FEM1dDirSol16.m - liniowy MES solver dla -au''+cu=f z war brzeg. Dirichleta - dowolna siatka
  test1dfem16.m kilka testów - liniowy MES solver dla -au''+cu=f z war brzeg. Dirichleta - dowolna siatka
• Lab 8 MRS dla r. parabolicznego tj. u_t-u_{xx}=f zdyskretyzowane MRS wzgledem x - then uzyjmy lsode()
  • Rozwiaz danym skryptem (lub piszac swoj) u_t-u_{xx}=f z zerowymi war brzegowymi Dirichelta i war poczatkowym u0(x)=sin(x), peak tzn funkcja rowna 100 w 1 pkcie siatki a zero w pozostalych (czy rozwiazanie sie wygladzi dla t>0?)
  • zmien kod tak by dzialal dla doowlnych niezerowych warunkow brzegowych tzn u_t-u_{xx}=f u(t,a)=ga(t) u(t,b)=gb(t) u(0,x)=u0(x) dla dowolnych f, ga,gb,u0.
  • zmien kod tak by dzialal dla niezerowych war brzegowych Neumanna tzn u'(t,a)=al(t);u'(t,b)=be(t) lub mieszanych np u(t,a)=al(t);u'(t,b)=be(t)
  • zmien kod by dziala dla dyskretyzacji MESem u_t-u_{xx}=f(t,x) z warunakmi Dirichleta na siatce dowolnej (uzywajac np macierzy ze skryptow z poprzedniego labu)
  • zaimplementuj jawny, niejawny schemat Eulera, Cranka-Nicholsan dla u_t-u_{xx}=f u(t,a)=ga(a) u(t,b)=gb(b) u(0,x)=u0(x) Testuj dla roznych wartosci parametrow dyskretyzacji np ustalajac h i biorac male lub duzy krok po czasie np dla f=0 and ga=gb=0. Testuj rzad zbieznosci metoda polowienia kroku wgzledem tylko h lub \tau. (np. dla u(t)=exp(-t)sin(x) a=0, b=\pi)
  
  FD1dtest.m -testy MRS dla u_t-u_{xx}=0 u(0)=u(pi)=0 u(0,x)=sin(x) - u(t,x)=exp(-t)sin(x) rozwiazanie
  FD1dtest16.m -testy MRS dla u_t-u_{xx}=f(t,x) u(a)=al(t) u(b)=be(t) u(0,x)=sin(x) - u(t,x)=exp(-t)sin(x) rozwiazanie

Skrypty m-pliki octave'a