function [U, Umax, Xmax] = theBigFish(alpha)


% Cours de méthodes numériques

% MATLAB : Problème 8

% Author : Lange Joseph

% Date : 20/12/2011

% Version : 2





%fonction qui "calcule la solution du

%problème du cube. Le paramètre alpha contiendra le nombre d'intervalles requis pour la discrétisation

%d'un quart d'une arête du cube. Le nombre de noeuds sur une arête sera donc n = 4*alpha+1. La

%fonction retournera l'ensemble U des inconnues dans un tableau a 3 dimensions. La valeur maximale

%sera continue dans Umax, tandis que le vecteur de taille de trois Xmax contiendra les coordonnées du

%noeud où ce maximum est atteint."(voir énoncé)

%J'ai optimisé la méthode donnée dans MyNiceProgram.m de 2 façons :

%1-On crée les diagonales de la sparse matrix avant de la créer

%2-On peut en fait décomposer le problème en différents problèmes

%identiques indépendant vu la symétrie du problème. On peut en fait

%décomposer le cube en six petites pyramides tel que chacune d'entre elle

%n'a aucune raison d'avoir une solution différente de toutes les autres

%pyramides. Cela permet donc d'économiser du temps lors de la résolution du

%système d'équations linéaire. Il y a moyen d'économiser de la mémoire

%dans mon programme en changeant la façon dont l'index est fait de manière

%à diminuer par 36 la taille de la matrice A en n'indexant que les points dans la pyramide.

%Une fois qu'on a trouvé la solution pour une petite pyramide on fait des copies de cette pyramide à

%l'arrache dans tous les sens pour retrouver la solution complète :)



na=4*alpha+1;



n = na^3;

h = 2/(na-1);



uns = ones(n, 1);

zero = zeros(n, 1);



diag0 = uns;

diagim = zero;

diagip = zero;

diagjm = zero;

diagjp = zero;

diagkm = zero;

diagkp = zero;



B = zero;





for i = 2:na-1

    for j = 2:na-1

        for k=2:na-1

            if i<= j && k<=i

                index = i + (j-1)*na +(k-1)*na^2;

                diag0(index) = 6;

                diagip(index+1) = -1;

                diagim(index-1) = -1;

                diagjp(index+na) = -1;

                diagjm(index-na) = -1;

                diagkp(index+na^2) = -1;

                diagkm(index-na^2) = -1;

                if i>=(na+1)/2 && j>=(na+1)/2 && k>=(na+1)/2,

                    B(index) = 10;

                else

                    B(index) = 0;

                end

                %Les if en dessous servent à faire la jonction en la

                %pyramide qu'on est en train de calculer et les autres

                %pyramides(on connait pas les points dans les autres

                %pyramides donc on prend les points dans cette pyramide ci

                %qui correspondent aux points qu'on veut dans les autres

                %pyramides)

                %En gros, on dit que sur le bord de cette pyramide, on

                %remet cette même pyramide dans un autre sens

                if i==j

                    %Pour les i

                    diagip(index+1) = 0;

                    diagjp(index+na) = -2;

                    %Pour les j

                    diagjm(index-na) = 0;

                    diagim(index-1) = -2;

                end

                if k==i

                    %Pour les k

                    diagkp(index+na^2)=0;

                    diagip(index+1) = -2;

                    %Pour les i

                    diagim(index-1) = 0;

                    diagkm(index-na^2) = -2;

                end

                if i==j && k==i

                    diagip(index+1) = 0;

                    diagkm(index-na^2) = -3;

                    diagjp(index+na) = -3;

                end

            end

        end

    end

end





A = spdiags([diagkm diagjm diagim diag0 diagip diagjp diagkp],[-na^2 -na -1 0 1 na na^2], n, n);



B = B.*(h*h);



%Solve linear system

U = A\B;





U = reshape(U, [na na na]);



%On recolle les petits bouts

%On reconstitue la première moitié

for i = 2:na-1

    for j = 2:na-1

        for k=2:na-1

            if k >= i

                U(i, j, k) = U(k, j, i);

            end

            if i>= j

                U(i, j, k) = U(j, i, k);

            end

        end

    end

end

%On reconstitue l'autre moitié

for i = 2:na-1

    for j = 2:na-1

        for k=2:na-1

            if k >= i

                U(i, j, k) = U(k, j, i);

            end

            if i>= j

                U(i, j, k) = U(j, i, k);

            end

        end

    end

end





[Umax, position] = max(U(:));



[X1 X2 X3] = ind2sub(size(U),position);

Xmax = [X1-1 X2-1 X3-1]*h-[1 1 1];





end