0

I am trying to implement different numerical methods in MATLAB without the use of the built-in function, such as gradient or del2. This is my code so far:

clear all
close all
x = [-1:0.1:1];
y = [-2:0.1:2];
vel = @(x,y) x+exp(-((x-x(1)).^2+(y-y(1)).^2));
nx = length(x);
ny = length(y);
derivx = zeros(nx-1,ny-1)
% The partial derivative with respect to x
for ii = 1:nx-1
    for jj = 1:ny-1
        derivx(ii,jj) = (vel(ii+1,jj) - vel(ii,jj))./(x(jj+1,ii)-x(jj,ii));
    end
end
% The partial with respect to y
derivy = zeros(ny-1,nx-1)
for ii = 1:ny-1
    for jj = 1:nx-1
    derivy(ii,jj) = (vel(ii+1,jj) - vel(ii,jj))./(y(jj+1,ii)-y(jj,ii));
    end
end

This code doesn't work with the error message saying matrix indices exceeded. 

Index in position 1 exceeds array bounds (must not exceed 1).
Error in untitled6 (line 13)
    derivx(ii,jj) = (vel(ii+1,jj) - vel(ii,jj))./(x(jj+1,ii)-x(jj,ii));

And how would i proceed to calculate the second order partials with repet to x and y (not the mixed)?

Thank you for your help in advance!

geggy
  • 13
  • 2
  • In the definition of `vel` you’re likely mixing up the vector `x` and the local variable `x`. Try changing the names of the local variables, for example `@(x1,y1)...`. – Cris Luengo Feb 17 '19 at 16:40

1 Answers1

0

To answer the question you've asked, the issue is with this: (x(jj+1,ii)-x(jj,ii). x is a vector, but you're treating it like a matrix. However, I think there are deeper issues in your code. For one, the way you are treating vel is rather unusual. You've written as a function of x and y, where x and y are presumably vectors (or matrices), but you only ever call it with scalars. If I had to make a guess, I would assume that you want to write vel as:

x = [-1:0.1:1];
y = [-2:0.1:2]'; % Note the transpose here
vel = x+exp(-((x-x(1)).^2+(y-y(1)).^2));

This will construct vel as a 2D matrix where each element (a,b) is the value of vel evaluated at x=a and y=b. Once you've done this, you can actually do away with the doubly-nested for-loops (which are almost never a good idea in MATLAB):

derivx = (vel(2:end,1:end-1) - vel(1:end-1,1:end-1)./(x(2:end)-x(1:end-1));
MrAzzaman
  • 4,734
  • 12
  • 25