Matlab多重积分的两种实现【从六重积分到一百重积分】

问题

今天被问了一个问题：

$$\begin{array}{l}\mu =\int \int \int \int \int \int f(x_1,x_2,x_3,x_4,x_5,x_6)d{x}_{1}d{x}_{2}d{x}_{3}d{x}_{4}d{x}_{5}d{x}_6\\ {\sigma }^2=\int \int \int \int \int \int {\left[f(x_1,x_2,x_3,x_4,x_5,x_6)-\mu \right]}^{2}d{x}_{1}d{x}_{2}d{x}_{3}d{x}_{4}d{x}_{5}d{x}_6\end{array}$$
约束：这个函数f的返回值是一个向量。

要求：用Matlab解决。

方案一

function ret = integral6(fun, lb, ub)
% int x2
fx1 = @(x2, x3, x4, x5, x6)integral(@(x1)fun(x1, x2, x3, x4, x5, x6), lb(1), ub(1), 'ArrayValued', 1);

% int x2
dfx2 = @(x2, x3, x4, x5, x6)cell2mat(arrayfun(@(xx2)fx1(xx2, x3, x4, x5, x6), x2, 'UniformOutput', 0));
fx2 = @(x3, x4, x5, x6)integral(@(x)dfx2(x, x3, x4, x5, x6), lb(2), ub(2), 'ArrayValued', 1);

% int x3
dfx3 = @(x3, x4, x5, x6)cell2mat(arrayfun(@(xx3)fx2(xx3, x4, x5, x6), x3, 'UniformOutput', 0));
fx3 = @(x4, x5, x6)integral(@(x)dfx3(x, x4, x5, x6), lb(3), ub(3), 'ArrayValued', 1);

% int x4
dfx4 = @(x4, x5, x6)cell2mat(arrayfun(@(xx4)fx3(xx4, x5, x6), x4, 'UniformOutput', 0));
fx4 = @(x5, x6)integral(@(x)dfx4(x, x5, x6), lb(4), ub(4), 'ArrayValued', 1);

% int x5
dfx5 = @(x5, x6)cell2mat(arrayfun(@(xx5)fx4(xx5, x6), x5, 'UniformOutput', 0));
fx5 = @(x6)integral(@(x)dfx5(x, x6), lb(5), ub(5), 'ArrayValued', 1);
% int x6
dfx6 = @(x6)cell2mat(arrayfun(@(xx6)fx5(xx6), x6, 'UniformOutput', 0));
fx6 = integral(@(x)dfx6(x), lb(6), ub(6), 'ArrayValued', 1);

ret = fx6;
end

看着很玄乎……其实都没啥，一层一层剥洋葱。这个解法唯一的意义在于学习几个函数的用法：arrayfun, cell2mat, integral，以及@定义函数的语法。因为计算时间幂级的，根本没用。

能够在有限时间得到结果的是大名鼎鼎的Monte Carlo法。

Monte Carlo法

关于这个方法的数学，这里不赘述，需要知道的，搜索了sci-hub；那些不知道怎么搜索的，大概率也不需要知道。

function ret = integral6mc(fun, lb, ub, N)
% fun是被积分的函数
% lb和ub是积分变量的范围，每个都是六维数组
% N MC采样的数目，一般取个几千、几万试一下差不多就行
V = prod(abs(ub-lb)); % 计算超立方体的体积，向量化计算每一个维度的长度（绝对值），求所有维度的乘积
n = length(lb); % 维数
sample = (ub-lb) .* rand(N, n) + repmat(lb,N,1); %产生一个Nxn的随机数组，然后转换到lb~ub的范围，repmat函数是复制矩阵，把行向量复制程Nxn

sample_arguments = num2cell(sample, 1);% 把上面的Nxn矩阵换成按列排列的cell（n个元素）
results = cell2mat(arrayfun(fun, sample_arguments{:}, 'UniformOutput', 0));%调用被积函数，被积函数的参数有n个，把cell展开（{:}操作），这里arrayfun得到的是cell，再合并成mat，就是N个结果的向量
ret = sum(results) * V / N; % 这是MC的核心算法，乘以体积除以样本数

这个方法就快多了……

而且，这个函数可以自动适应积分重数，被积函数可以是向量也可以是标量，只要fun的参数个数，lb、ub的大小能够适配就行。

是男人就积分100重啊！！

蒙特卡洛求积分之第二篇