function MBsphere
% La fonction 'MBsphere' simule le mouvement brownien sur la sphère (ou sa
% variante quand on omet le terme de dérive pour compenser l'effet du terme
% d'Itô). Elle trace l'évolution de la distance à l'origine en fonction du
% temps.

% 'T' est la durée de simulation.
T = 4;

% 'dt' est le pas de discrétisation, et 'N' le nombre de pas à effectuer.
dt = 1/5120;
N = ceil(T/dt);

% On alloue la mémoire pour les valeurs successives du vecteur X.
les_Xt = zeros(3,N+1);
% On initialise le vecteur-colonne 'X' représentant X pour t=0, et on
% enregistre la valeur correspondante dans les_Xt.
X = [1;0;0];
les_Xt (:,1) = X;
% Méthode d'Euler stochastique. Le calcul de la matrice P(X_t) est effectué
% grâce à la sous-routine 'P' codée ci-dessous. Attention au fait qu'ici X
% suit une formalisme de vecteur-colonne, et qu'il faut donc transposer par
% rapport aux formules de l'énoncé !
for i = 1:N
    X = X + sqrt(dt) * P(X')' * randn(3,1) - 4 * X * dt;
    % Remplacer la ligne précédente par la suivante pour traiter la
    % question 3.
    %X = X + sqrt(dt) * P(X')' * randn(3,1);
    les_Xt (:,i+1) = X;
end

% 'les_Rt' est l'évolution de la distance de X_t à l'origine.
les_Rt = sqrt (sum (les_Xt.^2));
% On trace l'évolution de la distance à l'origine. Pour qu'on voie bien
% l'axe des abcisses, je fais aussi plotter une ligne sur cet axe.
plot (dt*(0:N), les_Rt, [0,T], [0,0], 'black');

return


% Sous-routine 'P'.
    function PX = P (X)
        % La fonction 'P' calcule la matrice 3*3 'PX' que la formule (3)
        % associe au 3-vecteur-ligne 'X'.
        PX = X'*X - eye(3);
        return
    end

end