通俗易懂的monteCarlo积分方法(八)

计算任意重积分

一.基本问题
二.问题的分析
三.matlab实现代码
- 1.基本的求解函数intCal.m
- 2.多重积分计算的GUI实现
- - (1)基本架构
  - (2)结果显示
四.结论拓展
五.直接调用代码

用matlab设计计算一般多重积分的算法和GUI界面

一.基本问题

先考虑计算多重积分的问题:
∫Ωf(x→)dΩ\int_{\Omega}f(\overrightarrow x)d \Omega ∫Ωf(x)dΩ
其中：x→=(x1,x2,...xn)T\overrightarrow x = (x_1,x_2,...x_n)^Tx=(x1,x2,...xn)T ，且Ω:{(x1,x2,...xn)∣A(x→)<0}\Omega:\{(x_1,x_2,...x_n)|A(\overrightarrow x)<0\}Ω:{(x1,x2,...xn)∣A(x)<0}，其中AAA是(x1,x2,...xn)(x_1,x_2,...x_n)(x1,x2,...xn)总的约束方程，为p×1p\times 1p×1维的向量。且对每个xi(i=1,2,3...n)x_i(i = 1,2,3...n)xi(i=1,2,3...n)都有下界约束：a→:[a1,a2,.ai,..an]\overrightarrow a:[a_1,a_2,.a_i,..a_n]a:[a1,a2,.ai,..an],同时都是上界约束:

b→:[b1,b2,.bi,..bn]\overrightarrow b:[b_1,b_2,.b_i,..b_n]b:[b1,b2,.bi,..bn]。在这里我们求这个多重积分的积分值计算函数(用Matlab求解)同时绘制GUI界面，同时估计其精度。

二.问题的分析

Input:Input:Input:f(x1,x2,..xn)f(x_1,x_2,..x_n)f(x1,x2,..xn),a:(a1,a2,...an)Ta:(a_1,a_2,...a_n)^Ta:(a1,a2,...an)T,b:(b1,b2,...bn)Tb:(b_1,b_2,...b_n)^Tb:(b1,b2,...bn)T,A:(A1(x1,x2,...xn),A2(x1,x2...xn),...Ap(x1,x2,..xn))TA:(A_1(x_1,x_2,...x_n),A_2(x_1,x_2...x_n),...A_p(x_1,x_2,..x_n))^TA:(A1(x1,x2,...xn),A2(x1,x2...xn),...Ap(x1,x2,..xn))T

Output:IntValue,errorOutput:IntValue,errorOutput:IntValue,error

用MonteCarloMonteCarloMonteCarlo方法计算多重积分的值：

考虑产生在[a→,b→][\overrightarrow a,\overrightarrow b][a,b]内的均匀随机变量ζ→\overrightarrow \zetaζ，其密度函数为g(x→)=1Ωg(\overrightarrow x) = \frac{1}{\Omega}g(x)=Ω1,其中Ω=∏i=1n(bi−ai)\Omega = \prod_{i=1}^n(b_i-a_i)Ω=∏i=1n(bi−ai)。

那么其数学期望可以表示为:
Ef(x→)g(x→)=∫Ωf(x→)g(x→)g(x→)dΩ=∫Ωf(x→)dΩE\frac{f(\overrightarrow x)}{g(\overrightarrow x)} =\int_{\Omega}\frac{f(\overrightarrow x)}{g(\overrightarrow x)}g(\overrightarrow x)d \Omega = \int_{\Omega}f(\overrightarrow x)d \Omega Eg(x)f(x)=∫Ωg(x)f(x)g(x)dΩ=∫Ωf(x)dΩ
此时由大数定律：
∀ϵ>0,limN→∞P{∣1N∑i=1Nf(ζ→i)g(ζ→i)−∫Ωf(x→)dΩ∣<ϵ}=1\forall \epsilon>0,lim_{N\rightarrow \infty}P\{|\frac{1}{N}\sum_{i=1}^N\frac{f(\overrightarrow \zeta_i)}{g(\overrightarrow \zeta_i)}-\int_{\Omega}f(\overrightarrow x)d \Omega|<\epsilon\}=1 ∀ϵ>0,limN→∞P{∣N1i=1∑Ng(ζi)f(ζi)−∫Ωf(x)dΩ∣<ϵ}=1
上式构成了我们计算多重积分的基础。再由中心极限定理：
∀ϵ>0,limN→∞P{∣1N∑i=1Nf(ζ→i)g(ζ→i)−∫Ωf(x→)dΩ∣<λασN}=22π∫0λαe−t22dt=1−α\forall \epsilon>0,lim_{N\rightarrow \infty}P\{|\frac{1}{N}\sum_{i=1}^N\frac{f(\overrightarrow \zeta_i)}{g(\overrightarrow \zeta_i)}-\int_{\Omega}f(\overrightarrow x)d \Omega|<\frac{\lambda_{\alpha}\sigma}{\sqrt{N}}\}=\frac{2}{\sqrt{2\pi}}\int_0^{\lambda_{\alpha}}e^{-\frac{t^2}{2}}dt = 1-\alpha ∀ϵ>0,limN→∞P{∣N1i=1∑Ng(ζi)f(ζi)−∫Ωf(x)dΩ∣<Nλασ}=2π2∫0λαe−2t2dt=1−α
上面的式子构成了我们依概率估计精度的基础。

三.matlab实现代码

1.基本的求解函数intCal.m

function [intVal,error,xPoint,yPoint] = intCal(a,b,p,f,A,numbers)
%UNTITLED 此处显示有关此函数的摘要
%   有p个约束条件
[m1 n1] = size(a);
[m2 n2] = size(b);
if (m1 ~= m2)&&(n1~=n2)disp('Error of dimension of a and b!');return;
end
sum1 = 0;
sum2 = 0;
zeta = zeros(m1,1);
Omega = prod(b-a);
for k = 1:numbersfor i = 1:m1zeta(i) = unifrnd(a(i),b(i));endif sum(A(zeta)<0) == psum1 = sum1 + f(zeta);sum2 = sum2 + f(zeta)^2;endxPoint(k) = k;yPoint(k) = Omega*sum1/k ;
end
intVal = Omega*sum1/numbers ;
for k = 1:numberszPoint(k) = intVal;
end
sigma = sqrt(sum2/numbers-(sum1/numbers)^2);
error = 3*sigma/sqrt(numbers);%以99%的概率估计误差
hold on
plot(xPoint,yPoint,'b');
plot(xPoint,zPoint,'r');
legend('迭代值','实际值');
xlabel('迭代次数');
ylabel('积分值');
title('积分迭代曲线');
hold off
end

我们求以下三重积分的值做一下实验：
∭Vx2dV(V:{(x,y,z)∣(x2+y2<z<1)})\iiint_{V}x^2dV(V:\{(x,y,z)|(\sqrt{x^2+y^2}<z<1)\}) ∭Vx2dV(V:{(x,y,z)∣(x2+y2<z<1)})

%Input:
>> a = [-1;-1;0];
>> b = [1;1;1];
>> p = 2;
>> f = @(x)(x(1)^2);
>> numbers = 1000;
>> A =@(x)[x(3)-1;sqrt(x(1)^2+x(2)^2)-x(3)];
%Output:
>> intCal(a,b,p,f,A,numbers)ans =0.1545

该积分值与实际值π20\frac{\pi}{20}20π相差无几。

2.多重积分计算的GUI实现

(1)基本架构

(2)结果显示

四.结论拓展

已知了重积分的计算公式之后，如果要计算曲面积分的公式也可用广义的牛顿斯托克斯公式进行转化后计算：
∫∂Ωω=∫Ωdω\int_{\partial \Omega}\omega = \int_{\Omega}d\omega ∫∂Ωω=∫Ωdω

五.直接调用代码

clc,clear;
aCell =  inputdlg('请输入各个变量的下界限(列向量)','变量下界',1,{'-1;-1 ;0'});
bCell =  inputdlg('请输入各个变量的上界限(列向量)','变量上界',1,{'1;1 ;1'});
a = str2num(['[',aCell{1},']']);
b = str2num(['[',bCell{1},']']);
m1 = size(b,1);
m2 = size(a,1);
if (m1 ~= m2)&&(n1~=n2)disp('Error of dimension of a and b!');return;
end
%以下是计算多重积分的值
numbers =  inputdlg('请输入计算次数','计算次数',1,{'5000'});
numbers = str2num(numbers{1});
f = inputdlg('请输入被积函数','被积函数',1,{'x(1)^2'});
f = str2func(['@(x)(',f{1},')']);
A = inputdlg('请输入基本约束条件(A<=0)','约束条件',1,{'x(3)-1;sqrt(x(1)^2+x(2)^2)-x(3)'});
p = size(strfind(A{1},';'),2) + 1;
endA = A{1};
endA = endA(end);
if endA == ';'p = p - 1;
end
A = str2func(['@(x)[',A{1},']']);
sum1 = 0;
sum2 = 0;
zeta = zeros(m1,1);
Omega = prod(b - a);
for k = 1:numbersfor i = 1:m1zeta(i) = unifrnd(a(i),b(i));endif sum(A(zeta)<0) == psum1 = sum1 + f(zeta);sum2 = sum2 + f(zeta)^2;endxPoint(k) = k;yPoint(k) = Omega*sum1/k ;sum1Value(k) = sum1;sum2Value(k) = sum2;sigmaValue(k) = sqrt(sum2Value(k)/k -(sum1Value(k)/k)^2);errorValue(k) = 3*sigmaValue(k)/sqrt(k);maxPoint(k) = yPoint(k) +  errorValue(k);minPoint(k) = yPoint(k) -  errorValue(k);
end
intVal = Omega*sum1/numbers ;
for k = 1:numberszPoint(k) = intVal;
end
sigma = sqrt(sum2/numbers-(sum1/numbers)^2);
error = 3*sigma/sqrt(numbers);
step = inputdlg('请设置置信区间显示步长','显示步长',1,{'50'});
step = str2num(step{1});
hold on
h = figure(1);
plot(xPoint,yPoint,'b');
plot(xPoint,zPoint,'r');
plot(xPoint(1:step:numbers),maxPoint(1:step:numbers),'o','color','y');
plot(xPoint(1:step:numbers),minPoint(1:step:numbers),'o','color','c');
legend('迭代值','实际值','最大置信值','最小置信值');
xlabel('迭代次数');
ylabel('积分值');
title(['积分值是:',num2str(intVal),' 有99.7%的概率认为置信区间是:',...'[',num2str(intVal - error),',',num2str(intVal + error),']']);
hold off