wnum2

%O co tu chodzi? Mamy konkretne równanie do rozwišzania i konkretnš metodę
%której mamy użyć. Jest to silna analogia do labki o kwadraturach, bo
%równania różniczkowe rozwišzuje się poprzez całkowanie. Te metody sš
%iteracyjne tzn: Mamy dyskretnš dziedzinę czasu, zaczynamy od warunku
%poczštkowego który jest dany i po kolei wyznaczamy co ma się dziać w
%chwili t+h, t+2h, t+3h... t+n*h, gdzie h to nasza długoœć kroku
%całkowania.


clear all;
close all;
clc;


%Mamy nasze równanie w postaci dy/dt="coœ". To "coœ" wedle oznaczeń to jest
%nasze fn i tak należy je podstawić do naszych danych wzorów metod. Sam się
%trochę plštałem w tych oznaczeniach.

%Otwarta metoda Eulera (dy/dt=2yt/t^2+1)

y(1)=1;%Warunek poczštkowy. Pamiętamy że Matlab indeksuje od 1

h=0.1; %Krok calkowania
t=0:h:10; %Dziedzina czasu z krokiem co h

N=length(t);
for i=1:N-1
   y(i+1)=y(i)+h*2*(y(i)*t(i))/(t(i).^2+1); %yn+1=yn+h*fn, gdzie fn to nasze dyn/dtn
end

y_dok=(t.^2+1); %Nasze analityczne rozwišzanie dokładne, mamy podane, że y0=1

%Rysujemy
figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Otwarta metoda Eulera')
clear all; %Czyœcimy pamięć żeby nam się nie mieszało w następnych metodach

%Zamknięta metoda Eulera

%Czym się różni zamknięta metoda od otwartej? W zamkniętej po prawej
%stronie wzoru metody mamy człon z n+1. Tzn że nasze rozwišzanie dla
%kolejnej chwili czasowej niejako "zależy od samego siebie". Jest to
%delikatna trudnoœć tych metod, ale za to sš one chyba dokładniejsze. W
%każdym razie, za fn+1 podstawianmy naszš funkcję. Zatem z równania
%yn+1= yn+h*fn+1, gdzie  fn+1 = 2y(n+1)t(n+1) / t(n+1)^2+1, trzeba
%wyznaczyć na kartce yn+1. Przykład czegoœ takiego dla Adamsa-Multona r.2i3
%dołšczyłem do niniejszych rozwišzań.

y(1)=1;%Warunek poczštkowy

h=0.1; %Krok calkowania
t=0:h:10;
N=length(t);

y_dok=(t.^2+1);
for i=1:N-1
   y(i+1)=y(i)/(1-(2*h*t(i))/(t(i)^2+1)); %Wyznaczone yn+1 wstawiamy w naszš pętlę
end                                       %Nasza metoda wtedy ładnie zadziała

figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Zamknięta metoda Eulera')

clear all;

%Runge Kutty

y(1)=1;%Warunek poczštkowy

h=0.1; %Krok calkowania
t=0:h:10;
N=length(t);

%Tutaj cała trudnoœć polega na nie poplštaniu się w oznaczeniach.
%Fn to nasza funkcja, h/2 *(...) jest problematyczne.
%Drugi człon oznacza tyle, że do naszego fn zamiast tn podstawiamy
%tn+1, a zamiast yn podstawiamy yn+h*fn, gdzie fn to znowu nasza funkcja.
%Lekka incepcja i można się pogubić, ale przynajmniej nie trzeba z tego
%potem nic wyznaczać.

y_dok=(t.^2+1);
for i=1:N-1 %Wstaw sobie tš linijkę w wolframa i zobacz co mam na myœli, bo można oczoplšsu dostać
   y(i+1)=y(i)+ h/2*(2*y(i)*t(i)/(t(i)^2+1)+(2*t(i+1)*(y(i)+(2*h*y(i)*t(i))/(t(i)^2+1))/(t(i+1)^2+1)));
end

figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Runge Kutty')

%Gear r.2

%Wchodzimy w metody które wymagajš więcej niż jednego poprzednika. Tzn
%potrzebujemy dwóch lub więcej następujšcych po sobie punktów startowych,
%żeby to zadziałało. Tzw. metody wielokrokowe

%I tutaj zaczynajš się jaja. Byłem na konsultacjach u Opalskiego i
%dowiedziałem się, że dla podanych nam kroków całkowania(h=0.5,0.1,0.01,0.001)
%ta metoda Geara jest niestabilna, tzn pierwiastki równanania charekterystyczngo
%leżš poza obszarem stabilnoœci. I prawdopodobnie mieliœmy to zauważyć metodš
%prób i błędów, że skrócenie kroku całkowania generuje coraz większy błšd.

clear all;
h=1; %Krok calkowania dla którego to działa
y(1)=1;%Warunek poczštkowy
y(2)=y(1)*(h^2+1);%Drugi warunek poczštkowy uzależniony od długoœci naszego kroku
t=0:h:10;
N=length(t);

y_dok=(t.^2+1);
for i=2:N-1
   y(i+1)= ((4/3)*y(i)-(1/3)*y(i-1))/(1-(4/3)*t(i+1)/(t(i+1)^2+1));
end

figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Gear r.2')

%Adams-Bashforth r.2, któru również jest wielokrokowy i również œrednio
%działa dla zadanych nam długoœci kroków. Ewentualnie pomyliłem się w
%zapisie gdzieœ

clear all;
h=2; %Krok calkowania, dla którego to w sumie już nie ma sensu, ale wykres jest w miarę
y(1)=1;%Warunek poczštkowy
y(2)=y(1)*(h^2+1);
t=0:h:10;
N=length(t);

y_dok=(t.^2+1);
for i=2:N-1
    y(i+1) = y(i) + (h*3*((2*y(i)*t(i))/(t(i).^2+1))+ (2*y(i-1)*t(i-1))/(t(i-1).^2+1))/2;
end

figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Adams-Bashforth r.2')

%Adams-Multon r.2

clear all;
h=0.1;
y(1)=1;
t=0:h:10;
N=length(t);

y_dok=(t.^2+1);
for i=1:N-1 %Tutaj załšczyłem o co chodzi z tym wyznaczaniem yn+1, gdyby ktoœ miał wštpliwoœci
    y(i+1) =(y(i)+ (h*y(i)*t(i))/(t(i)^2+1)) / (1- h*t(i+1)/(t(i+1)^2+1));
end

figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Adams-Multon r.2')

%Adams-Multon r.3
%Ta metoda mimo że jest wielokrokowa to bardzo ładnie działa dla naszego
%przypadku nawet dla kroku 0.5
clear all;
h=0.5;
y(1)=1;
y(2)=y(1)*(h^2+1);
t=0:h:10;
N=length(t);

y_dok=(t.^2+1);
for i=2:N-1
    y(i+1)=(y(i)+(h/12)*(((16*y(i)*t(i))/(t(i)^2+1))-((2*y(i-1)*t(i-1))/(t(i-1)^2+1))))/(1-(10*h*t(i+1))/(12*(t(i+1)^2+1)));
end

figure
plot(t,y_dok,'-g')
hold on;
plot(t,y,'--r')
legend('dokładne','wyliczone')
title('Adams-Multon r.3')

%Długoœć kroku całkowania(h) możemy sobie zmieniać i patrzeć co się wtedy
%dzieje. Co z błędami? W zasadzie interesuje nas tylko błšd maksymalny
%danej metody, czyli największy jaki jest ona w stanie popełnić. Jeżeli
%nawet on jest w miarę niski, to jesteœmy szczęœliwi. Jak go policzyć?
%Mamy teraz w pamięci wartoœci z ostatniej metody więc robimy
max_err_am3=max(abs(y_dok-y))
%I mamy nasz błšd maksymalny. Teraz możnaby caaały ten długi kod wrzucić w
%czterokrotnš pętlę dla zadanych kroków całkowania, pozbierać te błędy w tablicę i
%wykreœlić ich zależnoœć od długoœci kroku przy pomocy loglog(), ale to już
%łatwe.




%Zostało jeszcze wyznaczenie kroku całkowania, dla którego Euler otwarty
%robi się niestabilny. O stabilnoœci jest w ksišżce w rozdziale 9.3.2.
%Euler otwarty jest stabilny, gdy lambdy leżš wewnštrz koła jednostkowego
%o œrodku (-1,j0). Nasze równanie jest w postaci y'=lambda*y.
%Nasza lambda=2t/(t^2+1), a warunek stabilnoœci to |h*lambda+1|<=1
%Widzimy że potęga mianownika jest wyższa niż licznika, więc wydaje mi się,
%że ta metoda będzie zawsze stabilna dla tego przykładu, ale mogę się mylić.



%No to teraz dla przećwiczenia można się pobawić z drugim przykładem czyli
%dy/dt=(4-2y)/3
%ale postępowanie jest identyczne więc już nie będę tego wszystkiego
%przepisywał po raz drugi :D