Speeding up integration within nested for loops

25 ビュー (過去 30 日間)
Aravindh Shankar 2022 年 6 月 20 日

I have to find the elements of the following matrix, each element of which is
where K itself involves an integral
where f and q are known (I've left them out so as not to overload the question with math); they involve trigonometric functions of theta and square roots of theta,x,x'.
To find J I tried two approaches:
1) Numerical integration: I defined K(x,x') as K = @(x,x') integral(@theta ..Integrand(theta,x,x').., 0, pi)
and J = @(x,limit1,limit2) integral(@(x') K(x,x'), limit1, limit2) and evaluate them as
for i = 1:numel(x)
for j = 1:numel(x)
J_arr(i,j) = (1/(x(j+1) - x(j)))*J(x(i),x_arr(j),x_arr(j+1));
end
end
This works but the issue is that it is too slow for my purpose - the array x has 2047 elements so this script is carrying out both integrals 2047*2047 times. Is there a faster way to do this integration numerically? (Have thought along the lines of integral2, arrayfun so far).
2) Symbolic integration: Define K and J as above, but try to symbolically integrate them before the loops so that inside the loop I would just be substituting values for x_i and the limits. The issue with this is that the functions f and q in the definition of K are sufficiently complicated that Matlab has trouble finding the analytic symbolic integral, even for the first inner (theta) integration.
Any help in speeding up the numerical calculation/implementing a symbolic integration would be much appreciated!
Thanks

サインインしてコメントする。

採用された回答

Johan 2022 年 6 月 20 日
One trick I know to increase speed of calculation is to vectorize your code. I'm not sure that would work for you but you could try using the ArrayValued option of the integral function. Also it looks like the constant factor in Jij can be computed as an array once and multiplied with the result of the integral afterward.
x = rand(10,1);
Kint = @(x1,x2) integral(@(theta) sin(theta+x2+x1), 0, pi,'ArrayValued',true); %I put in a filler function here
Jint = @(x1,limit1,limit2) integral(@(x2) Kint(x1,x2), limit1, limit2,'ArrayValued',true);
prefactor = ones(numel(x),1).*(1./(x(2:end) - x(1:end-1)))'; %Compute prefactor for all ij
J_arr = zeros(numel(x),numel(x)-1); %initialize result array
for iJ = 1:numel(x)-1
J_arr(:,iJ) = Jint(x, x(iJ), x(iJ+1));
end
J_arr.*prefactor
ans = 10×9
1.5154 1.2476 1.0782 1.3314 1.3181 0.7666 1.1930 1.3678 1.4834 0.9472 0.5979 0.3974 0.7129 0.6969 0.0473 0.5474 0.7552 0.9151 0.9196 0.5679 0.3666 0.6839 0.6678 0.0159 0.5177 0.7263 0.8877 0.5684 0.1938 -0.0125 0.3196 0.3030 -0.3635 0.1486 0.3628 0.5382 1.4513 1.1705 0.9956 1.2592 1.2454 0.6764 1.1160 1.2966 1.4191 0.5364 0.1604 -0.0459 0.2869 0.2703 -0.3964 0.1158 0.3300 0.5064 0.2857 -0.0973 -0.3023 0.0330 0.0164 -0.6450 -0.1375 0.0756 0.2576 1.4131 1.1251 0.9472 1.2165 1.2025 0.6240 1.0707 1.2545 1.3808 0.7013 0.3336 0.1283 0.4563 0.4399 -0.2241 0.2864 0.4994 0.6703 1.6477 1.4118 1.2564 1.4840 1.4719 0.9648 1.3574 1.5176 1.6167
Hope this helps,
5 件のコメント3 件の古いコメントを表示3 件の古いコメントを非表示
Johan 2022 年 7 月 5 日
The integral function seems to have trouble reaching its default tolerance value for your function, I reduced the tolerance to 1e-4 I don't know if that is satisfactory for you or not.
With larger tolerance, your code is still faster for nested than for arrayed for small input size, however, I see the nested way seems to increase quadraticaly with input size whereas the arrayed increases somewhat linearly:
A small optimisation I did was un-nest your functions in the Kint integration. I noticed Matlab tends to slow down then you have function calling function I'm not sure why. Its much less readable and you should be extra careful the equation is right but that leads to faster computation.
step = [0.451,0.251,0.151,0.1];
timed = zeros(2,2);
N = zeros(2,1);
for i_step = 1:length(step)
X = -0.999:step(i_step):1;
%evaluate only once because of computation time limit
tic
nested(X);
timed(i_step, 1) = toc;
tic
arrayed(X);
timed(i_step, 2) = toc;
% Version I used on my computer for the graph above
% funnest = @() nested(X);
% funarr = @() arrayed(X);
% timed(i_step, 1) = timeit(funnest);
% timed(i_step, 2) = timeit(funarr);
N(i_step) = size(X,2);
end
plot(N,timed,'-*')
legend({'nested','arrayed'})
xlabel('X array size')
ylabel('Computing time (s)')
axis square
function nested(x)
%Constants
mstar = 1;
kappa = 1;
g_v = 1;
el = 1.38*1e4;
r_s = 5;
n = mstar^2*el^4/(pi*kappa^2*r_s^2);
k_F = sqrt(2*pi*n);
E_F = (k_F^2)/(2*mstar);
q_TF = 2*g_v*el^2*mstar/kappa;
Kint = @(x_1,x_2) (integral(@(theta) ((mstar/(2*pi^2)).*...
(2*pi*el^2./(kappa.*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - (4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta)))))).*...
(1 - ((sqrt(2*pi*el^2*n*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) -...
(4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/(kappa*mstar))).^2 ...
./(((sqrt(2*pi*el^2*n*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) -...
(4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/(kappa*mstar))).*...
sqrt(1 + ((sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - (4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/q_TF))).*...
(((sqrt(2*pi*el^2*n*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - (4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*...
cos(theta))))/(kappa*mstar))).*sqrt(1 + ((sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - ...
(4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/q_TF))) + E_F*abs(x_1) + E_F.*abs(x_2))))))...
,0,pi,'ArrayValued',true,'RelTol',0,'AbsTol',1e-4));
Jint = @(x_1,lim1,lim2) (integral(@(x_2) Kint(x_1,x_2),lim1,lim2,'ArrayValued',true,'RelTol',0,'AbsTol',1e-4));
J_arr = zeros(numel(x)-1);
for i = 1:length(x)-1
x_i_bar = (x(i) + x(i+1))/2;
for j = 1:numel(x)-1
J_arr(i,j) = (1./(x(j+1) - x(j))).*(Jint(x_i_bar,x(j),x(j+1)));
end
end
end
function arrayed(x)
%Constants
mstar = 1;
kappa = 1;
g_v = 1;
el = 1.38*1e4;
r_s = 5;
n = mstar^2*el^4/(pi*kappa^2*r_s^2);
k_F = sqrt(2*pi*n);
E_F = (k_F^2)/(2*mstar);
q_TF = 2*g_v*el^2*mstar/kappa;
Kint = @(x_1,x_2) (integral(@(theta) ((mstar/(2*pi^2)).*...
(2*pi*el^2./(kappa.*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - (4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta)))))).*...
(1 - ((sqrt(2*pi*el^2*n*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) -...
(4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/(kappa*mstar))).^2 ...
./(((sqrt(2*pi*el^2*n*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) -...
(4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/(kappa*mstar))).*...
sqrt(1 + ((sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - (4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/q_TF))).*...
(((sqrt(2*pi*el^2*n*(sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - (4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*...
cos(theta))))/(kappa*mstar))).*sqrt(1 + ((sqrt((2*mstar*E_F*(x_1 + x_2 + 2)) - ...
(4*mstar*E_F*sqrt(x_1 + 1)*sqrt(x_2 + 1).*cos(theta))))/q_TF))) + E_F*abs(x_1) + E_F.*abs(x_2))))))...
,0,pi,'ArrayValued',true,'RelTol',0,'AbsTol',1e-4));
Jint = @(x_1,lim1,lim2) (integral(@(x_2) Kint(x_1,x_2),lim1,lim2,'ArrayValued',true,'RelTol',0,'AbsTol',1e-4));
J_arr = zeros(numel(x)-1);
x_i_bar = (x(1:end-1) + x(2:end))/2;
for j = 1:numel(x)-1
J_arr(:,j) = (1./(x(2:end)-x(1:end-1))).*Jint(x_i_bar,x(j),x(j+1));
end
end
Aravindh Shankar 2022 年 7 月 5 日

As a future step I will try parallelizing using parfor once I have server access to see if that improves speed as well.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeSpecial Values についてさらに検索

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by