MATLAB: Code’s optimization and speeding up (Legendre’s polynomials and functions)

for loopMATLABoptimizationvectorizationvectorizing

Hi. I would like to improve this function to reduce the run time and optimize it.
    function [Pl,Plm] = funlegendre(gamma)
        Plm = zeros(70,71);
        Pl = zeros(70,1);
        P0 = 1;
        Pl(1) = gamma;
        Plm(1,1) = sqrt(1-gamma^2);  
        for L=2:70
            for M=0:L
                    if(M==0)
                        if (L==2)
                        Pl(L) = ((2*L-1)*gamma*Pl(L-1)-(L-1)*P0)/L;
                        else
                        Pl(l) = ((2*l-1)*gamma*Pl(l-1)-(l-1)*Pl(l-2))/l;
                        end
                    elseif(M<L)
                        if(L==2)
                            if(M==1)
                            Plm(L,M) = (2*L-1)*Plm(1,1)*Pl(L-1);  
                            else
                            Plm(L,M) = (2*M-1)*Plm(1,1)*Plm(L-1,m-1); 
                            end
                        else
                            if(M==1)
                            Plm(L,M) = Plm(L-2,m) + (2*L-1)*Plm(1,1)*Pl(L-1);  
                            else
                            Plm(L,M) = Plm(L-2,M) + (2*L-1)*Plm(1,1)*Plm(L-1,M-1); 
                            end                 
                        end
                    elseif(M==L)
                        Plm(L,M) = (2*L-1)*Plm(1,1)*Plm(L-1,L-1);
                    else
                        Plm(L,M) = 0;
                    end   
            end
        end
        Pl = sparse(Pl);
        Plm = sparse(Plm);
        end

Best Answer

I've simplified the loop contents by moving all repeated calculations outside. Instead of checking inside the loop for M==0, M<L, M==L (there is not possible ELSE case!) is not necessary, if the M==0 case is move before the M-loop and the M==L case behind the loop. Integer indices are helpful also:
function [Pl, Plm] = funlegendreJan(gamma)
i1 = uint8(1);
i2 = uint8(2);
Plm = zeros(70,71);
Pl = zeros(70,1);
P0 = 1;
Pl(1) = gamma;
Plm(1,1) = sqrt(1 - gamma * gamma);
Plm11 = Plm(1, 1);
% L == 2:
Pl(2)     = (3*gamma*Pl(1) - P0) * 0.5;
Plm(2, 1) = 3 * Plm11 * Pl(1);
Plm(2, 2) = 3 * Plm11 * Plm11;
% L > 2:
for L = 3:70
   L1  = L - 1;
   L2  = L - 2;
   iL  = uint8(L);
   iL1 = uint8(L1);
   iL2 = uint8(L2);
   Pl(iL) = ((L + L1) * gamma * Pl(iL1) - L1 * Pl(iL2)) / L;
   C = (L + L1) * Plm11;
   Plm(iL, i1) = Plm(iL2, i1) + C * Pl(iL1);
   for iM = i2:iL1
      Plm(iL, iM) = Plm(iL2, iM) + C * Plm(iL1, iM - i1);
   end
   Plm(iL, iL) = C * Plm(iL1, iL1);
end
%Pl = sparse(Pl);
%Plm = sparse(Plm);
You did not tell us in your former thread, that Pl and Plm are sparse. This is not helpful afaics: Pl is full at all and Plm is small and >50% full. Therefore SPARSE matrices will waste computing time only.
The above version is 58% faster than the original. The non-sparse arrays will accelerate the rest of the simulation also.
Related Question