
% - - - - - - - - -- - - - -- - - 
% vConstantin Greif
% Stochastic interpretation
% L_k used is falcone + lions
% Semi Lagragian Schemes
% - - - - - - - - - - -- - -- - 

clc
clear all
close all



%% Example 1:


f_fct = @(alp_t, t,x,y) 0;
%sin(x)*sin(y)*((1)*(2-t) -1) -2*(2-t)*cos(x)*cos(y)*sin(x+y)*sin(x+y) ;
%g_fct = @(x,y) 0;
sig_fct = @(alp_t,x,y) [sqrt(2)*sin(x+y), sqrt(2)*sin(x+y) ; 2, 0];
mu_fct = @(alp_t,x,y) [2, 0];

u_sol = @(t,x,y) (2-t)*sin(x)*sin(y);
%u_sol = @(t,x,y) exp(t)*exp(sqrt(0.5)*x)*exp(sqrt(0.5)*y); 

h = 1/8; % time step
h_t = h;
T = 1;
p = 2; % dimensions of sig

x_vec = -1:h:1;
y_vec = -1:h:1;
t_vec = 0:h_t:T;

n_t = length(t_vec);
n_x = length(x_vec);
n_y = length(y_vec);

Lamb = [1,-1];
curr = 0;
%mini = 0;
u_mat = zeros(n_x,n_y);

% A_M = zeros(n_t,2);
% A_M(:,1) = -1;
% A_M(:,2) = 0;
alpha = 0; % this valule will be used as first test control 
%[alpha,u_num] = stoch_SL_scheme(f_fct,g_fct,sig_fct,b_fct,A_M,x_vec,t_vec);

%% Example 2:


f_fct = @(lam1,lam2,t,x,y) (1-t)*sin(x)*sin(y)-2*lam1*lam2*(2-t)*cos(x)*cos(y);
%sin(x)*sin(y)*((1)*(2-t) -1) -2*(2-t)*cos(x)*cos(y)*sin(x+y)*sin(x+y) ;
%g_fct = @(x,y) 0;
sig_fct = @(lam1,lam2,x,y) [sqrt(2)*lam1, sqrt(2)*lam2 ; 0, 0];
mu_fct = @(alp_t,x,y) [0, 0];

u_sol = @(t,x,y) (2-t)*sin(x)*sin(y);
%u_sol = @(t,x,y) exp(t)*exp(sqrt(0.5)*x)*exp(sqrt(0.5)*y); 

h = 1/20; % time step
h_t = h;
T = 1;
p = 2; % dimensions of sig

x_vec = -1:h:1;
y_vec = -1:h:1;
t_vec = 0:h_t:T;

n_t = length(t_vec);
n_x = length(x_vec);
n_y = length(y_vec);

Lamb = -1:0.1:1;
curr = 0;
%mini = 0;
u_mat = zeros(n_x,n_y);

% A_M = zeros(n_t,2);
% A_M(:,1) = -1;
% A_M(:,2) = 0;
alpha = 0; % this valule will be used as first test control 
%[alpha,u_num] = stoch_SL_scheme(f_fct,g_fct,sig_fct,b_fct,A_M,x_vec,t_vec);

%% Example 3:


f_fct = @(alp1, alp2,t,x,y) 0;
%sin(x)*sin(y)*((1)*(2-t) -1) -2*(2-t)*cos(x)*cos(y)*sin(x+y)*sin(x+y) ;
%g_fct = @(x,y) 0;
sig_fct = @(alp1,alp2,x,y) [sqrt(sqrt(2)), 0 ; 0, sqrt(sqrt(2))];
mu_fct = @(alp_t,x,y) [0, 0];

u_sol = @(t,x,y) exp(t)*exp(sqrt(0.5)*x)*exp(sqrt(0.5)*y);
%u_sol = @(t,x,y) exp(t)*exp(sqrt(0.5)*x)*exp(sqrt(0.5)*y); 

h = 1/20; % time step
h_t = h;
T = 1;
p = 2; % dimensions of sig

x_vec = -1:h:1;
y_vec = -1:h:1;
t_vec = 0:h_t:T;

n_t = length(t_vec);
n_x = length(x_vec);
n_y = length(y_vec);

Lamb = [1];
curr = 0;
%mini = 0;
u_mat = zeros(n_x,n_y);

% A_M = zeros(n_t,2);
% A_M(:,1) = -1;
% A_M(:,2) = 0;
alpha = 0; % this valule will be used as first test control 
%[alpha,u_num] = stoch_SL_scheme(f_fct,g_fct,sig_fct,b_fct,A_M,x_vec,t_vec);


%%
u_til = zeros(n_t,n_x,n_y);
for i = 1:n_x
    for j = 1:n_y
        for k = 1:n_t
            if (i == 1)||(i == n_x) || (j == 1) || (j == n_y) || ( k == 1)
                u_til(k,i,j) = u_sol(t_vec(k),x_vec(i),y_vec(j));
            end
        end
    end
end

for k = 1:n_t-1 % time steps
    
    for i = 2:n_x-1
        
        for j = 2:n_y-1
        
            x = x_vec(i);
            y = y_vec(j);
            t = t_vec(k+1);
            for ii = 1:n_x
                for jj = 1:n_y
                  u_mat(ii,jj) = u_til(k,ii,jj);  
                end
            end
            kay = sqrt((p+1)*h_t);
            
            for la = 1:length(Lamb)
                
                lamb = Lamb(la);
                lamb2 = ((-1)^la)*sqrt(1-lamb^2);
                sig = sig_fct(lamb,lamb2,x,y);
                mu = mu_fct(lamb,x,y);
                
                error = bilin_interpol(u_mat,x_vec,y_vec,x,y) - u_til(k,i,j);
            
                L = ( bilin_interpol(u_mat,x_vec,y_vec,x + (kay^2)*mu(1),y + (kay^2)*mu(2) ) - bilin_interpol(u_mat,x_vec,y_vec,x,y)  )/(kay^2) + ...
                    ( bilin_interpol(u_mat,x_vec,y_vec,x + (kay)*sig(1,1),y + (kay)*sig(2,1) ) - 2*bilin_interpol(u_mat,x_vec,y_vec,x,y) + bilin_interpol(u_mat,x_vec,y_vec,x - (kay)*sig(1,1),y - (kay)*sig(2,1) ) )/(2*kay^2) + ...
                    ( bilin_interpol(u_mat,x_vec,y_vec,x + (kay)*sig(1,2),y + (kay)*sig(2,2) ) - 2*bilin_interpol(u_mat,x_vec,y_vec,x,y) + bilin_interpol(u_mat,x_vec,y_vec,x - (kay)*sig(1,2),y - (kay)*sig(2,2) ) )/(2*kay^2);
            
                curr = h_t*f_fct(lamb,lamb2,t,x,y) + h_t*L + u_til(k,i,j);
                
                if la == 1
                    mini = curr;
                end
                
                if curr < mini
                    mini = curr;
                end
                
            
            end
        
        u_til(k+1,i,j) = mini;
        
        end
    end
end




%% plot
u_plot = zeros(n_x,n_y);
for i = 1:n_x
    for j = 1:n_y
   u_plot(i,j) = u_til(n_t,i,j); 
    end
end
subplot(2,1,1)
%u_plot = u_num(n_t,:,:);
surf(x_vec,y_vec,u_plot');
xlabel('x');
ylabel('y');
xlim([-1 ,1]);
ylim([-1, 1]);
zlim([0, 12]);

%u_sol = @(x,y) (2-1)*sin(x)*sin(y);
u_sol = @(x,y) exp(1)*exp(sqrt(0.5)*x)*exp(sqrt(0.5)*y);   % sol
subplot(2,1,2);
ezsurf(u_sol,[-1, 1, -1, 1]);
