Speeding up integration within nested for loops
10 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Aravindh Shankar
il 20 Giu 2022
Modificato: Aravindh Shankar
il 5 Lug 2022
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
0 Commenti
Risposta accettata
Johan
il 20 Giu 2022
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
Hope this helps,
5 Commenti
Johan
il 5 Lug 2022
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
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Applications in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!