Частные производные первого и второго порядка функции с двумя переменными в MATLAB

Question

Я пытаюсь реализовать различные численные методы в MATLAB без использования встроенной функции, такой как градиент или del2. Это мой код до сих пор:

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

Этот код не работает с сообщением об ошибке о превышении индексов матрицы.

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));

И как мне рассчитать партиалы второго порядка с повторением по x и y (не смешанным)?

Заранее спасибо за помощь!

Answer 1

Чтобы ответить на вопрос, который вы задали, проблема заключается в следующем: (x(jj+1,ii)-x(jj,ii), x это вектор, но вы рассматриваете это как матрицу. Тем не менее, я думаю, что в вашем коде есть более глубокие проблемы. С одной стороны, как вы относитесь к vel довольно необычно. Вы написали как функцию x а также y, где x а также y предположительно векторы (или матрицы), но вы когда-либо называете это только скалярами. Если бы мне пришлось сделать предположение, я бы предположил, что вы хотите написать vel как:

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));

Это построит vel в качестве 2D-матрицы, где каждый элемент (a,b) является значением vel оценивается в x=a а также y=b, Как только вы это сделаете, вы можете покончить с дважды вложенными циклами for (которые в MATLAB почти никогда не бывают хорошими):

derivx = (vel(2:end,1:end-1) - vel(1:end-1,1:end-1)./(x(2:end)-x(1:end-1));