Problem with the implementation of the quadruple sum

Hello, I'm trying write code for quadruple sum:

where x and y are matrix NxN. Output matrix also is NxN. I wrote something like that:

function v = RCPM(x,y,a,b,N,mi,V,n)
  m = n - 1;
  v = zeros(N,N);
  for i = 1:m
      disp(i);
      for j = 1:m
          for k = 1:m
              for l = 1:m
                  s = mi^((a(j) - a(k + 1))^2 + (b(i) - b(l+1))^2) ...
                      * cos(V * (x * (a(j) - a(k + 1)) + ...
                      y * (b(i) - b(l+1))));
                  v = v + s;
              end
          end
      end
  end
end

But looks like it don't work well. Could anyone have a look at whether I implemented the above formula well?

[James Knowles]

My advice would be to firstly, remove the 'function' related commands. Then secondly, ensure every variable is defined correctly. For example you are using 'a' and 'b' as arrays but have not defined them as such in the above code. Then once you have it working efficiently in this manner you can read the 'function' commands if necessary.

[Walter Roberson]

I see no reason to remove the function command. function is a good thing there.

[Walter Roberson]

It is not obvious to me that the formula definitely intends algebraic matrix multiplication like you implemented by using the * operators. It seems plausible to me that element-by-element multiplication should be used instead -- but I cannot tell for sure either way.

Answer 1

[Walter Roberson]

You are iterating 1:m for each variable, which is m iterations each. However, your formula have 0 to M, which is M+1 iterations each. You define m as N-1 . It is not at all obvious that m = M+1, which would require that M = N-2

You code a(j) - a(k + 1) and b(i) - b(l+1), corresponding to alpha_lx - alpha_lx1 and beta_ly - beta_ly1 . Anyone reading this is going to assume that you got confused with the variable name lx1 as thinking that should be lx+1 but it is a completely different variable name. None of the subscripts should have +1 in them.

It is confusing that you have used i, j, k, l as the variable names controlling the loops and not ly, lx, lx1, ly1 (or Ly, Lx, Lx1, Ly1 for readability)

(a(lx)-a(lx1)) and (b(ly)-b(ly1)) occur twice in the expression, so it would probably make sense to compute them into separate variables.

For efficiency, pull expressions out of the inner loops. For example V * (x * (a(j) - a(k + 1)) does not need to be recalculated for each iteration of for l since it does not involve l