matlab - ネストされたtrapz二重統合

Question

quad2dネストされたtrapzループで呼び出しをバイパスする方法があるかどうか知りたいです。私の問題についてもう少し詳しく説明します。現在、この方法で二重積分の計算を実行しています。

clc, clear all, close all
load E_integral.mat

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);
    int_k3(ii) = 2*quad2d(phi11,-100,100,-100,100);
end
toc

は次のようにPHI11定義されます

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

次の方法でE_integral.mat取得できます。

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')

ここでの問題は、ネストされた関数を使用して、見落としquad2dてhandle functionより実用的なアプローチを支持することは可能でしょうか?trapz

これまでのところ、次のコードを試しましたが、期待どおりの結果が得られませんでした。

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