Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Replies: 0

 Francesco Perrone Posts: 39 Registered: 5/2/12
Posted: May 6, 2013 6:36 PM

I'd like knowing if there is any way to bypass the call for quad2d with a nested trapz loop. I'll discuss my problem with some more detail: currently, I perform the calculation of a double integral this way:

clc, clear all, close all

c = 1.476;
gamma = 3.0;

beta_int = (c*gamma)./(k_int.*sqrt(E_integral));

figure, loglog(k_int,beta_int,'r','LineWidth',2.0), grid on;

k1 = (.01:.1:100);
k2 = .01:.1:100;
k3 = -100:.1:100;

int_k3 = zeros(size(k2));
int_k3k2 = zeros(size(k1));

tic
for ii = 1:numel(k1)
phi11 = @(k2,k3) PHI11(k1(ii),k2,k3,k_int,beta_int);
end
toc

where PHI11 is defined as

function phi11 = PHI11(k1,k2,k3,k_int,beta_int)
k = sqrt(k1.^2 + k2.^2 + k3.^2);
ksq = k.^2;
k1sq = k1.^2;
fourpi = 4.*pi;
beta = exp(interp1(log(k_int),log(beta_int),log(k),'linear'));
k30 = k3 + beta.*k1;
k0 = sqrt(k1.^2 + k2.^2 + k30.^2);
k0sq = k0.^2;
k04sq = k0.^4;
Ek0 = (1.453.*k04sq)./((1 + k0sq).^(17/6));

C1 = (beta.*k1sq.*(k0sq - 2.*(k30.^2) + beta.*k1.*k30))./(ksq.*(k1.^2 + k2.^2));
C2 = ((k2.*k0sq)./((k1.^2 + k2.^2).^(3/2))).*atan2((beta.*k1.*sqrt(k1.^2 + k2.^2)),(k0sq - k30.*k1.*beta));
xhsi1 = C1 - (k2./k1).*C2;
xhsi1_sq = xhsi1.^2;
phi11 = (Ek0./(fourpi.*k04sq)).*(k0sq - k1sq - 2.*k1.*k30.*xhsi1 + (k1.^2 + k2.^2).*xhsi1_sq);
end

and E_integral.mat can be obtained this way:

clc,clear all,close all

k_int = .001:.01:1000;

Ek = (1.453.*k_int.^4)./((1 + k_int.^2).^(17/6));

E = @(k_int) (1.453.*k_int.^4)./((1 + k_int.^2).^(17/6));

E_integral = zeros(size(k_int));

for ii = 1:numel(k_int)
E_integral(ii) = integral(E,k_int(ii),Inf);
end

save('E_integral','k_int','E_integral')

Now the question is: would it be possible to overlook quad2d and the handle function in favor of a more practicle approach, by using a nested trapz function?

So far, I've tried the following piece of code, which has not yielded to the expected results:

int_k33 = zeros(size(k2));
S_11 = zeros(size(k1));
tic
for ii = 1:1
for jj = 1:numel(k2)
int_k33(jj) = trapz(k3,PHI11(k1(ii),k2(jj),k3,k_int,beta_int));
end
S_11(ii) = 4*trapz(k2,int_k33);
end
toc