Matrix increasing more than it should during iteration

I have an array (which I have called mu) and I am performing some operations on it to create a new array mu_new. I am then repeatedly iterating. I would like to make sure that the sum of mu and mu_new remains the same throughout my operations, so I have created a running sum called der to keep track of the changes between the two arrays. My problem is, after just 2 iterations the value I have for der is not the same as the value I get when I subtract the sums of the two arrays. Can anyone see why this would be the case?

Note: I removed some operations in the code to simplify it for debugging, so it is not a problem that der doesn't equal 0 here. I just want to see why der is not equal to the difference in sums of the two arrays. There are obviously some prespecified objects such as M, zeta, gx, etc.

M = 25;
mu=zeros(2*M+1,2*M+1,2*M+1);
mu(2*M+1,2*M+1,2*M+1)=1;
mu_new=mu;
zeta = rand(size(mu));
gx = rand(size(mu));

%%% iterate code below. On the second iteration der-sum(mu_new-mu,'all') is
%%% significantly differnet from 0 

for i=1:2
    der=0;
    for x=(M+1):(2*M+1)
        for y=(M+1):(2*M+1)
            for z=(M+1):(2*M+1)
        
                
                x1=M+1-(x-M-1);
                y1=M+1+(y-M-1)-(x-M-1);
                z1=M+1+(z-M-1)-(x-M-1);
        
                
                x2=M+1+(x-M-1)-(y-M-1);
                y2=M+1-(y-M-1);
                z2=M+1+(z-M-1)-(y-M-1);
        
                x3=M+1+(x-M-1)-(z-M-1);
                y3=M+1+(y-M-1)-(z-M-1);
                z3=M+1-(z-M-1);
        
                leave=gx(x,y,z)+zeta(x,y,z)+gx(x1,y1,z1)+zeta(x1,y1,z1)...
                    +gx(x2,y2,z2)+zeta(x2,y2,z2)+gx(x3,y3,z3)+zeta(x3,y3,z3);
        
                mu_new(x,y,z)=(1-leave)*mu(x,y,z);
                der=der-leave*mu(x,y,z);
        
                if x<2*M+1 && y<2*M+1 && z<2*M+1
                    mu_new(x+1,y+1,z+1)=mu(x+1,y+1,z+1)...
                        +(gx(x,y,z)+zeta(x,y,z))*mu(x,y,z);
                    der=der+(gx(x,y,z)+zeta(x,y,z))*mu(x,y,z);
                end
        
                if x==(M+1) && y<2*M+1 && z<2*M+1
                   mu_new(x1+1,y1+1,z1+1)=mu(x1+1,y1+1,z1+1)...
                       +(gx(x1,y1,z1)+zeta(x1,y1,z1))*mu(x,y,z);
                   der=der+(gx(x1,y1,z1)+zeta(x1,y1,z1))*mu(x,y,z);
                else
                    mu_new(x-1,y,z)=mu(x-1,y,z)...
                        +(gx(x1,y1,z1)+zeta(x1,y1,z1))*mu(x,y,z);
                    der=der+(gx(x1,y1,z1)+zeta(x1,y1,z1))*mu(x,y,z);
                end
        
            end
        end
    end
    der-sum(mu_new-mu,'all')
    mu_new=mu;
end

ans = 0

ans = 0

%der-sum(mu_new-mu,'all')

[Torsten]

I get 0 as result (see above).

Maybe your arrays or functions gx and zeta return Inf or NaN values ?

[Mahmood Haddara]

It is 0 the first time you run it but after iterating once (updating mu and rerunning) you will see it becomes nonzero (see above)

I checked and there is no Inf or NaN.

[Torsten]

Works. See above.

[Mahmood Haddara]

Very strange. Thanks for your help!

[Torsten]

In your original code you set

mu = mu_new

instead of

mu_new = mu

at the end of your code.

In this case,

der-sum(mu_new-mu,'all')

was not 0 in the second run.

I must admit that I cannot understand why.

[Mahmood Haddara]

Very interesting. Thanks for spotting that!

Answer 1

[Yogesh]

Hi @Mahmood Haddara

With the information you have provided above, I have verified your code for M in range [1, 100] and its working as expected. If there is a problem in successive runs make sure you are resetting all variables to required values as they will acquire some value in workspace.

Hope that helps!

[Mahmood Haddara]

Thanks!

[Mahmood Haddara]

Hi Yogesh,

As Torsten mentioned above, if we set

mu = mu_new

At the end instead of

mu_new = mu

We get der not equal to 0 on the second run. Do you have any idea why this is the case?