%                   Numerik LAK: Einschrittverfahren für ODE
%


% exakte Lösung

x = 1:5;
y_ex= x.^2+2*x+2 - 2*(x+1).*log(x+1);

% numerische Lösungen

bsp = @bsp_2;
h=.2;opt=1;  y = init_value(bsp,0,2,5,h,opt);

y1 = y(6:5:end)       % num. Lösung an x = 1,2,3,4,5
e1 = y_ex' - y1

% Halbieren der Schrittweite

h=.1;opt=1;  y = init_value(bsp,0,2,5,h,opt);

y2 = y(11:10:end);
e2 = y_ex' - y2

h=.05;opt=1;  y = init_value(bsp,0,2,5,h,opt);

y3 = y(21:20:end);
e3 = y_ex' - y3

% Rate

[e3./e2 e2./ e1]

% Heunverfahren

h=.2;opt=2;  y = init_value(bsp,0,2,5,h,opt);

y1 = y(6:5:end);       % num. Lösung an x = 1,2,3,4,5
e1 = y_ex' - y1

% Halbieren der Schrittweite

h=.1;opt=2;  y = init_value(bsp,0,2,5,h,opt);

y2 = y(11:10:end);
e2 = y_ex' - y2

h=.05;opt=2;  y = init_value(bsp,0,2,5,h,opt);

y3 = y(21:20:end);
e3 = y_ex' - y3

% Rate

[e3./e2 e2./e1]

% Runge-Kutta

h=.2;opt=3;  y = init_value(bsp,0,2,5,h,opt);

y1 = y(6:5:end);       % num. Lösung an x = 1,2,3,4,5
e1 = y_ex' - y1

% Halbieren der Schrittweite

h=.1;opt=3;  y = init_value(bsp,0,2,5,h,opt);

y2 = y(11:10:end);
e2 = y_ex' - y2

h=.05;opt=3;  y = init_value(bsp,0,2,5,h,opt);

y3 = y(21:20:end);
e3 = y_ex' - y3

% Rate

[e3./e2 e2./e1]


% Beispiel 3

clear
y_ex = exp(-30*1.5)/3;

bsp = @bsp_3;
h=.1;
opt=1;  y1 = init_value(bsp,0,1/3,1.5,h,opt);y_ex-y1(end)
opt=2;  y2 = init_value(bsp,0,1/3,1.5,h,opt);y_ex-y2(end)
opt=3;  y3 = init_value(bsp,0,1/3,1.5,h,opt);y_ex-y3(end)


% Beispiel 4

clear

bsp = @bsp_4;
x0 = [80 30];

% Euler

h=1;opt=1;  y1 = init_value(bsp,0,x0,12,h,opt); plot(y1(:,1),y1(:,2),'.r'), hold on
h=.5;opt=1;  y1 = init_value(bsp,0,x0,12,h,opt); plot(y1(:,1),y1(:,2),'.b')
h=.25;opt=1;  y1 = init_value(bsp,0,x0,12,h,opt); plot(y1(:,1),y1(:,2),'.g')

% Runge Kutta

h=1;opt=3;  y1 = init_value(bsp,0,x0,12,h,opt); plot(y1(:,1),y1(:,2),'.r'), hold on
h=.5;opt=3;  y1 = init_value(bsp,0,x0,12,h,opt); plot(y1(:,1),y1(:,2),'.b')
h=.25;opt=3;  y1 = init_value(bsp,0,x0,12,h,opt); plot(y1(:,1),y1(:,2),'.g')

%
% Ende von lecture1.m
%

function y = init_value(dgl,x0,y0,xf,h,opt);
% y = init_value(dgl,x0,y0,xf,h,opt);
% opt = 1   Euler's method
% opt = 2 Heun's method
% opt = 3 Runge-Kutta, 4th order
%
% y0 initial value (rowvector)
% x0 initial time
% xf final time
% h stepsize
% dgl function handle pointing to specific example
% output:
% y solution at x = i*h, i=0,..., n

% number of steps and initialize
n = ceil((xf-x0)/h);
[mm,m] = size(y0);
y= zeros(n+1,m);
y(1,:) = y0;
x  =x0;

if opt ==  1
    for i=1:n
        y(i+1,:) = y(i,:) + h*feval(dgl,x,y(i,:));
        x= x+h;
    end
elseif opt == 2
     for i=1:n
         z1 = feval(dgl,x,y(i,:));
         x= x+h;
         z2 = feval(dgl, x, y(i,:) + h*z1);
         y(i+1,:) = y(i) + .5*h* (z1 + z2);
     end
elseif opt ==3
    h2 = h/2;
    for i=1:n
        z1 = feval(dgl,x,y(i,:));
        x= x+h2;
        z2 = feval(dgl, x, y(i,:) + h2*z1);
        z3 = feval(dgl, x, y(i,:) + h2*z2);
        x= x+h2;
        z4 = feval(dgl, x, y(i,:) + h*z3);
         y(i+1,:) = y(i,:) + h* (z1 + 2*(z2+z3)+z4)/6;
    end
end
         
%
% Ende von init_value.m
%


function f = bsp_2(x,y);
%
f = (y+x^2-2)/(x+1);

function f = bsp_3(x,y)
%
f = -30*y;

function f = bsp_4(x,y)
%
f = [.25*y(1)-.01*y(1)*y(2), -y(2)+.01*y(1)*y(2)];