Trying to define an array as the derivative of the terms of another array gives just zeros.

I have this code. Most of it is working fine. The problem is when I'm trying to define the array dpx and dpy. Instead of giving me the correct values it spits out all zeros once we get to the second time iteration. It works fine on the first though. I'm not sure if I might just be indexing something incorrectly. Since I don't know what is causing the mistake, I'll attach the full code:

N = 4;
Nx = 3;
Ny = 3;
h = 0.2;
t = 0.00005;
n_max = 20000;
Nx_max = 2*Nx+1;
Ny_max = 2*Ny+1;
z = 2/3;
gamma = 4/3;
d = 0.2;
f = zeros(Nx_max,Ny_max,5);
sf=zeros(Nx_max, Ny_max,5,Nx_max);
sf1=zeros(Nx_max,Ny_max,5);
p = zeros(Nx_max,Ny_max, 5);
dxx = zeros(Nx_max, Ny_max);
dxy = zeros(Nx_max, Ny_max);
dyy = zeros(Nx_max, Ny_max);
Ix = zeros(Nx_max,Ny_max, 5);
Iy = zeros(Nx_max,Ny_max,5);
Vx = zeros(Nx_max,Ny_max,5);
Vy = zeros(Nx_max,Ny_max, 5);
lp = zeros(Nx_max,Ny_max,5);
lf = zeros(Nx_max,Ny_max,5);
dpx = zeros(Nx_max,Ny_max,5);
dpy = zeros(Nx_max,Ny_max,5);
plf = zeros(1,5);
pf = zeros(1,5);
t2=2;   %indexing so that we only save last few iterations
t1=1;%also used for indexing

for n =1:5
    for x  = 1:Nx_max;
        for y = 1:Ny_max;         

              %Defining terms and initial value of function

              dxx(x,y) = (h^z)*(x^2+y^2)^(1/3)* ((d-1+z)-z*x^2/(x^2+y^2));
              dyy(x,y) = (h^z)*(x^2+y^2)^(1/3)* ((d-1+z)-z*y^2/(x^2+y^2));
              dxy(x,y) = (h^z)*(x^2+y^2)^(1/3)* (-z*x*y/(x^2+y^2));
              p(x,y,1) = x^2+3*y;
              f(x,y,1) = y^3+x^4 ;

          end
      end

%This next loop was pulled out because I needed terms like p(x+1) which isn't known until the whole x loop is run through at least once (I think)

for x = 1:Nx_max
    for y = 1:Ny_max
        dpx(1,y,t1) = 0;
        dpx(Nx_max,y,t1) = 0;
        dpx(x,1,t1) = 0;
        dpx(x,Nx_max,t1) = 0;
        dpy(1,y,t1) = 0;
        dpy(Nx_max,y,t1) = 0;
        dpy(x,1,t1) = 0;
        dpy(x,Nx_max,t1) = 0;
        if x>1 && x<Nx_max   %So that we don't run into an index being zero
           dpx(x,y,t1) = p(x+1,y,t1) - p(x-1,y,t1);
        end
        if y> 1 && y<Ny_max    %Same as the x reasoning
           dpy(x,y,t1) = p(x,y+1,t1)-p(x,y-1,t1);
        end
        Ix= f.*dpx/2;
        Iy= f.*dpy/2;
    end
end

          Vx = (N-1)*h*(convn(Ix,dxx)+convn(Iy,dxy));
          Vy = (N-1)*h*(convn(Ix,dxy)+convn(Iy,dyy));

%This section seems to have the same issue as the the dpx and dpy's

          for x = 2:Nx_max-1
              for y = 2:Ny_max-1
                  lp(x,y,t1) = (1/(4*h^2))*dxx(x,y)*p(x+1,y+1) + (dxx(x,y)/h^2 - Vx(x,y)/(2*h) - x/(2*gamma))* p(x+1,y,t1) - ...
                      (1/(4*h^2))*dxy(x,y)*p(x+1,y-1,t1) + (dyy(x,y)/(h^2) - Vy(x,y,t1)/(2*h) - y/(2*gamma))*p(x,y+1,t1) - ...
                      (2*dxx(x,y) + 2*dyy(x,y))*(1/h^2)*p(x,y) +(dyy(x,y)/(h^2) + Vy(x,y)/(2*h) +y/(2*gamma))*p(x,y-1) - ...
                      (dxy(x,y)/(4*h^2))*p(x-1,y+1) +(dxx(x,y)/(h^2) +Vx(x,y)/(2*h) +x/(2*gamma))*p(x-1,y,t1) + (dxy(x,y)/(4*h^2))*p(x-1,y-1,t1);
                  lf(x,y,t1)= (1/(4*h^2))*dxx(x+1,y+1)*f(x+1,y+1,t1) + (dxx(x+1,y)/h^2 + Vx(x+1,y,t1)/(2*h) + (x+1)/(2*gamma))* f(x+1,y,t1) - ...
                      (1/(4*h^2))*dxy(x+1,y-1)*f(x+1,y-1,t1) + (dyy(x,y+1)/(h^2) + Vy(x,y+1,t1)/(2*h) +(x+1)/(2*gamma))*f(x,y+1,t1) - ...
                      (2*dxx(x,y) + 2*dyy(x,y))*(1/h^2)*f(x,y,t1) +(dyy(x,y-1)/(h^2) - Vy(x,y-1,t1)/(2*h) -(y-1)/(2*gamma))*f(x,y-1,t1) - ...
                      (dxy(x-1,y+1)/(4*h^2))*f(x-1,y+1,t1) +(dxx(x-1,y)/(h^2) -Vx(x-1,y,t1)/(2*h) -(x-1)/(2*gamma))*f(x-1,y,t1) + (dxy(x-1,y-1)/(4*h^2))*f(x-1,y-1,t1);

                  lp(1,y,t1) = 0;
                  lp(Nx_max,y,t1)= 0;
                  lp(x,1,t1) = 0;
                  lp(x,Ny_max,t1) = 0;
                  lf(1,y,t1) = 0;
                  lf(Nx_max,y,t1)= 0;
                  lf(x,1,t1) = 0;
                  lf(x,Ny_max,t1) = 0;

              end
          end
          plf = sum(sum(p.*lf));
          pf = sum(sum(p.*f));

%Next we use all the values found to get the next time iteration

          for x = 1:Nx_max
              for y = 1:Ny_max
                  p(x,y,t2)= t*(lp(x,y,t1) - plf(t1)*p(x,y,t1)/pf(t1)) + p(x,y,t1);
                  f(x,y,t2) = t*(lf(x,y,t1) - plf(t1)*f(x,y,t1)/pf(t1)) + p(x,y,t1);
              end
          end

      if t2>5 %tells it to loop back to first spot in the array. This should conserve on memory.
          t2=1;
      end
  end

Answer 1

[Geoff Hayes]

The following code is where the dpx and dpy matrices are set

 if x>1 && x<Nx_max   %So that we don't run into an index being zero
     dpx(x,y,t1) = p(x+1,y,t1) - p(x-1,y,t1);
 end
 if y> 1 && y<Ny_max    %Same as the x reasoning
     dpy(x,y,t1) = p(x,y+1,t1)-p(x,y-1,t1);
 end

At the end of the script/function, if you print the results for dpx and dpy, then like you said only dpx(:,:,1) and dpy(:,:,1) have non-zero values.

As per the code, the third dimension in each of these arrays uses the tl variable which is initialized to 1 but is never incremented. So each iteration overwrites the last set of data and the matrices do not get set properly.

The trick will be how to decide when to increment t1. Maybe an outer loop is missing to iterate over t1 from 1:5?

NOTE that there may be a similar problem with the p and f matrices. They are both defined to be 7x7x5 but only the first 7x7 matrix is set:

 p(x,y,1) = x^2+3*y;
 f(x,y,1) = y^3+x^4 ;

Perhaps n should be used instead of 1?

[sasha]

Yes, thanks for that I caught that myself a moment ago.

There should be one last line that reads:

t1 = t2
if t2<5
   t2 = t2+1
else
   t2 = 1
end

Right now I have the max n as 5, but I'm going to possibly need very many more iterations. Hence the t's