% Řešení obyčejných diferenciálních rovnic

% symbolické řešení - dsolve

% y''(t) = a, y(0) = 1, y'(0) = 2
y = dsolve('D2y = a', 'y(0) = 1', 'Dy(0) = 2') % = (a*t^2)/2 + 2*t + 1
a = 1.2
t = 0:100;
plot(t, eval(y))

% y' = -x, x' = -y, x(0) = pi/2, y(0) = pi/2
[x, y] = dsolve('Dy = -x', 'Dx = -y', ...
                'x(0) = pi/2', 'y(0) = pi/2') % x = pi/(2*exp(t)), y = pi/(2*exp(t))

% y' = sin(yt), y(0) = 0.1
%% y = dsolve('Dy = sin(y*t)', 'y(0) = 0.1') % explicit solution could not be found

% numerické řešení - ode23, ode45, ...

% y' = sin(yt)
fce_inline = inline('sin(y*t)', 't', 'y') % inline deklarace funkce
[t, y] = ode45(fce_inline, [0, 10], 0.1)
plot(t, y, '*-')

%deklarace funkce v souboru
%function dy = fce_file(t, y)
%dy = sin(y*t)
[t, y] = ode45(@fce_file, [0, 10], 0.1)

% y' = -x, x' = -y, x(0) = pi/2, y(0) = pi/2

%function dz = fce2(t, z)
%dz = [-z(2); -z(1)];

[t, z] = ode45(@fce2, [0, 10], [pi/2, pi/2])
plot(t, z(:,1)) % vykreli reseni x(t)

% y''(t) = a, y(0) = 1, y'(0) = 2, t in <0, 100>
% rovnici 2. řádu nutno převést na soustavu 1. řádu
% y1 = y, y2 = y' ==> y1' = y2, y2' = y''

%function dy = fce3(t, y, a)
%dy = [y(2); a];

[t, z] = ode45(@(t, z) fce3(t, z, 1.2), [0, 100], [1, 2])
subplot(2, 1, 1)
plot(t, z(:,1))
title('y(t)')
subplot(2, 1, 2)
plot(t, z(:,2))
title('dy(t)')
clf

% y'' + 5y' + 6y = exp(-x), y(0) = 0, y'(0) = 0.1, x in <0, 10>
% y1 = y, y2 = y' ==> y1' = y2, y2' = exp(-x) - 5y2 - 6y1

%function dy = fce4(x, y)
%dy = [y(2);
%      exp(-y) - 5 * y(2) - 6 * y(1)];

[x, y] = ode45(@fce4, [0, 10], [0, 0.1])
plot(x, y)
legend('y', 'dy')

% 2y''' + y'' = 0, y(0) = 0, y'(0) = 0, y''(0) = 0.332, t in <0, 5>
% y1 = y, y2 = y', y3 = y'' ==> y1' = y2, y2' = y3, y3' = -1/2 y3

%function dy = fce5(t, y)
%dy = [y(2);
%      y(3);
%      -1/2 * y(3)];

[t, y] = ode45(@fce5, [0, 5], [0, 0, 0.332])
plot(t, y)
legend('y', 'dy', 'd2y')

% nastavení parametrů řešiče ODE - odeset
[t, z] = ode45(@(t, z) fce3(t, z, 1.2), [0, 100], [1, 2], ...
               options) % 53 časových kroků
options = odeset('AbsTol', 1e-1)
[t, z] = ode45(@(t, z) fce3(t, z, 1.2), [0, 100], [1, 2], ...
               options) % 41 časových kroků

% funkce ODE
%   ode45 - 4. řád
%   ode23 - 2. řád
%   ode113 - proměnný řád
%   ode15s - proměnný řád, stiff úlohy
%   ode23s - 2. řád, stiff úlohy
%   ...

% příklad:
% kmitání 1D lineárního harmonického oscilátoru
%   rovnice: m * x''(t) + b * x'(t) + k * x(t) = 0
%   počáteční podmínky: x(0) = x0, x'(0) = 0
%   časový interval: <0, 40>
%
% převod na soustavu rovnic 1. rádu:
% x1 = x, x2 = x'
%  ==> x1' = x2, x2' = -(b * x2 + k * x1) / m

%function dx = fce_osc(t, x, m, b, k)
%dx = [x(2);
%      -(b * x(2) + k * x(1)) / m];

[t, x] = ode45(@(t, x) fce_osc(t, x, 1, 0.4, 2), [0, 40], [1, 0])
plot(t, x)
legend('x', 'dx')