wnum

clear all;
close all;
clc;
%O co tu chodzi? Matlab ma zbiór rowziązywaczy układów równań
%różniczkowych. Tzw funkcji ODE. Różnią się one trochę, jedne działają
%wolniej a dokładniej, inne odwrotnie. My zajmiemy się ode45 i ode15s. Ta
%pierwsza jest właśnie wolniejsza ale dokładniejsza, ta druga szybsza, ale
%mniej dokładna. Taka funkcja ode zwraca nam dyskrętną dziedzinę czasu
%którą sobie dobrała dla danego przypadku(bo jest na tyle mądra) oraz
%oczywiście gotowe rozwiązanie.

tspan=[0 30*1e-3]; %Potrzebujemy przedziału czasowego
u0=[0;1]; %I wartości początkowych dla naszej metody ode

[t_45,u_45]=ode45('odefun',tspan,u0); %A co to odefun? Wyjasniłem w m-pliku o tej nazwie.

%t_45 to nasza dziedzina czasu, a u_45 to macierz, która zawiera dwie kolumny tj. jedną i
%druga funkcję której szukaliœmy czyli u1(t) i u2(t). Nazwy zmiennych dla
%zachowania porządku.

[t_15s,u_15s]=ode15s('odefun',tspan,u0);

%No dobra, zatem mamy już nasze policzone rozwiązania przybliżone, a co z
%dokładnym? No cóż, w przypadku układu równań sprawa nie jest taka prosta.
%Jak spojrzymy na to jak każą nam policzyć u2 dokładnie, to można się
%zagubić. Nie umiem szczególnie układów równań rózniczkowych, ale najwyrażniej
%jakoś da się to zrobić, żeby wyliczyć stałe c1 i c2, od których nasze
%rozwiązanie dokładne zależy. Można to zrobić przy pomocy wartości
%własnych(Było coś takiego na algebrze) macierzy A. Co to macierz A?
%Wspominałem w pliku z funkcją odefun że nasz układ możemy zapisać w
%postaci du/dt=A*u+b. Pamiętamy metody rozwiązywania układów równań
%liniowych za pomocą macierzy? Było na algebrze i pierwszej labce z wnum.
%
%   Jeżeli mamy układ:
%     y1'=2y1+3y2 + 5
%     y2'=4y1+6y2 + 8
%   To nasza macierz A=[2,3;4,6], a b=[5;8].
%
%Trzeba więc doprowadzić nasz układ do takiej postaci, żeby widzieć wyraźnie
%co stoi w obu równaniach przy u1 i u2, a po lewej stronie mieć tylko du1/dt
%lub du2/dt. Całe szczęście, że proszą nas tylko o zajęcie się funkcją u2.

R1=1e3;
R2=R1;
C1=1e-8;
C2=2e-6;
E=2;
%Nasza macierz A
A=[-2/(R1*C1), 1/(R1*C1); 1/(R2*C2), -1/(R2*C2)];

wart_wl=eig(A); %Do liczenia takich wartoœci własnych mamy ładnš funkcję eig()

lam1=wart_wl(1);%Mamy nasze lambda1 i lambda2
lam2=wart_wl(2);

c1=(lam2/(lam1-lam2)) * (E - u0(1)* (1+lam2*R2*C2)/(lam2*R2*C2)+ (u0(2)*1)/(lam2*R2*C2));
c2= -E-c1+u0(2);

%No i mamy nasze rowzišzanie dokładne. Potrzebujemy dwóch wersji, żeby
%póŸniej ładnie móc policzyć błędy. Tzn żeby wymiary macierzy się zgadzały.
u2d_45=c1*exp(lam1*t_45) + c2*exp(lam2*t_45) + E;
u2d_15s=c1*exp(lam1*t_15s) + c2*exp(lam2*t_15s) + E;

u2_45=u_45(:,2); % : znaczy "wszystko". Czyli bierzemy całš drugš kolumnę
u2_15s=u_15s(:,2);

figure
plot(t_45,u2d_45,'-g');
hold on;
grid on;
plot(t_45,u2_45,'--r');
legend('dokładne','wyliczone');
title('ode45');

figure
plot(t_15s,u2d_15s,'-g');
hold on;
grid on;
plot(t_15s,u2_15s,'--r');
legend('dokładne','wyliczone');
title('ode15s');

%===========b)=======================


%A jak policzyć błędy? Tak jak w zadaniu 1
max_err1_45=max(abs(u2d_45-u2_45))
max_err1_15s=max(abs(u2d_15s-u2_15s))
%No i widzimy, że ode_15s generuje troszkę większy błšd

%Tutaj mamy zobaczyć, że ode jest mšdre i samo na bieżšco dobiera długoœć
%kroku tak jak mu wygodnie.
diff1=diff(t_45);%diff() liczy nam różnicę po kolei między 1 a 2, 2 a 3, 3 a 4 elementem
diff2=diff(t_15s);

figure
semilogy(t_45(1:length(diff1)),diff1)
title('ode45');
xlabel('czas');
ylabel('długoœć kroku');

figure
semilogy(t_15s(1:length(diff2)),diff2)
title('ode15s');
xlabel('czas');
ylabel('długoœć kroku');
%===========c)=======================



for N=1:6
options=odeset('RelTol',1/(10^N));
%Funkcjom ode możemy zadawać dodatkowe parametry w taki oto sposób

[t_45,u_45]=ode45('odefun',tspan,u0,options);
[t_15s,u_15s]=ode15s('odefun',tspan,u0,options);

%Ode mogš zmieniac iloœć kroków czasowych, więc rozwišzania dokładne trzeba
%policzyć na nowo dla danej dziedziny czasu, w celu póŸniejszego wyznaczenia
%błędów.

u2d_45=c1*exp(lam1*t_45) + c2*exp(lam2*t_45) + E;
u2d_15s=c1*exp(lam1*t_15s) + c2*exp(lam2*t_15s) + E;

max_err_45(N)=max(abs(u_45(:,2)-u2d_45));
max_err_15s(N)=max(abs(u_15s(:,2)-u2d_15s));

t_len_45(N)=length(t_45);
t_len_15s(N)=length(t_15s);

%Nie wiem jaki zamysł mieli autorzy, ale gołym okiem różnicy nie widzę.
figure
plot(t_45,u_45(:,2))
str=sprintf("Ode45 dla RelTol= %1.0d",1/(10^N));
title(str);
figure
plot(t_15s,u_15s(:,2))
str=sprintf("Ode15s dla RelTol= %1.0d",1/(10^N));
title(str);
end

%Max błšd w zależnoœci od tolerancji
figure
semilogx([1e-1,1e-2,1e-3,1e-4,1e-5,1e-6],max_err_45)
title('Max błšd od RelTol dla Ode45')
xlabel('RelTol');
ylabel('Max błšd')

figure
semilogx([1e-1,1e-2,1e-3,1e-4,1e-5,1e-6],max_err_15s)
title('Max błšd od RelTol dla Ode15s')
xlabel('RelTol');
ylabel('Max błšd')

%Liczba kroków w zależnoœci od tolerancji
figure
semilogx([1e-1,1e-2,1e-3,1e-4,1e-5,1e-6],t_len_45)
title('Iloœć kroków od RelTol dla Ode45')
xlabel('RelTol');
ylabel('Iloœć kroków')

figure
semilogx([1e-1,1e-2,1e-3,1e-4,1e-5,1e-6],t_len_15s)
title('Iloœć kroków od RelTol dla Ode15s')
xlabel('RelTol');
ylabel('Iloœć kroków')

%=========d)============

%Metodami Eulera to bawiliœmy się już w zadaniu 1, ale teraz wyższa szkoła
%jazdy, bo musimy za ich pomocš rozwišzać układ równań. Pamiętamy, że nasz
%układ możemy zapisć w postaci du/dt=A*u+b gdzie A i b wyznaczyliœmy już wyżej:

A=[-2/(R1*C1), 1/(R1*C1); 1/(R2*C2), -1/(R2*C2)];
b=[E/(R1*C1);0];


%Euler otwarty (yn+1=yn+h*fn), gdzie fn=du/dt=A*u+b, a to już ładnie mamy

h=1e-6;
t=0:h:30e-3;
u0=[0;1];
u_eul_o=u0;
%Jedyna trudnoœć względem pojedyńczego równania to zapis, bo lecimy po
%kolumnach macierzy, a nie elementach tablicy.
for i=1:length(t)-1
    u_eul_o(:,i+1)=u_eul_o(:,i)+h*(A*u_eul_o(:,i)+b);
end
figure
plot(t,u_eul_o(2,:));
title('u2 - Euler Otwarty')

%Euler zamknięty (yn+1=yn+h*fn+1)

h=1e-6;
t=0:h:30e-3;
u0=[0;1];
u_eul_z=u0;

%Tutaj analogicznie jak w eulerze zamkniętym, trzeba wyznaczyć
%yn+1 z równania metody po podstawieniu, że fn+1=A*u(n+1)+b.

%Tu poprzednio był niedziałajšcy euler zamknięty, ale został naprawiony

%Kolega podrzucił działajšce rozwišzanie równania macierzowego, funkcję
%eye(), która robi nam macierz jednostkowš oraz pomnożenie przez macierz
%odwrotnš zamiast dzielenia. Sam bym na to nie wpadł

for i=1:length(t)-1
%Tak było
%   u_eul_z(:,i+1)=(u_eul_z(:,i)+h*b)\(1-A*h);
%A to już działa:
    u_eul_z(:,i+1)=inv((eye(2)-A.*h))*(u_eul_z(:,i)+h.*b);
end
figure
plot(t,u_eul_z(2,:));
title('u2 - Euler zamknięty')

%No to teraz trzebaby oba Eulery powrzucać w pętlę, pozliczać max_err w iteracjach,
%tak jak wczeœniej, tj max_err=max(abs(u_dok-u)) i wykreœlić loglog(max_err,[1e-6,2e-6...10e-6])
%Ale to już każdy raczej zrobi.

%=========e)========
%Podmienić te  R1,R2,C1,C2,E w tym pliku oraz odefun.m to nie problem ;)

%================================Ciekawostka==================================

%A jak ze stabilnoœciš eulera otwartego? W tym przypadku h*lambda musi należeć do [0,2],
%wpisujšc komendę 2/eig(A) dostajemy granicznš długoœć kroku dla której
%euler otwarty jest stabilny czyli ~0.98e-5. Wstawiajšc długoœć kroku np 1e-5
%zobaczymy piękny efekt niestabilnoœci metody eulera :)