%priklad 2.1, globalne Lipschitz
clear all;close all;clc;

fun = @(u,V) [V(1)*cos(u)];
exa = @(x) (exp(sin(x)));
T = [0 1];
Y0 = [1];

steps = zeros(8,1);
approx = zeros(8,1);
errs = zeros(8,1);
p = zeros(8,1);
estimate = zeros(8,1);
times = zeros(8,1);
M = 1.2748;
L = 1;
%chyby
for j = 1:8
    steps(j) = 1/(2^(j));
    [t Y] = euler(fun,Y0,T,steps(j));    
    f1 = feval(exa,t(end));    
    approx(j) = Y(end);
    times(j) = t(end);
    errs(j) = abs(Y(end) - f1);
    estimate(j) = ((M*steps(j))/(2*L))*(exp(L*1) - 1);  %vzorec (7)  
end

%numericky rad konvergence
p(1) = 0;
for j = 2:8
    p(j) = log(errs(j-1)/errs(j))/log(steps(j-1)/steps(j));
end

%tabulka
for k = 1:8
  fprintf('hi = 1/%d\n y = %5.4f, t = %5.4f, e = %5.4f, \n e_hat = %5.4f, p = %5.4f\n',...
  2^k,approx(k),times(k),errs(k),estimate(k),p(k));
  fprintf('--\n')
  end



%graf pro h = 1/32
s = linspace(0,1,500);

[t_n Y_n] = euler(fun,Y0,T,1/32);
figure
plot(t_n,Y_n,'--b','linewidth',2);hold on;xlabel('t');ylabel('y');
plot(s,exa(s),'r','linewidth',2);legend('numericke','presne','Location','northwest');
title('h = 1/32');text(0.5, 1.5,['y(t) = exp(sin t)'],'color','red');
hold off



exact = feval(exa,t(end));exact = feval(exa,t(end));
fprintf('presna hodnota v t = 1 je %f\n',exact)
%%
%================================================================================
%priklad 2.2
close all;clear all;clc

fun = @(t,Y) [Y(1)/t - (Y(1)/t)^2];
exa = @(x) (x./(1+log(x)));
T = [1 2];
Y0 = [1];

M = 1;
L = 3;

steps = zeros(8,1);
approx = zeros(8,1);
errs = zeros(8,1);
p = zeros(8,1);
estimate = zeros(8,1);
times = zeros(8,1);
%chyby
for j = 1:8
    steps(j) = 1/(2^(j));
    [t Y] = euler(fun,Y0,T,steps(j));
    f1 = feval(exa,t(end));    
    approx(j) = Y(end);
    times(j) = t(end);
    errs(j) = abs(Y(end) - f1);
    estimate(j) = ((M*steps(j))/(2*L))*(exp(L*1) - 1);
        
end

%numericky rad konvergence
p(1) = 0;
for j = 2:8
    p(j) = log(errs(j-1)/errs(j))/log(steps(j-1)/steps(j));
end

%tabulka
for k = 1:8
  fprintf('hi = 1/%d\n y = %5.4f, t = %5.4f, e = %5.4f, \n e_hat = %5.4f, p = %5.4f\n',...
  2^k,approx(k),times(k),errs(k),estimate(k),p(k));
  fprintf('--\n')
  end



%graf pro h = 1/32
s = linspace(1,2,500);

[t_n Y_n] = euler(fun,Y0,T,1/32);
figure
plot(t_n,Y_n,'--b','linewidth',2);hold on;xlabel('t');ylabel('y');
plot(s,exa(s),'r','linewidth',2);legend('numericke','presne','Location','northwest');
title('h = 1/32');text(1.4, 1.1,['y(t) = t/(1+ln t)'],'color','red');
hold off


exact = feval(exa,2);
fprintf('presna hodnota v t = 2 je %f\n',exact)

%%
%===========================================================================
%priklad 2.3 - v spocten v Pythonu. float 16 v moji verzi Octave neni. 

%priklad_2.3.py
%Dekuji za pochopeni. 

%%
%===========================================================================
%priklad 3.1 

close all;clear all;clc
%implicitni motivace
F = @(t,Y)(-20*Y(1));
Y0 = [1];
T = [0 1];
x = linspace(0,1,1000);
ex = @(t) exp(-20*t);

% forward h = 1/7
[t1 Y1] = euler(F,Y0,T,1/7);

%forward h = 1/8
[t2 Y2] = euler(F,Y0,T,1/8);

%forward h = 1/9
[t3 Y3] = euler(F,Y0,T,1/9);
figure;
plot(x,ex(x),'linewidth',2);hold on;
plot(t1,Y1,'linewidth',2);hold on;
plot(t2,Y2,'linewidth',2);hold on;
plot(t3,Y3,'linewidth',2);
legend('presne','h = 1/7','h = 1/8','h = 1/9','location','northwest');
axis([0 1 -10 10])



%backward h = 1/7
[b1 B1] = implicit_euler(F,Y0,T,1/7);

%backward h = 1/8
[b2 B2] = implicit_euler(F,Y0,T,1/8);

%backward h = 1/9
[b3 B3] = implicit_euler(F,Y0,T,1/9);

figure;plot(x,ex(x),'linewidth',2);hold on;
plot(b1,B1,'linewidth',2);hold on;
plot(b2,B2,'linewidth',2);hold on;
plot(b3,B3,'linewidth',2);
legend('presne','h = 1/7','h = 1/8','h = 1/9');hold off

%explicit stability h < 1/10
[e1 E1] = euler(F,Y0,T,1/11);
[e2 E2] = euler(F,Y0,T,1/12);
[e3 E3] = euler(F,Y0,T,1/20);
figure;plot(x,ex(x),'linewidth',2);hold on;
plot(e1,E1,'linewidth',2);hold on;
plot(e2,E2,'linewidth',2);hold on;
plot(e3,E3,'linewidth',2);
legend('presne','h = 1/11','h = 1/12','h = 1/20');hold off;

%compare when explicit stable
for k = 1:4
  h = 1/(20 + 10*k);
  [t Y]=implicit_euler_newton(F,Y0,T,h);
  [s S] = euler(F,Y0,T,h);
  
  maxErrIm(k) = max(abs(ex(t) - Y));
  maxErrEx(k) = max(abs(ex(s) - S)); 
  
  end
  %tabulka chyb 
  
  [maxErrEx' maxErrIm']

%%
%==================================================================
%priklad 4.5, srovnani dvou odhadu celkove chyby

clear all;close all;clc

F = @(t,Y) [Y(1)/t - (Y(1)/t)^2];
exa = @(x) (x./(1+log(x)));
T = [1 2];
Y0 = [1];

e = ones(8,1);
e_hat = ones(8,1);
h = ones(8,1);

M = 1;
L = 3;

%e_hat_old je odhad 2.11
%e_hat_new je odhad 4.19

for i = 1:8
    h(i) = 1/(2^i);
    [t_n1 Y_n1] = euler(F,Y0,T,h(i));
    [t_n2 Y_n2] = euler(F,Y0,T,h(i)/2);
    e(i) = abs(Y_n1(end) - exa(2));
    e_hat_new(i) = 2*abs(Y_n1(end) - Y_n2(end));  
    e_hat_old(i) = ((M*h(i))/(2*L))*(exp(L*1) - 1);
    fprintf('h_i = 1/%d, e_hat_old = %5.4f, e_hat_new = %5.4f\n',2^i,e_hat_old(i),e_hat_new(i));
end

[t Y] = adapt(F,Y0,T);
figure;plot(t,Y,'x','linewidth',2);hold on;
plot(t,exa(t),'linewidth',2);
legend('adapt','presne','location','southeast');title('h0 = 0.02, tol = 0.001')
%%
%=====================================================================
%priklad 4.6

clear all;close all;clc
%efekt promenlivosti kroku
F = @(x, X) [sin(x) - X(1)];
T = [0 10];
Y0 = 1;
[tn Yn] = adapt(F,Y0,T,0.2,0.01);
figure;plot(tn,Yn,'x','linewidth',2);

%%
%========================================================================
%priklad 5.7, system rovnic
clear all;close all;clc

F = @(t,X) [2*X(2) - 4*t;...
            -X(1) + X(3)-exp(t) + 2;...
            X(1) - 2*X(2) + X(3) + 4*t];
T = [0 1];
Y0 = [-1 0 2];

%exact
W = @(t) [-cos(2*t); sin(2*t) + 2*t; cos(2*t) + exp(t)];


%tabulka w/backward
N = 8;
e_2 = zeros(N,1);
e_inf = zeros(N,1);
h = zeros(N,1);
approx = zeros(N,3);

    
    for i = 1:N
        
        h(i) = 1/(2^i);
        [t_nt Y_nt] = implicit_euler_newton(F,Y0,T,h(i));
        e_2(i) = (sum((Y_nt(end,:) - W(1)').^2))^0.5;
        e_inf(i) = max(abs(Y_nt(end,:) - W(1)'));  
        approx(i,:) = Y_nt(end,:);
        disp(i);
        disp(Y_nt(end,:));
        
    end
%numerical rate
p_2(1) = 0;
p_inf(1) = 0;
for j = 2:N
    p_2(j) = log(e_2(j-1)/e_2(j))/log(h(j-1)/h(j));
    p_inf(j) = log(e_inf(j-1)/e_inf(j))/log(h(j-1)/h(j));
end


%tabulka 5.6
for k = 1:N
  fprintf('hi = 1/%d\n e_2 = %5.4f, p_2 = %5.4f, \n e_inf = %5.4f, p_inf = %5.4f\n',...
  2^k,e_2(k),p_2(k),e_inf(k),p_inf(k));
  fprintf('--\n')
  end
  

%graf pro h = 1/64
 t = 0:0.002:1;
 %the 3 components of the exact solution
 W1 = @(t) [-cos(2*t)];
 W2 = @(t) [sin(2*t) + 2*t];
 W3 = @(t) [cos(2*t) + exp(t)];
 
  %initial conditions
 F = @(t,X) [2*X(2) - 4*t;...
     -X(1) + X(3)-exp(t) + 2;...
     X(1) - 2*X(2) + X(3) + 4*t];
 T = [0 1];
 Y0 = [-1 0 2];
 %backward
 [t_n2 Y_n2] = implicit_euler_newton(F,Y0,T,1/64);
 
 figure
 %jeden obrazek
 
 plot(t,W1(t),'b-','linewidth',2);hold on
 plot(t_n2,Y_n2(:,1),'--','color',[0.25882, 0.92157, 0.95686],'linewidth',2);axis([0 1.02 -1 3]);
 text(0.89,0.52,['1. slozka'],'color','blue');
 hold on;

 plot(t,W2(t),'-','color',[0.77255, 0.95686, 0.25882],'linewidth',2);hold on
 plot(t_n2,Y_n2(:,2),'--','color',[0.27059,0.95686,0.25882],'linewidth',2);
 text(0.05,0.06,['2. slozka'],'color','g');
 hold on
 
 plot(t,W3(t),'r-','linewidth',2);hold on
 plot(t_n2,Y_n2(:,3),'--','color',[0.95686,0.44314,0.25882],'linewidth',2);hold on;
 text(0.01,2.2,['3. slozka'],'color','red');
 
 %chyba
err(1) = abs(Y_n2(end,1) - W1(1));
err(2) = abs(Y_n2(end,2) - W2(1));
err(3) = abs(Y_n2(end,3) - W3(1));

disp(err);
%===================================================================
%priklad 6.8
clear all;close all;clc


F = @(t,X) [32*X(1) + 66*X(2) + (2/3)*t + (2/3);...
-66*X(1) - 133*X(2) - (1/3)*t - (1/3)];

Y0 = [(1/3) (1/3)];
T = [0 0.5];
h = 0.025;

[t Y] = implicit_euler_newton(F,Y0,T,h);
[h H] = heun(F,Y0,T,h);

%exact solution
G1 = @(t,X) ((2/3)*t + (2/3)*exp(-t)-(1/3)*exp(-100*t));
G2 = @(t,X) (-(1/3)*t - (1/3)*exp(-t) + (2/3)*exp(-100*t));
z = 0:0.002:0.5;

%implicit graph
figure;plot(t,Y(:,1),'--','color',[0.2, 0.2, 1],'linewidth',2);hold on;
plot(z,G1(z),'-','color',[0.4, 0.4, 1],'linewidth',2);hold on;
text(0.2, 0.6,['1. slozka'],'color',[0.2, 0.2, 1])
plot(t,Y(:,2),'--','color',[0, 0.6, 0],'linewidth',2);hold on;
plot(z,G2(z),'-','color',[0.4, 0.8, 0],'linewidth',2);
text(0.1, -0.23,['2. slozka'],'color',[0, 0.4, 0])
title('Implicitni h = 0.025')


%heun graph
figure;plot(h,H(:,1),'--','color',[0.2, 0.2, 1],'linewidth',2);hold on;
plot(z,G1(z),'-','color',[0.4, 0.4, 1],'linewidth',2);hold on;
title('Heun, h_1 = 0.025');text(0.1, -436,['1. slozka'],'color',[0.2, 0.2, 1])

%smaller h
h = 0.00625;
[t2 Y2] = implicit_euler_newton(F,Y0,T,h);
[h2 H2] = heun(F,Y0,T,h);

figure;plot(t2,Y2(:,1),'--','color',[0.2, 0.2, 1],'linewidth',2);hold on;
plot(z,G1(z),'-','color',[0.4, 0.4, 1]);hold on;
plot(t2,Y2(:,2),'--','color',[0, 0.6, 0],'linewidth',2);hold on;
plot(z,G2(z),'-','color',[0.4, 0.8, 0]);hold on;
text(0.05, 0.58,['1. slozka'],'color',[0.2, 0.2, 1]);
text(0.08, -0.22,['2. slozka'],'color',[0, 0.4, 0]);
plot(h2,H2(:,1),'-.','linewidth',2,'color',[0.95, 0.259, 0.27]);hold on;
plot(h2,H2(:,2),'-.','linewidth',2,'color',[0.95, 0.725, 0.258])
title('h = 0.00625')

%differences in graphs not well visible 
%do errors

n = length(t2);
W = @(t) [(2/3)*t + (2/3)*exp(-t)-(1/3)*exp(-100*t);...
-(1/3)*t - (1/3)*exp(-t) + (2/3)*exp(-100*t)];

presne = reshape(W(t2),length(t2),2);
errim = norm((Y2 - presne),inf)
errheun = norm((H2 - presne),inf)

%is implicit euler better?
errim < errheun
%no

%%
%============================================================================
%euler-cromer, kyvadlo
clear all;close all;clc

g = 9.81;
l = 1;
F = @(t,X) [X(2);-(g/l)*X(1)];
Y0 = [0.2 0];
h = 0.04;
tf = 10;

%exact solution
foo = @(t) [0.2*cos(sqrt(g/l).*t)];

%=============
%forward euler
[t_e Y_e] = euler(F,Y0,[0 tf],h);

%compare exolution w/numeric  
figure;plot(t_e,foo(t_e),'linewidth',2);hold on;plot(t_e,Y_e(:,1),'linewidth',2);
legend('presne','euler, h=0.04');

%backward euler
[b B] = implicit_euler_newton(F,Y0,[0 tf],h);

%compare exolution w/numeric
figure;plot(b,foo(b),'linewidth',2);hold on;plot(b,B(:,1),'linewidth',2);
legend('presne','implicit, h=0.04');

%euler cromer
i = 1;
g = 9.81;
l = 1;
Y(i,1) = 0.2;Y(i,2) = 0;

t(i) = 0;
tf = 10;
h = 0.04;
%exact solution
foo = @(t) [0.2*cos(sqrt(g/l).*t)];
%energy
E = @(t,Y) [0.5*Y(2)^2 + 0.5*(g/l)*Y(1)^2];
%initial energy
EK(i) = 0.5*Y(i,2)^2;
EP(i) = 0.5*(g/l)*Y(i,1)^2;
%t(end) might not be tf

while t(i) < tf
  Y(i+1,2) = Y(i,2) -  (g/l) * Y(i,1) * h;
  Y(i+1,1) = Y(i,1) + Y(i+1,2) * h;
  t(i+1) = t(i) + h;
  i = i + 1; 
  EK(i) = 0.5*Y(i,2)^2;
  EP(i) = 0.5*(g/l)*Y(i,1)^2;  
  end
t = t';
[t Y];


figure;plot(t,foo(t),'linewidth',2);hold on;plot(t,Y(:,1),'--','linewidth',2);
axis([0 tf]);legend('presne','Euler-Cromer')

%energy cromer
for q = 1:length(t)
  energy_c(q) = E(t(q),Y(q,:));
  end

figure;plot(t,EK,'linewidth',2);hold on;
plot(t,EP,'linewidth',2);hold on;
plot(t,energy_c,'linewidth',2);xlabel('cas');ylabel('energie');
axis([0 tf]);legend('kinetic','potential','total','location','northeast');
fprintf('pocatecni celkova energie %f\n',energy_c(1));
fprintf('koncova celkova energie %f\n',energy_c(end));



%==============================================================================
%priklad 7.9, volny pad s odporem vzduchu
clear all;close all;clc

%fyzikalni konstanty
g = 9.81; %m/sec^2
b = 0.2;
cor = 0.5;

%bariera
G = @(u,W) [-W(1)];


%pocatecni uloha
F = @(t,X) [X(2); -9.81 - b*X(2)];
Y0 = [1 0];

%reseni pomoci Euler + bisekce
h = 0.5;
tf = 1;
tol = 0.01;
[t Y] = bar_puleni(F,G,Y0,tf,h,tol,cor);

figure;plot(t,Y(:,1),'x','linewidth',2);hold on;plot(t,Y(:,1));
set (gca,'xtick',t);
set (gca,'xticklabel',{});
axis([0 tf 0 1]); title('h = 0.5, t v [0,1], tol = 0.01');

%%
%==========================================================================
%priklad 7.10
clear all;close all;clc

%fyzikalni konstanty
cor = 1;

%bariera
G = @(u,W) [(W(1)-(1/8))*(W(1)+(1/8))];

%pocatecni uloha
F = @(t,X) [X(2); 2];
Y0 = [0 -0.8568];
tf = 1;

%reseni pomoci Euler + bisekce
h = 0.04;
tol = 0.001;
[t Y] = bar_puleni(F,G,Y0,tf,h,tol,cor); 

figure;plot(t,Y(:,1),'x','linewidth',2);hold on;plot(t,Y(:,1));hold on;
set(gca,'ytick',[-0.15, -0.125, -0.1, -0.05, 0, 0.05, 0.1, 0.125, 0.15]);
set(gca,'yticklabel',{'0.15', '-1/8', '-0.1', '-0.05', '0', '0.05', '0.1','1/8','0.15'});
title('h = 0.04, tol = 0.001, cor = 1')



%clanek example 2.1 s v = -1.76579
Y0 = [0 -1.76579];
[t2 Y2] = bar_puleni(F,G,Y0,tf,h,tol,cor);

hold on;
plot(t2,Y2(:,1),'x','linewidth',2);hold on;plot(t2,Y2(:,1));
hold on;plot([0 1],[0.125 0.125],'g');
hold on;plot([0 1],[-0.125 -0.125],'g');
legend('v0 = -0.8568','','v0 = -1.76579','','location','northeastoutside')

%%
%================================================================

%priklad 7.12

clear all;close all;clc

%fyzikalni konstanty
cor = 1;

%bariera
G = @(u,W) [(W(1)-(3/8))*(W(1)+(3/8))];

%pocatecni uloha
F = @(t,X) [X(2);6*t];
Y0 = [0 -1];
tf = 1;

%reseni pomoci Euler + bisekce
h = 0.04;
tol = 0.001;
cor = 1;

[t Y] = bar_puleni(F,G,Y0,tf,h,tol,cor);


%druha volba pocatecni podminky
Y0 = [0 -1.218];

[t2 Y2] = bar_puleni(F,G,Y0,tf,h,tol,cor);

%graficke srovnani
figure;plot(t,Y(:,1),'x','linewidth',2);hold on;plot(t,Y(:,1));hold on;
set (gca,'ytick',[-0.4, -0.375, -0.2, 0, 0.2, 0.375, 0.4])
set (gca,'yticklabel',{'-0.4','-3/8','-0.2','0','0.2','3/8','0.4'});
title('h = 0.04, tol = 0.001, cor = 1')
hold on;plot(t2,Y2(:,1),'x','linewidth',2);hold on;
plot(t2,Y2(:,1));hold on;
legend('v0 = -1','','v0 = -1.218','','location','northeastoutside')
plot([0 1],[-0.375 -0.375],'g');
hold on;plot([0 1],[0.375 0.375],'g');

%%
%==============================================================
