What is wrong in script?

for m1 =0:M
    T2=0;
    for m2=0:M
        T3=0;
        for s1=0:m1/2
            T4=0
            for j1=0:m1-2*s1
                T5=0
                for s2=0:(M-m1)/2
                    T6=0;
                    for j2=0:(M-m1-2*s2)
                        T6=T6+(1/(po^2))^j2;
                    end
                    T5=T5+T6*((2*1i)/(delta^(0.5)))^(M-m1-2*s2);
                end
                T4=T4+T5*(1/(po^2))^j1;
            end
            T3=T3+T4*((2*1i)/(delta^(0.5)))^(m1-2*s1);
        end
        T2=T2+T3*(factorial(M)/(factorial(m2)*factorial(M-m2)));
    end
    T1=T1+T2*(factorial(M)/(factorial(m1)*factorial(M-m1)));
end

[Torsten]

Just out of interest: Where do all these horrible sums stem from that you posted during the last days and weeks ?

What is meant if the sum goes to m1/2 or (M-m1)/2 ? Is it the usual treatment that summing up to 5/2, e.g., means summing for 0,1 and 2 ?

Further, you didn't give values for variables and loop limits. So nobody can tell what you mean by "What is wrong in Script ?" because nobody can test it.

[Athira T Das]

@Torsten It is an assignment.

M=3

w=0.2;
z=linspace(0:10000);
k = 5.9275e+06;
p=(k^2*z).^(-(3/5));

[Torsten]

And what about the loop limits ?

[Jan]

@Athira T Das: Please share the important information with the readers: Why do you assume, that the script is wrong?

[Torsten]

Is this the code you try to run ?

M = 3;
w = 0.2;
z = 0:10000;
k = 5.9275e+06;
Po = (k^2*z).^(-(3/5));
Delta = 1i*k./(2*z) + 1/w^2 + 1./Po.^2;

po = Po(1);
delta = Delta(1);
T1 = 0;
for m1 =0:M
    T2=0;
    for m2=0:M
        T3=0;
        for s1=0:m1/2
            T4=0;
            for j1=0:m1-2*s1
                T5=0;
                for s2=0:(M-m1)/2
                    T6=0;
                    for j2=0:(M-m1-2*s2)
                        T6=T6+(1/(po^2))^j2;
                    end
                    T5=T5+T6*((2*1i)/(delta^(0.5)))^(M-m1-2*s2);
                end
                T4=T4+T5*(1/(po^2))^j1;
            end
            T3=T3+T4*((2*1i)/(delta^(0.5)))^(m1-2*s1);
        end
        T2=T2+T3*(factorial(M)/(factorial(m2)*factorial(M-m2)));
    end
    T1=T1+T2*(factorial(M)/(factorial(m1)*factorial(M-m1)));
end
T1

T1 = 0

[Athira T Das]

@Torsten yes.

[Torsten]

What about the value for k ?

sqrt(0.545) or 5.9275e6 ?

What about the upper looplimit ? Like this: loop over integers up until 2 ?

k = 0;
for j = 0:5/2
  k = k + 1;
end
k

k = 3

[Athira T Das]

@Torsten

k = 5.9275e+06

And the loop limits depends on the value of M.

I am just confused about whether the for loops for the multiplication of summation is right or wrong?

[Torsten]

k = 5.9275e+06

And why do you write

rho = (0.545*z)^(-3/5) 

instead of

rho = (5.9275e+06^2*z)^(-3/5) 

in your formula ?

[Athira T Das]

Ok. But I am confused whether my for loops are working properly or not?

Answer 1

[Torsten]

I would have progrmmed it this way, but the result seems to be the same.

My guess is that there is a mathematical trick to extremely simplify the expression (something like sum_{i=0}^(i=n} nchoosek(n,i) * p^i * (1-p)^(n-i) = 1), but I don't see it at the moment.

M = 3;
w = 0.2;
z = 0:10000;
k = 5.9275e+06;
Po = (k^2*z).^(-(3/5));
Delta = 1i*k./(2*z) + 1/w^2 + 1./Po.^2;
po = Po(1);
delta = Delta(1);

X = 0;
for m1 = 0:M
  faktor = nchoosek(M,m1);
  sum1 = 0.0;
  for s1 = 0:m1/2
    faktor1 = (2*1i/sqrt(delta)).^(m1-2*s1);
    sum_inner = 0.0;
    for j1 = 0:m1-2*s1
      sum_inner = sum_inner + (1/po^2)^(j1);
    end
    sum1 = sum1 + faktor1*sum_inner;
  end
  sum2 = 0.0;
  for s2 = 0:(M-m1)/2
    faktor2 = (2*1i/sqrt(delta)).^(M-m1-2*s2);
    sum_inner = 0.0;
    for j2 = 0:M-m1-2*s2
      sum_inner = sum_inner + (1/po^2)^(j2);
    end
    sum2 = sum2 + faktor2*sum_inner;
  end  
  X = X + faktor*sum1*sum2;
end
X = X*2^M

X = 0

[Athira T Das]

@TorstenThis is the most simplified version of my actual equation. In your code the summation over m2 is missing.

[Torsten]

Since all the expressions don't depend on m2, you can take the sum out and it gives

sum_{k=0}^{k=M} nchoosek(M,k) = 2^M

[Athira T Das]

@Torsten

If any terms depends on m2, can I write the code like

for m1 = 0:M
    faktor = nchoosek(M,m1);
  for m2 = 0:M
    faktor0 = nchoosek(M,m2);
    .......................................
        ....................................
         X = X + faktor0.*sum1.*sum2;
  end
  X = X*faktor
end

[Torsten]

I think something like this:

M = 6;
w = 0.2;
z = 0:10000;
k = 5.9275e+06;
Po = (k^2*z).^(-(3/5));
Delta = 1i*k./(2*z) + 1/w^2 + 1./Po.^2;
po = Po(1);
delta = Delta(1);

X = 0;
for m1 = 0:M
  faktor_m1 = nchoosek(M,m1);
  sum_m2 = 0.0;
  for m2 = 0:M
    faktor_m2 = nchoosek(M,m2);
    sum1 = 0.0;
    for s1 = 0:m1/2
      faktor1 = (2*1i/sqrt(delta)).^(m1-2*s1);
      sum_inner = 0.0;
      for j1 = 0:m1-2*s1
        sum_inner = sum_inner + (1/po^2)^(j1);
      end
      sum1 = sum1 + faktor1*sum_inner;
    end
    sum2 = 0.0;
    for s2 = 0:(M-m1)/2
      faktor2 = (2*1i/sqrt(delta)).^(M-m1-2*s2);
      sum_inner = 0.0;
      for j2 = 0:M-m1-2*s2
        sum_inner = sum_inner + (1/po^2)^(j2);
      end
      sum2 = sum2 + faktor2*sum_inner;
    end 
    sum_m2 = sum_m2 + faktor_m2*sum1*sum2;
  end
  X = X + faktor_m1*sum_m2;
end
X

X = 2048