Simple Math With Complex Arrays

I seem to be unable to do some basic arithmetic with complex numbers in matlab.

When I define the arrays,

 B=complex(cos(A),sin(A));
 C=complex(cos(A),sin(A));

I get complex arrays. A is real. When I then define a pair of vectors,

i=2:M-1;
j=2:N-1;

To make the following calculation,

C(i,j) = (B(i-1,j)-B(i+1,j))/(2*dx)

or this,

C(i,j)= complex(real((B(i-1,j)-B(i+1,j))/(2*dx)),imag((B(i-1,j)-B(i+1,j))/(2*dx)));

C(i,j) is real. This is just subtraction. Does there exist some way to do basic arithmetic with complex arrays? Or do I need to break these arrays up and then define my own algebra?

[James Tursa]

Show an entire running code snippet with A, M, N, dx etc all defined that replicates your problem. That way we can run exactly what you are running.

To answer your question, you don't need to do anything special to do complex arithmetic in MATLAB. All of the operators etc are already coded to handle complex variables.

[Matthew Reed]

There is something bizarre going on here. When I use a simplified version, it works,

A(1:250,1:250)=20;

[M,N]=size(A);

dx=1;

B=complex(cos(A),sin(A));
C=complex(cos(A),sin(A));
i=2:M-1;
j=2:N-1;
C(i,j) = (B(i-1,j)-B(i+1,j))/(2*dx);

As one expects, these lines of code, copied individually from the .m file without renaming, work.

I1 is a double precision array of real numbers produced by another bit of code, choose any 250x250 array. It never changes.

[M,N]=size(I1);
dx=1;
Phinot = I1; 
i=2:M-1;
j=2:N-1;
Phinot2=complex(cos(Phinot),sin(Phinot));
GradxPhi2=complex(cos(Phinot),sin(Phinot));
GradxPhi2(i,j)= (Phinot2(i-1,j)-Phinot2(i+1,j))/(2*dx);

If, however, I run the entire previous body of the code, it works through these lines:

GradxPhi2=complex(cos(Phinot),sin(Phinot));
GradyPhi2=complex(cos(Phinot),sin(Phinot));

At this point GradxPhi2 is complex. But as soon as it runs this line:

 GradxPhi2(i,j)= (Phinot2(i-1,j)-Phinot2(i+1,j))/(2*dx);

it turns back into a real array. This only seems to happen when I've run the other, irrelevant code in the .m file.

There isn't enough room to post all of it. The breaking point is line 108.

[Matthew Reed]

This is a self-contained example where the code fails and produces a real GradxPhi2 from an imaginary Phinot2.

lambda=.000000532;
I1(1:250,1:250)=0;
I2(1:250,1:250)=.1;

dx=1;
dz=1;
NA=.1;
k=2*pi/lambda;
[M,N]=size(I1);
L=dx*M;
W=dx*N;

Phinot = I1;
Phinot(:,:) = 0; %an assumed phase before iteration.
DeltaI = (I1-I2)/dz;

%note the index grids
 kkkk=2:M-1;
  nnnn=2:N-1;

   GradxI1 = 0.*I1;
   GradyI1 = 0.*I1;
   GradxPhi= 0.*Phinot;
   GradyPhi= 0.*Phinot;
   GradsquaredPhi = 0.*Phinot;

   corrector = dx*max(max(DeltaI))/(2*sqrt(2)*NA);

   Phinot2 = complex(cos(Phinot),sin(Phinot));

   Phinot2=complex(cos(Phinot),sin(Phinot));
   GradxPhi2=complex(cos(Phinot),sin(Phinot));
   GradyPhi2=complex(cos(Phinot),sin(Phinot));

   GradxPhi2(kkkk,nnnn)= (Phinot2(kkkk-1,nnnn)-Phinot2(kkkk+1,nnnn))/(2*dx);

Answer 1

[Matthew Reed]

This is part of a bizarre feature that Matlab uses to save memory.

This 'bug' will always occur when an arrays values are all 0+0*i.

For instance, when troubleshooting a loop, if that loop, internally, uses complex algebra to do things like take cross products, it will look like it's malfunctioning if the initialization values in the loop are 0+0*i. It will, however, start producing complex numbers, and the complex algebra will work properly, on the second iteration if the first iteration's output is nonzero.

[James Tursa]

So, you have discovered the value of creating a small working subset of code that reproduces the problem. Had you posted that example at the outset, we could have seen right away the 0+0i issue and advised you what the issue was. MATLAB routinely converts the results of operations and function calls to strict real if the imaginary part is all 0's.

Answer 2

[James Tursa]

So, MATLAB has functions i and j defined as the imaginary sqrt(-1) values. I notice that you are using i and j as indexing. This will shadow these MATLAB i and j functions. If you have any code that uses i and j assuming they are the sqrt(-1) value but instead get your indexing vectors, this will of course will not work. I would recommend changing your indexing variable names to something else, like x and y or maybe I and J so that there isn't a potential conflict. And rewrite any code you wrote that uses i and j for sqrt(-1) to instead use 1i and 1j.

[Matthew Reed]

Just did that, used matlab's nice little variable replacer in the bulk code. It did not fix the problem. i and j as vectors in the minimal example worked just fine.