使用匿名函数进行多维集成时出错

Question

我正在尝试使用以下算法执行多维积分：

y= @(a,b,c) a+b+c; %function - 3 dim
y1=@(a,b) integral(@(c) y(a,b,c),-20,20); % L,H - limits for  c`
y2=@(a) integral(@(b) y1(a,b),-20,20); % L,H - limits for b
y3=integral(@(a) y2(a),-20,20); % L,H - limits for a

但这会产生以下错误：

Error using integralCalc/finalInputChecks (line 515)
Output of the function must be the same size as the input. 
If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
                finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in @(a)integral(@(b)y1(a,b),-20,20)


Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

有人可以帮助我理解和纠正我的错误，或者提出更好的方法吗？

附注
integral3函数我知道，但是我需要这个方法，因为我以后要尝试4,5,6....维度。

Answer 1

我不确定这是否适用于任何可能的情况，但对于这个简单的情况来说效果很好。只需使用符号数学。

syms a b c
y=a+b+c;
y1=int(y,c,-20,20)
y2=int(y1,b,-20,20)
y3=int(y2,a,-20,20)

但是，请小心创建变量。不要动态创建yn！

Answer 2

要了解为什么您会收到此错误，让我们使用“常规”函数重写您的代码：

function q42536274
  y3(-20,20);
end

function out = y(a,b,c)
  out = a + b + c;
  % The result of the above is some [1x150 double] vector. No problem here.
end

function out = y1(a,b,L,H)
  out = integral(@(c)y(a,b,c),L,H);
  % The result of the above is the scalar [-1.421085471520200e-14]. Problem!
end

function out = y2(a,L,H)
  out = integral(@(b)y1(a,b,L,H),L,H);
end

function out = y3(L,H)
  out = integral(@(a)y2(a,L,H),L,H);
end

这是发生错误时工作区的样子：

现在我们可以看到 MATLAB 抱怨的是什么：fx 和 x 中的元素数量不同！在这种情况下，MATLAB 应该如何进行数值积分？ 0-阶近似？这是模棱两可的。

我们需要告诉 MATLAB 如何摆脱这种混乱局面。我们可以这样做的一种方法是：

function out = y1(a,b,L,H)
  out = ones(size(a))*integral(@(c)y(a,b,c),L,H);
  % Now the result of the above is also a [1x150 double] vector. Yey!
end

function out = y2(a,L,H)
  out = ones(size(a))*integral(@(b)y1(a,b,L,H),L,H);
  % Same as above...
end

因此我们得到的输出-2.2737e-11 与0 的正确答案“非常接近”。