% Exercise sheet 3
% SDE dX_t = A X_t dt + B dB_t 
% A = A2 here (replace in line 25 with A1, A2, A4 to check these)

clear
close all

X(:,1) = [1; 2]; % initial value
Y(:,1) = [1; 2];
t(1) = 0;
dt = 0.01; % timestep
N = 1000;

A1 = [1 0; 0 -1];
A2 = [-1 1; 1 -1];
A3 = [0 1; -1 0];
A4 = [-1 0.5; 0.5 -1];

B = [1 0; 0 1];

for n = 2:N
    r1 = normrnd(0,dt); %iid gaussian rv's
    r2 = normrnd(0,dt);
    r = [r1; r2];
    X(:,n) = X(:,n-1) + A2*X(:,n-1)*dt + B*r; %euler method for SDE
    t(n) = t(n-1) + dt;
    Y(:,n) = expm(A2*t(n))*X(:,1); %deterministic solution
end

eig(A1)
eig(A2)
eig(A3)
eig(B)

plot(t,X(1,:))
hold on
plot(t,Y(1,:),'r');
xlabel('t')
ylabel('X1')
%title('X1')
hold off;
figure
plot(t,X(2,:))
hold on
plot(t,Y(2,:),'r');
xlabel('t')
ylabel('X2')
hold off;