%% Rozwiązywanie równań różniczkowych
% opracował dr inż. Robert JAKUBOWSKI, Politechnika Rzeszowska, KSiSL
%% Rozwiązywanie równań różniczkowych z wykorzystaniem metod numerycznych
% W poprzednim wykładzie zaprezentowane zostało podejście do rozwiązywania równań różniczkowych z wykorzystaniem Symbolic toolbox. Tym razem zaprezentuję podejście do rozwiąywania równań różniczkowych z wykorzystaniem funkcji ode - ordinary differential equations.
% Jednym z narzedzi do rozwiązywania rónań i układów równań różniczkowych jest funkcja ode45. Zagadnienia związane z rozwiąywaniem równań różniczkowych w technice najczęściej stosowane są do zgadnień w dziedzinie czasu - związane z ruchem w czasie. Funkcję ode45 można użyć w następujący sposób:
% [t,y] = ode45(@odefun,tspan,y0)
% gdzie:
% t - wektor czasu
% y - wektor rozwiąznia funkcji dla określonego czasu t
% odefun - nazwa funkcji zawierającej równanie różniczkowe
% tspan - przedział czasu
% y0 - wartość y dla początkowej chwili czasu t0
% Przedsyawmy to na przykładzie. Poszukać rozwiązania równania różniczkowego opisanego równaniem:
% \frac{\textrm{dy}}{\textrm{dt}}=\sin \left(t\right)
% w przedziale czasu tspan=[0 , 4pi), gdy dla poczatkowej chwili czasu t0 y0=0.
% W celu rozwiązania równania przgotowuje się funkcje  zmiennej y i t. Taką funkcję przedsyawiono poniżej w funkcjech pomocniczych pod nazą odefun1.
% Zdefinujmy zatem przedział czasu tspan oraz wartość y0 dla poczatkowej chwili czasu t0

tspan=[0 4*pi]
y0=-1;

% Wywołajmy funkcję ode45 otrzymujemy rozwiązanie, czyli przebieg zmiennej y w czasie t:

[t,y] = ode45(@odefun1,tspan,y0);

% Narysujmy tą funkcję na wykresie

plot(t,y)
xlabel('t');
ylabel('y')
grid on

%
% Tym razem rozwiążmy zadanie, gdy równanie różniczkowe wygląda następująco
% \frac{\textrm{dy}}{\textrm{dt}}=y\left(1-\frac{y}{K}\right)
% gdzie k=100
% Przedstawmy rozwiązanie równania w przedziale czasu tspan=[0,10]s, gdy dla poczatkowej chwili czasu t0 y0=1.
% W celu rozwiązania równania przgotowuje się funkcje zmiennej y i t. Taką funkcję przedsyawiono poniżej w funkcjech pomocniczych pod nazą odefun2.
% Zdefinujmy zatem przedział czasu tspan oraz wartość y0 dla poczatkowej chwili czasu t0

tspan=[0 10]
y0=1;

% Wywołajmy funkcję ode45:

[t,y] = ode45(@odefun2,tspan,y0);

% Wynik przedstawmy na wykresie

plot(t,y)
xlabel('t')
ylabel('y')
grid on

% W rozwiązaniu otrzymaliśmy wektor czasu i odpowiadające wartości y. W ten sposób gdy chcemy znaleźć wartość y dla konkretnej chwili czasy t, możemy mieć problem. Dane są wartości y dla określonych wartości czasu t. Na przykład w tym przypadku dostaliśmy wektor czasu t o n elementach:

n=length(t)
disp(t')

% Gdy interesuje nas wynik dla prezentowanej chwili czasu t np. dla t(25), to możemy otrzymać wartość y w sposób następujący

t(25)
y(25)

% Przedstawmy otrzymany punkt na wykresie

hold on
plot(t(25),y(25),'pr')
hold off

% Problem pojawia się, gdy chcemy otrzymać rozwiązanie dla chwili czasu t, która nie jest zawarta w wektorze t. W takim przypadku, możemy wywołać funkcję ode45 nieco inaczej, przypisując jej tylko jedną wartość wyjściową:

rozw = ode45(@odefun2,tspan,y0)

% W tym przypadku zmienna rozw jest strukturą zawierającą rozwiązanie, które potem w sposób łatwy z użyciem funkcji deval można przeliczyć dla dowolnej chwili czasu t, dla której wcześniej wyznaczyliśmy rozwiązanie. Naprzykład teraz wyznaczmy rozwiązanie dla t=2 i przedstawimy rozwiązanie na wykresie

% Wyznaczanie wartości y od t=2 dla wcześniej wyznaczonego rozwiązania jako
% struktura rozw
y2=deval(rozw,2)
% Powtórzenie wykresu na podstawie wcześniej wyznaczonych punktów
plot(t,y)
hold on
xlabel('t')
ylabel('y')
plot(2,y2,'pr')
% Wyznaczenie wartości y od t=2.2222
y2_222=deval(yodt,2.222)
plot(2.222,y2_222,'pc')
% Wyznaczenie wartości y dla trzech wartośći chwili czasu t
y4_8=deval(yodt,[4,6,8])
plot([4,6,8],y4_8,'gp')
grid on
hold off

%
%% Zadanie do wykonania w ramach ćwiczeń
% Wyznaczyć wartość rozwiązania równania w przedziale czasu t=0:2pi
% \frac{\textrm{dx}}{\textrm{dt}}=2\sin \left(x\right)-\cos \left(t^2 \right)
% początkową wartość rozwiązania w chwili czasu t=0 przyjąć x0=1
% Przedstawić rozwiązanie graficznie linią niebieską ciągłą
% Następnie zmienić wartość poczatkową x0=0, dla chwili t0, rozwiązanie przedstawić na tym samym wykresie linią przerywaną w kolorze czerwony,
% Uwaga : Funkcjach na dole w Funkcje pomocnicze należy przygotować stosowną funkcję zawirającą równanie



%
%
%% Funkcje pomocnicze
% Funkcja, przygotowana do rozwiązania równania różniczkowego za pomocą funkcji ode45, zawierająca równanie
% \frac{\textrm{dy}}{\textrm{dt}}=\sin \left(t\right)
% Funkcja odefun1 (nazwa może być dowolna). Po nzawie funkcji zmienne niezależna  t i zmiennea zależna y, która w tym wypadku nie występuje w równaniu

function dydt = odefun1(t,y)

dydt = sin(t);
end

% Funkcja odefun2 przygotowana do rozwiązania równania różniczkowego za pomocą funkcji ode45, zawierająca równanie
% \frac{\textrm{dy}}{\textrm{dt}}=y\left(1-\frac{y}{K}\right)
% Funkcja odefun2. Po nzawie funkcji zmienne niezależna  t i zmiennea zależna y

function dydt = odefun2(t,y)
K = 100;
dydt = y*(1 - y/K);
end

% Funkcja, przygotowana do rozwiązania równania różniczkowego
% \frac{\textrm{dx}}{\textrm{dt}}=2\sin \left(x\right)-\cos \left(t^2 \right)

