Numeryczne równania różniczkowe

• przeczytaj help dla lsode i narysuj wykres rozwiązania 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 są 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 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 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+k*h) w miejsce x_{n+k} itd w schemat x_{n+k}-\Phi(t+n*h,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)-h*f(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 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 i jego testami etc
  midpoint17.m - schemat midpoint (wersja 2017/18)
  midpoint.m - implementacja schematu midpoint - stara wersja
  testmidpoint.m - testy rzędu zbieżności schematu midpoint oraz brak stabilności dla dx/dt=-x z x(0)=1 na długim odcinku czasu (im mniejsze h tym później) - wykorzystuje midpoint.m
  taylor17.m - schemat Taylora
  testlte.m - testy rzędu lokalnego błędu (LTE) dola Eulerów i midpoint

• Napisz m-plik ze schematem Adamsa-Bashforda 2giego rzędu -
  x(n+1)=x(n)+(h/2)*[3*f(n)-f(n-1)] z f(n)=f(t(n),x(n))
  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ść trajektorii? Powtórz zadanie dla y'=-sin(y) z y(0)=0 y'(0)=1.
• Zaimplementuj schemat predyktor korektor - predyktorem jawnym Eulerem, korektorem niejawnym Eulerem - rozwiązując równania nieliniowe metoda iteracji prostych tzn x^{k+1}=x_n+h*f(x^k,t^k+h) (jeśli x^{k+1}=x^k to x^k=x_{n+1} dla niejawnego Eulera) za wziąć x^0=x_n+hf_n (predyktor jawny Euler) - napisz ten schemat biorąc ilość nieliniowych iteracji M jako opcje - testuj dla M=1,2,3 - można napisać tez wersje z warunkiem stopu zależnym od błędu np zatrzymujemy się jak
  |x^{k+1}-x_n-h*f(x^k,t^k+h)|<= h lub <= x_n*h
• Napisz 2krokowy schemat Adamsa-Moultona
  x(n+1)=x(n)+(h/12)*(5*f(n+1)+8*f(n)-f(n-1)) z f(n)=f(t(n),x(n))
  i przetestuj rząd zbieżności - do policzenia x(1) weź schemat Heuna.
• Zaimplementuj 2krokowy schemat predyktor - korektor: za predyktor 2krokowy Adams-Bashford i za korektor 2-krokowy Adams-Moulton. Przetestuj rząd zbieżności.
• Powtórz testy rzędu lokalnego błędu i zbieżności 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) można brać dokładne wartości rozwiązania lub obliczyć je schematem Heuna (rzędu 2) - czy wyniki w obu przypadkach będą różne jakościowo?
  

  Skrypty

  adamsb2s17.m - explicit Adams-Bashforth scheme (of order 2)
  testadamsb2s17.m - tests of explicit Adams-Bashforth scheme (of order 2)
  AdamsB2step.m - jawny Adams-Bashforth 2 krokowy
  Taylor.m - jawny Taylor (rząd 2)
  modEuler.m - jawny zmodyfikowany Euler (Runge-Kutta rząd 2)
  testAB2start.m - testy wartości startowe (x^h_1) dla 2krokowego Adams-Bashforth
  testTaylor.m - testy - Taylor (rząd 2)
  testmodEul.m - testy - jawny zmodyfikowany Euler (Runge-Kutta rząd 2)
  LTEtests.m - testy rzędu lokalnego błędu kilku schematów np schematów Eulera

• Porównaj 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 użyciem lsode z opcją stiff i raz z non-stiff (lsode_options("integration method","non-stiff")). Czas obliczeń uzyskamy za pomocą: (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.
  5. -y''+(1+t*t)y=sin(x)+ (1+x*x)*sin(x) z y(0)=0 i y(1)=sin(1) (rozwiązanie y(x)=sin(x)).
  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().).
• 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 punktów. 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 dyskretnych 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 błąd z rozwiązanie m dokładnym (które można policzyć - jest postaci u(x)=cexp(t)+dexp(-t) a współczynniki spełniają układ r. liniowych c+d=1;c+d*exp(-2b)=exp(-b)) Porównaj z metoda strzałó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} .

1. Policz rząd 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 był warunek typu Neumanna to bierzemy różnice w tył). 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). Przetestuj teraz dla rozwiązania takiego, że u''(a)=0 np u=sin(x) z a=0.
2. Napisz funkcje rozwiązująca zadanie - eps*u''+u'=f z u(a)=\alpha,u(b)=\beta - eps stała >0 modyfikując kod z poprzednich zadań. Pochodna u' w schemacie MRS aproksymujemy różnica 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 - błąd 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. Rząd schematu. Policz eksperymentalnie rząd lokalnego błędu dla dyskretyzacji MRS -\triangle u=f w [0,1]^2 , u=0 na brzegu dla znanego rozwiązania u(x)=\sin(\pi*x)\sin(\pi*y) w normach dyskretnych L^2 i maksimum (można narysować wykres błędu punktowego L_hu-f_h)
6. Klasyczny problem - rząd błędu 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. Wpływ współczynnika 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 środka kwadratu a zero w pozostałych punktach (przybliżenie delty Diraca w (0,0)). na siatce N punktów dla N=40,80,160 pktów w każdym kierunku dla różnych c=0,1,10,100,1000 . Narysuj rozwiązanie MRS. Czy widać większe 'rozmycie' rozwiązania dla c=0?
8. Rząd błędu dla mało regularnego rozwiązania 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. Niedokładny 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. Stabilność . Policz normy indukowane przez dyskretne L^1_h, L^2_h i max - macierzy A+cI i jej odwrotnej dla różnych c (A macierz dyskretyzacji MRS 2wymiarowego Laplacianu - na siatce rownomiernej na kwadracie) dla N=10,20,40,80,160.- UWAGA! jak to zrobić z wykorzystaniem funkcji octave'a norm()?

• Rozwiąż danym skryptem (lub pisząc swój) u_t-u_{xx}=f z zerowymi war brzegowymi Dirichelta i war początkowym u0(x)=sin(x), peak tzn funkcja równa 100 w 1 pkcie siatki a zero w pozostałych (czy rozwiązanie się wygładzi dla t>0?)
• zmień kod tak by działał dla dowolnych niezerowych warunków 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.
• zmień kod tak by działał 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)
• zmień kod by działał dla dyskretyzacji MESem u_t-u_{xx}=f(t,x) z warunkami Dirichleta na siatce dowolnej (używając np macierzy ze skryptów z poprzedniego labu)
• zaimplementuj jawny, niejawny schemat Eulera, Cranka-Nicholson dla u_t-u_{xx}=f u(t,a)=ga(a) u(t,b)=gb(b) u(0,x)=u0(x) Testuj dla różnych wartości parametrów dyskretyzacji np ustalając h i biorąc małe lub duży krok po czasie np dla f=0 and ga=gb=0. Testuj rząd zbieżności metodą połowienia kroku względem tylko h lub \tau. (np. dla u(t)=exp(-t)sin(x) a=0, b=\pi)

Skrypty m-pliki octave'a