An elementwise multiplication is performed by the .* operator, while * is a matrix multiplication. This should solve your problem already. Please post your trial to vectorize the code instead of letting us create it from scratch.
An optimization of code depends on the input values. It matters if Nx_max is 10 or 1e7. It matters if dxy and gamma are arrays or functions. So please post a version of your code, which can be started by copy&paste. Otherwise suggestions contains a lot of guessing.
But if guessing is required, consider the standard tips fopr loop optimization at first:
- avoid repeated calculations
- pre-allocate the output
E.g. (1/(2*h^2)) is calculated nearly Ny_max * Nx_max times. So better create a variable before the loop to carry the result. So befpre creating the vectorized version, measure the speed with a cleaner version like this:
hs = h^2;
h2 = 2 * h;
hs2 = 2 * hs;
c1 = 1/hs2;
g2 = 2 * gamma;
lp = zeros(Nx_max-1, Ny_max-1, t1);
for x = 2:Nx_max-1
c2 = (x-Nx-1)/g2;
for y = 2:Ny_max-1
lp(x,y,t1) = c1 * dxy(x,y) * p(x+1,y+1,t1) + ...
(dxx(x,y) / hs - vx(x,y) / h2 - c2) * p(x+1,y,t1) - ...
c1 * dxy(x,y) * p(x+1,y-1,t1) + ...
(dyy(x,y)/hs - vy(x,y)/h2 - (y-Ny-1)/g2) * p(x,y+1,t1) - ...
(2 * dxx(x,y) + 2 * dyy(x,y)) * p(x,y,t1) / hs + ...
(dyy(x,y)/hs + vy(x,y)/h2 +(y-Ny-1)/g2) * p(x,y-1,t1) - ...
dxy(x,y) / hs2 * p(x-1,y+1,t1) + ...
(dxx(x,y)/hs + vx(x,y) / h2 + c2) * p(x-1,y,t1) + ...
dxy(x,y) / hs2 * p(x-1,y-1,t1);
end
end
Then, in a next step, replace the loop over x by replacing all "x" by "(2:Nx_max-1)" and the operators "*" and "/" by their elementwise versions ".*" and "./". Test, if the result is still not changed except for rounding errors. Fix "(2:Nx_max-1)+1" to the faster "3:Nx_max", because in the 2nd version the index vector is not calculated exlicitly.
Finally you can try to replace the inner loop also, but perform a speed test afterwards. A complete vectorization is not necessarily the fastest version.
