векторизации MATLAB функция

Question

У меня двойное суммирование по м = 1: M и п = 1: N для полярной точке с координатами ро, фи, г:

Я написали векторизованные обозначения:

N = 10; 
M = 10; 
n = 1:N; 
m = 1:M; 

rho = 1; 
phi = 1; 
z = 1; 

summ = cos (n*z) * besselj(m'-1, n*rho) * cos(m*phi)'; 

Теперь мне нужно переписать эту функцию для приема векторов (столбцов) координат Rho, фи, г. Я попробовал arrayfun, cellfun, просто для цикла - они работают слишком медленно для меня. Я знаю о «советах и ​​трюках манипулирования массивами MATLAB», но как новичок MATLAB я не могу понять repmat и другие функции.

Может ли кто-нибудь предложить векторизованное решение?

Answer 1

Я думаю, что ваш код уже хорошо векторизация (для n и m). Если вы хотите, чтобы эта функция также принимала массив значений rho/phi/z, я предлагаю вам просто обрабатывать значения в цикле for, поскольку я сомневаюсь, что дальнейшая векторизация принесет значительные улучшения (плюс код будет сложнее читать).

Сказав, что в приведенном ниже коде я попытался выделить векторную часть, в которой вы вычисляете различные компоненты, которые умножаются на {row N} * { matrix N*M } * {col M} = {scalar}, делая один вызов функций BESSELJ и COS (я помещаю каждый из строк/матриц/столбец в третьем измерении). Их размножение еще done in a loop (ARRAYFUN быть точным):

%# parameters 
N = 10; M = 10; 
n = 1:N; m = 1:M; 

num = 50; 
rho = 1:num; phi = 1:num; z = 1:num; 

%# straightforward FOR-loop 
tic 
result1 = zeros(1,num); 
for i=1:num 
    result1(i) = cos(n*z(i)) * besselj(m'-1, n*rho(i)) * cos(m*phi(i))'; 
end 
toc 

%# vectorized computation of the components 
tic 
a = cos(bsxfun(@times, n, permute(z(:),[3 2 1]))); 
b = besselj(m'-1, reshape(bsxfun(@times,n,rho(:))',[],1)');    %' 
b = permute(reshape(b',[length(m) length(n) length(rho)]), [2 1 3]); %' 
c = cos(bsxfun(@times, m, permute(phi(:),[3 2 1]))); 
result2 = arrayfun(@(i) a(:,:,i)*b(:,:,i)*c(:,:,i)', 1:num);   %' 
toc 

%# make sure the two results are the same 
assert(isequal(result1,result2)) 

Я сделал еще один эталонный тест, используя функцию TIMEIT (дает более справедливые тайминги). Результат согласуется с предыдущим:

0.0062407 # elapsed time (seconds) for the my solution 
0.015677  # elapsed time (seconds) for the FOR-loop solution 

Обратите внимание, что по мере увеличения размера входных векторов, оба метода начнут иметь аналогичные тайминги (для-петли даже выигрывает в некоторых случаях)

Answer 2

Вам нужно создать две матрицы, скажем m_ и n_ так, что при выборе элемента i,j каждой матрицы вы получаете желаемый индекс для обоих m и n.

Большинство функций MATLAB принимают матрицы и векторы и вычисляют их результаты по элементам. Таким образом, чтобы получить двойную сумму, вы вычисляете все элементы суммы параллельно на f(m_, n_) и суммируете их.

В вашем случае (обратите внимание, что оператор .* выполняет поэлементное умножение матриц)

N = 10; 
M = 10; 
n = 1:N; 
m = 1:M; 

rho = 1; 
phi = 1; 
z = 1; 

% N rows x M columns for each matrix 
% n_ - all columns are identical 
% m_ - all rows are identical 
n_ = repmat(n', 1, M); 
m_ = repmat(m , N, 1); 

element_nm = cos (n_*z) .* besselj(m_-1, n_*rho) .* cos(m_*phi); 
sum_all = sum(element_nm(:));