嵌套数值积分答案

【问题标题】：Nested numerical integration嵌套数值积分
【发布时间】：2013-11-10 20:24:49
【问题描述】：

链接中的问题：可以解析积分，答案是 4，但是我有兴趣使用 Matlab 对它进行数值积分，因为它在形式上类似于我无法解析积分的问题。因为两个内积分中的函数是x、y和z和z的函数，所以出现了数值积分的困难。

【问题讨论】：

你好。 :-) 不相信我，是吗？也许有人想出了一个解决方案。一条评论：没有什么取决于 x，所以你可以去掉一个内积分。 \int_0^1 dx = 1。但不会改变问题。
我删除了我的答案，这不完全是一个解决方案。更好地指出作为注释，使用符号数学工具箱，定义 syms y z 你会得到 4 对于 int(z * exp(int(1 / (y + z), y, 0, 1)), z, 0, 2)，否则，integral2 在这里不适用。
@A.Donda 是的，我的目标是这样写来概括问题。我希望这里有人提出一个优雅的解决方案。
我实际上想出了一个答案。它可能相当于 Guddu 的，只是使用 Matlab 的集成函数 quad 实现的，因此更优雅，可能在数值上更精确。既然我现在意识到我在你之前的问题上写的完全错误，我宁愿删除那个答案。但也许您想删除整个问题？

【解决方案1】：

嗯，这很奇怪，因为在发布者之前的类似问题上，我声称无法做到这一点，现在看了 Guddu 的回答后，我意识到它并没有那么复杂。我之前写的，数值积分结果是一个数字而不是一个函数，这是正确的——但题外话：人们可以定义一个函数来评估每个给定参数的积分，这种方式实际上是具有数值积分的结果。

不管怎样，就这样吧：

function q = outer

    f = @(z) (z .* exp(inner(z)));
    q = quad(f, eps, 2);

end

function qs = inner(zs)
% compute \int_0^1 1 / (y + z) dy for given z

    qs = nan(size(zs));
    for i = 1 : numel(zs)
        z = zs(i);
        f = @(y) (1 ./ (y + z));
        qs(i) = quad(f, 0 , 1);
    end

end

我在评论中应用了我自己建议的简化，消除了 x。函数inner 计算y 上的内积分值作为z 的函数。然后函数 external 计算 z 上的外积分。我通过让积分从 eps 而不是 0 运行来避免 z = 0 处的极点。结果是

4.00000013663955

inner 必须使用for 循环来实现，因为赋予quad 的函数需要能够同时为多个参数值返回其值。

【讨论】：

感谢您的回答
@Guddu 谢谢你的回答。

【解决方案2】：

这绝不是优雅的。希望有人能比我更好地使用matlab函数。我已经尝试过蛮力的方式只是为了练习数值积分。我试图通过利用它也乘以 z 的事实来避免 z=0 处的内部积分中的极点。我得到 3.9993。有人必须通过使用比梯形规则更好的东西来获得更好的解决方案

function []=sofn
clear all

global x y z xx yy zz dx dy

dx=0.05;
x=0:dx:1;
dy=0.002;
dz=0.002;
y=0:dy:1;
z=0:dz:2;

xx=length(x);
yy=length(y);
zz=length(z);

s1=0;
for i=1:zz-1
    s1=s1+0.5*dz*(z(i+1)*exp(inte1(z(i+1)))+z(i)*exp(inte1(z(i))));
end
s1

end

function s2=inte1(localz)
global y yy dy

if localz==0
    s2=0;
else
s2=0;
for j=1:yy-1
    s2=s2+0.5*dy*(inte2(y(j),localz)+inte2(y(j+1),localz));
end
end

end

function s3=inte2(localy,localz)
global x xx dx

s3=0;
for k=1:xx-1
    s3=s3+0.5*dx*(2/(localy+localz));
end

end

【讨论】：