For loop not working properly

I have a problem with a loop not working correctly for some reason, or maybe it is not the loop itself, but rather something I did not pay attention to. I have been trying to solve the issue for hours now, to no avail.

Here's the code:

syms t;
syms m(t);
Nq=3;
L(1,1)=1;
L(2,1)=2;
L(3,1)=3;
w(1,1)=0.2;
w(2,1)=0.1;
w(3,1)=0.7;
G=1;

for x=1:1:2*Nq

      m0val(x)=w(1,1)*L(1,1).^(x-1) + w(2,1)*L(2,1).^(x-1) + w(3,1)*L(3,1 ).^(x-1);

end

ode=diff(m,t)==0;
cond = m(0)==m0val(1);
mf(1)=dsolve(ode,cond);

for x=2:1:2*Nq
    ode=diff(m,t)==(x-1)*G*mf(x-1);
    cond=m(0)==m0val(x);
    mf(x)=dsolve(ode,cond);

end

z=1;
for y=0.2:0.2:1
    for x=1:1:2*Nq
        mfv(z,x)=subs(mf(x),t,y);
    end

          for x=1:1:2*Nq+1
                  P(x,1)=eq(x,1);
          end
       ****for x=1:1:2*Nq+1
              if x~=2*Nq+1
                  P(x,2)=((-1).^(x-1))*mfv(z,x);
              else 
                  P(x,2)=0;
              end****
       end

          for y=3:1:2*Nq+1
              for x=1:1:2*Nq+2-y
                  P(x,y)=P(1,y-1)*P(x+1,y-2)-P(1,y-2)*P(x+1,y-1);
              end
          end
          alpha(1)=mfv(z,1);
          for x=2:1:2*Nq
              alpha(x)=P(1,x+1)/(P(1,x)*P(1,x-1));
          end
          a(1)=alpha(2);
          for x=2:1:Nq
              a(x)=alpha(2*x)+alpha(2*x-1);
          end
          for x=1:1:Nq-1
              b(x)=-(alpha(2*x+1)*alpha(2*x)).^0.5;
          end
          for x=1:1:Nq
              Jacobi(x,x)=a(x);
          end
          for x=1:1:Nq-1
              Jacobi (x+1,x)=b(x);
              Jacobi(x,x+1)=b(x);
          end
          [evec,eval]=eig(Jacobi);

                  for x=1:1:Nq
                  L(x,z+1)=eval(x,x);
                  w(x,z+1)=mfv(z,1)*evec(1,x).^2;
                  end

     z=z+1;
  end

The bit between * is where it's not working properly, because if I compute let's say P(2,2), it should be equal to (-1)^(1)*mfv(z,2) (where z=1 for the first run). It gives a value of 1 which is the value of P(1,2).

Any help would be really appreciated.

Answer 1

[per isakson]

Your code is sluggish. Preallocating variables, which the Code Analyzer proposes, would improve speed.

• I put your code in a function, because it is easier to debug functions
• I removed the "****"
• I set dbstop if error
• I started the function, which halted at an error in the function, sym/eig

Error using (line 40)
Input must not contain 'NaN' or 'Inf'.
Error in cssm (line 65)
        [evec,eval]=eig(Jacobi); 
K>> A
A =
[ 1,   0,   0]
[ 0, NaN, NaN]
[ 0, NaN, NaN]
K>>

I switch to "my function". Your code passed an input argument, which contains NaN

K>> Jacobi
Jacobi =
[ 1,   0,   0]
[ 0, NaN, NaN]
[ 0, NaN, NaN]
K>> z
z =
     1

.

"compute let's say P(2,2), it should be equal to (-1)^(1)*mfv(z,2) (where z=1 for the first run). It gives a value of 1 which is the value of P(1,2)"

Could it be that you are a victim of Matlab's smartness?

Below is an excerpt of your code.

• After the first loop P is a logical array.
• In the second loop you have P(x,2)=((-1).^(x-1))*mfv(z,x);, where LHS is logical and the RHS is <1x1 sym>. Matlab then converts the RHS to the logical value true, which Matlab display as 1. (Without warning.)

K>> ((-1).^(x-1))*mfv(z,x)
ans =
-27/10
K>> P(x,2)
ans =
     1
K>>

The excerpt of your code

        for x=1:1:2*Nq+1
            P(x,1)=eq(x,1);
        end
        for x=1:1:2*Nq+1
            if x~=2*Nq+1
                P(x,2)=((-1).^(x-1))*mfv(z,x);
            else
                P(x,2)=0;
            end
        end

[Edgard El Cham]

So how would I make the matrix P not a logical, bur rather containing numeric values?

[per isakson]

Replace

 P(x,1)=eq(x,1);

by

 P(x,1) = double( eq(x,1) );

is the short answer.

A good side-effect of preallocation of memory is the "declaration" of the class/type of the variables.