Clever Geek Handbook

Speed ā€‹ā€‹up nested loop

I am working on speeding up the following function, but with no results:

function beta = beta_c(k,c,gamma)
beta = zeros(size(k));
E = @(x) (1.453*x.^4)./((1 + x.^2).^(17/6));
for ii = 1:size(k,1)
    for jj = 1:size(k,2)
        E_int = integral(E,k(ii,jj),10000);
        beta(ii,jj) = c*gamma/(k(ii,jj)*sqrt(E_int));
    end
end
end

So far I've solved it like this:

function beta = beta_calc(k,c,gamma)
k_1d = reshape(k,[1,numel(k)]);
E_1d =@(k) 1.453.*k.^4./((1 + k.^2).^(17/6));
E_int = zeros(1,numel(k_1d));
parfor ii = 1:numel(k_1d)
E_int(ii) = quad(E_1d,k_1d(ii),10000);
end
beta_1d = c*gamma./(k_1d.*sqrt(E_int));
beta = reshape(beta_1d,[size(k,1),size(k,2)]);
end

I don't think it improved the performance. What do you think about this?

Could you shed some light?

I thank you in advance.

EDIT

I'm going to present some theoretical background related to my question. Beta is usually calculated as follows:

enter image description here

Therefore, in the given case of one-dimensional k array E_int can be calculated as

E = 1.453.*k.^4./((1 + k.^2).^(17/6));
E_int = 1.5 - cumtrapz(k,E);

or alternatively like

E_int(1) = 1.5;
for jj = 2:numel(k)
E =@(k) 1.453.*k.^4./((1 + k.^2).^(17/6));
E_int(jj) = E_int(jj - 1) - integral(E,k(jj-1),k(jj));
end

However, it k

is currently a matrix k(size1,size2)

.

+3


source to share


2 answers


Here's a different approach, parallel because it's easy to use with spmd

or parfor

. Instead of integral

consider quad

, see link for examples ...



+2


source


I like this question.

Problem: the function integral

accepts only scalars as integration constraints. Hence, it is difficult to vectorize the computation E_int

.



Key: there seems to be a lot of redundancy in integrating the same function over and over from k(ii,jj)

ad infinitum ...

Suggested solution: how about sorting the values k

from smallest to largest and integrating E_sort_int(si) = integral( E, sortedK(si), sortedK(si+1) );

with sortedK( numel(k) + 1 ) = 10000;

. Then the full value E_int = cumsum( E_sort_int );

(you only need to "undo" the sort and resize it to k

).

+2


source






More articles:


All Articles