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

clc
clear all
close all

% NEED TO CHANGE INTEROPL FOR MEHR DIMENSIONAL

%% Example 1:


f_fct = @(alp_t, x,y) (1 -x^2);
g_fct = @(x,y) 0;
sig_fct = @(alp_t,x,y) alp_t;
mu_fct = @(alp_t) 0;

u_sol = @(t,x,y) t*(1-x^2); 

h = 0.12; % time step
T = 1;
p = 2; % dimensions

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

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

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 v = 1:y_x
    u_til(1,i,v) = g_fct(x_vec(i),y_vec(v)); 
end
end

for j = 1:n_t-1 % time steps
    
    for i = 2:n_x-1
        
        for v = 2:n_y -1
        
        kay = sqrt((p+1)*h);
        L_app = (interpol_operator(u_til(j,:),x_vec,x_vec(i) + (kay^2)*mu_fct(alpha)) - interpol_operator(u_til(j,:),x_vec,x_vec(i)))/(kay^2) ...
            + (interpol_operator(u_til(j,:),x_vec,x_vec(i) + (kay)*sig_fct(alpha)) - 2*interpol_operator(u_til(j,:),x_vec,x_vec(i)) + interpol_operator(u_til(j,:),x_vec,x_vec(i) - (kay)*sig_fct(alpha)) )/(2*(kay^2)) ; % add the diffusion term for more dimensions
        min = f_fct(alpha,x_vec(i))*h + h*L_app + u_til(j,i);
        for M = 1:2
            alpha = A_M(j,M);
        kay = sqrt((p+1)*h);
        L_app = (interpol_operator(u_til(j,:),x_vec,x_vec(i) + (kay^2)*mu_fct(alpha)) - interpol_operator(u_til(j,:),x_vec,x_vec(i)))/(kay^2) ...
            + (interpol_operator(u_til(j,:),x_vec,x_vec(i) + (kay)*sig_fct(alpha)) - 2*interpol_operator(u_til(j,:),x_vec,x_vec(i)) + interpol_operator(u_til(j,:),x_vec,x_vec(i) - (kay)*sig_fct(alpha)) )/(2*(kay^2)); % add the diffusion term
        
            if min > (f_fct(alpha,x_vec(i))*h + h*L_app + u_til(j,i) )
            min = f_fct(alpha,x_vec(i))*h + h*L_app + u_til(j,i);
            end
        end
        
        u_til(j+1,i,v) = min;
        
        end
    end
end




% plot
subplot(2,1,1)
surf(t_vec,x_vec,u_til');
xlabel('t');
ylabel('x');

subplot(2,1,2);
ezsurf(u_sol,[0, 1, -1, 1]);
