0

我想知道是否有任何方法可以绕过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
4

1 回答 1

2

我不确定为什么要避免使用该quad2d函数，但可以将 2D 正交拆分为两个 1D 嵌套正交（即使该函数未分解）。

这是一种处理与两个trapz调用的集成的方法。关键是在一个大表中生成所有值的集合以及trapz两个维度的使用。

首先，我定义测试值：

f= @(x,y) sin(x.*cos(y));
N = 1000;

然后，我为 x 和 y 方向定义网格（类似于您的k2and k3）：

xpts1d = linspace(0,1,N+1);
ypts1d = linspace(0,1,N+1);

然后我评估所有点对中的函数 f：

xpts = xpts1d'*ones(1,N+1);
ypts = ones(N+1,1)*ypts1d;
values = f(xpts,ypts);

然后，可以使用两个嵌套调用来完成集成：

I = trapz(xpts1d,trapz(ypts1d,values,2),1);

这给出了正确的答案。

在我的答案的第一个版本中，我使用了该quad函数，它直接与函数句柄一起使用。“天真的”想法是简单地嵌套两个调用，例如：

I = quad( @(y) quad( @(x)f(x,y) ,0,1),0,1)

但这实际上不起作用。关键是内部quad调用不适用于向量值函数，而外部调用需要向量值函数。

为了让它工作，我们只需要改变对 , 的向量值版本的内部quad调用quadv：

I = quad( @(y) quadv( @(x)f(x,y) ,0,1),0,1)

现在它应该可以工作了，我们可以检查它是否确实给出了相同的答案

I = quad2d( f, 0,1,0,1)

如果你真的想使用这个trapz函数，我建议你在上面构建一个函数trapz，接受一个函数句柄并trapz在里面调用。

于 2013-05-14T13:45:58.440 回答