Calculating second centred difference for FDM - indices for iteration

I have an NxN matrix 'phi' that is a function of x and y, and is constantly getting updated over iteration (is a function of t as well). To calculate the Laplacian (sum of spatial second derivatives), I calculate the second centred difference, defined as :

I am attaching the codes I have been given in two cases: first when phi is just a function of x (1D) and second when it is a function of x and y (2D). In both cases, my doubt is why pind and nind are defined this way (pind and nind refer to the positive index - forward step in x or y dirn - and negative index - backward step - respectively). I believe it is just a question of me not understanding the syntax used.

pind=[2:N,N];
nind=[1,1:N-1];

for i=1:length(c)
    

   dphidx = (phi(pind)-phi(nind))/(2*h);
   d2phi = (phi(pind)+phi(nind)-2*phi)/(h^2);
end

code i) 1 dimensional example

pind=[2:N,1];
nind=[N,1:N-1];

%iteration begins
for iter = 1:Niter,
    
dphidx = (phi(:,pind)-phi(:,nind))/(2*h);
dphidy = (phi(pind,:)-phi(nind,:))/(2*h);
d2phi = (phi(:,pind)+phi(:,nind)+phi(pind,:)+phi(nind,:)-4*phi)/(h^2);
end

code ii) 2 dimensional example

I have checked test outputs to see how pind and nind change over the iterations, but I still don't understand why the definition is done like this. Can attach screenshots of that test output if required.

Answer 1

[darova]

i think something is wrong here

pind = 2 3 4 5
nind = 1 1 2 3

correct way for 1D would be:

for i = 2:length(phi)-1
    dphidx = (phi(i-1)-phi(i+1))/2/h;
    d2phidx = (phi(i-1)+2*phi(i)-phi(i+1))/h^2;
end

2d:

for i = 2:size(phi,1)-1
    for j = 2:size(phi,2)-1
        dphidx = (phi(i-1,j)-phi(i+1,j))/2/h;           % dphi/dx
        dphidy = (phi(i,j-1)-phi(i,j+1))/2/h;           % dphi/dy
        jj = j-1:j+1;
        dphidx1 = (phi(i-1,jj)-phi(i+1,jj))/2/h;
        d2phi = (dphidx1(3) - dphidx1(1))/2/h;          % d2phi/dx/dy
    end
end

[Pranav Krishnan]

Hey, thank you so much for your answer. But I don't fully understand what you've done to calculate d2phi..

• Why is it necessary to define a separate matrix dphidx1? What does it mean and how does it look, since passing jj to it gives the terms j-1, j and j+1 to it?
• What does it mean to call dphidx1(3)? How is one paramter being passed to an nxn matrix? Is it 3rd row? Why do 3 and 1 give d2phi?
• Does this account for boundary conditions? (intuitively it should, because every row starts from the second element, up till the second last; likewise for columns. Just confirming though, since I think the assumption in this program is that my NxN grid defined here is periodic. That is, it has continuous boundaries

[darova]

Here is the formula

dphidx1 is . It has length of 3 elements.

Here is a simple scheme: