任意重积分的计算直接调用代码(Matlab)

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:numbers
    for i = 1:m1
        zeta(i) = unifrnd(a(i),b(i));
    end
    if sum(A(zeta)<0) == p
        sum1 = sum1 + f(zeta);
        sum2 = sum2 + f(zeta)^2;
    end
    xPoint(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:numbers
    zPoint(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