Combining 3 for loops into 1 in Matlab

I have this section of code which is working fine but i don't like the way i have implemented those 3 separate for loops. Can anyone suggest to me that how can i merge them together to make it even more efficient and compact ? Thank you

    clear variables;
    close all;
    clc;

    ilambda=5;
    mu=2;

    n=2;
    jmax=n+1;
    P= sym('P',[jmax,jmax]);

    for j1 = 1:jmax
            for j2 = 2:jmax-1
                [c1,c2,c3,c4,c5]=coefficients(ilambda,mu,j1,j2);    
                if j1<jmax
                    E(j1, j2) = c1*P(j1, j2) - c2 * P(j1+1, j2) - c3 * P(j1, j2+1) - c4 * P(j1, j2-1);
                else
                    E(j1, j2) = c1 * P(j1, j2) - c3 * P(j1, j2+1) - c4 * P(j1, j2-1);
                end

            end
    end

        j2=1;
        for j1=1:jmax;
            [c1,c2,c3,c4,c5]=coefficients(ilambda,mu,j1,j2);    
            if (j1<jmax)
                E(j1, j2) = c1*P(j1, j2) - c2 * P(j1+1, j2) - c3 * P(j1, j2+1);
            else
                E(j1, j2) = c1*P(j1, j2) - c3 * P(j1, j2+1);
            end
        end

        j2=jmax;
        for j1=1:jmax;
            [c1,c2,c3,c4,c5]=coefficients(ilambda,mu,j1,j2); 
            if (j1==1)
                E(j1, j2) = c1*P(j1, j2) - c2 * P(j1+1, j2) - c4 * P(j1, j2-1);
            elseif (j1==jmax)
                c1=((j1-1)*mu)+((j2-1)*mu);
                E(j1, j2) = c1 * P(j1, j2) - c4 * P(j1, j2-1) - c5 * P(j1-1, j2);
            else
                E(j1, j2) = c1*P(j1, j2) - c2 * P(j1+1, j2)- c4 * P(j1, j2-1)- c5 * P(j1-1, j2);
            end
        end

[Star Strider]

With the conditionals, it’s going to be difficult to vectorise your code. Otherwise, I would use ‘J1’ and ‘J2’ matrices derived from meshgrid(j1,j2). Your conditionals make that approach difficult.

Preallocating is always good, but the sym call is unnecessary. Instead, use:

P = zeros(jmax);
P = nan(jmax);

or something similar.

Note that giving one argument to those functions creates a (jmax x jmax) matrix of them.

[Aftab Ahmed Khan]

Ok done. If jmax is very large. Let 2601, then the following loop takes too much time. Can i implement it in a different way ?

      index=1;
      % Build up coefficient matrix
      for ix=1:jmax
          for iy=1:jmax
          [a,b]=coeffs(E(ix,iy));
          A(index,ismember(allvars,b))=a;
          index=index+1;
          end
      end

[Star Strider]

Unfortunately, I doubt it, because coeffs cannot be vectorised (to the best of my knowledge), since the Symbolic Math Toolbox isn’t intended for such calculations.

It would be nice for you to be able to get away from the Symbolic Math Toolbox for your calculations, but I can’t figure out how to numerically get the coefficients you want. The best I can come up with is that for a statement such as:

E(j1, j2) = c1*P(j1, j2) - c2 * P(j1+1, j2) - c4 * P(j1, j2-1);

creating a coefficient matrix ‘C’ (for want of originality this morning), you could define it as:

C(j1, j2)   =  c1;
C(j1+1, j2) = -c2;
C(j1, j2-1) = -c4;

and so for the others. That would probably work providing you don’t overwrite them with other values later. That could eliminate the later coeffs call and the computational overhead of the Symbolic Math Toolbox, since your ‘E’ matrix just seems to be a way to store the coefficients anyway.

This just a guess. You might code it as an experiment to see if it gives you the results you want.

[Aftab Ahmed Khan]

Hi, I think your guess is right. The way i am approaching this problem is like this but i can't map those coefficients in a right order. I want to map them like the table given below.

clear variables;
close all;
clc;

jmax=3;
mu=2;
ilambda=5;

      C=zeros(jmax^2);
          for j1=1:jmax
              for j2=1:jmax
                  switch j2
                      case 1                                                                  
                              c1=ilambda+((j1-1)*mu)+((j2-1)*mu);
                              c2=(((j1-1)+1)*mu);
                              c3=(((j2-1)+1)*mu);
                              c4=ilambda;
                              c5=ilambda;

                          if (j1<jmax)
                              C(j1,j2)=c1;
                              C(j1+1,j2)=c2;
                              C(j1,j2+1)=c3;
                          else
                              C(j1,j2)=c1;
                              C(j1,j2+1)=c3;
                          end

                      case jmax                                                               
                              c1=ilambda+((j1-1)*mu)+((j2-1)*mu);
                              c2=(((j1-1)+1)*mu);
                              c3=(((j2-1)+1)*mu);
                              c4=ilambda;
                              c5=ilambda; 
                          if (j1==1)
                              C(j1,j2)=c1;
                              C(j1+1,j2)=c2;
                              C(j1,j2-1)=c4;
                          elseif (j1==jmax)
                              c1=((j1-1)*mu)+((j2-1)*mu);
                              C(j1,j2)=c1;
                              C(j1,j2-1)=c4;
                              C(j1-1,j2)=c5;
                          else
                              C(j1,j2)=c1;
                              C(j1+1,j2)=c2;
                              C(j1,j2-1)=c4;
                              C(j1-1,j2)=c5;
                          end

                      otherwise
                              c1=ilambda+((j1-1)*mu)+((j2-1)*mu);
                              c2=(((j1-1)+1)*mu);
                              c3=(((j2-1)+1)*mu);
                              c4=ilambda;
                              c5=ilambda;
                          if j1<jmax
                              C(j1,j2)=c1;
                              C(j1+1,j2)=c2;
                              C(j1,j2+1)=c3;
                              C(j1,j2-1)=c4;
                          else
                              C(j1,j2)=c1;
                              C(j1,j2+1)=c3;
                              C(j1,j2-1)=c4;
                          end
                  end
              end
          end

[Star Strider]

The matrix you created looks as though it will work. It may seem awkward code, but if it maps your coefficients correctly, and is faster and more efficient than using the Symbolic Math Toolbox, that’s likely the best you can hope for. It solves the problem of getting the coefficients in the correct order, because you’re mapping them directly to your matrix. Conditional statements are difficult to make into efficient, vectorised code.

If it works, go with it!