求用matlab计算一个六重积分的具体程序,绘制仿真图

问题首先,说一下条件给的不清楚或可能有错的地方：1、按照公式中使用的符号,下面的常数d、h应为D、H,a1、a2应为R1、R2；2、公式中,第三重积分的下限确定是2*sqrt(R1^2-x)?【后面有讨论】3、求积分得到的结果就是一个具体的数,标题要求的【绘制仿真图】是什么意思? 参考代码下面给出我写的代码,供参考。% 常数定义 I=1; u0=4*pi*1e-007; D=0。15; H=11。4; R1=7。5; R2=14;   % 六元函数定义（注意把里面第二重积分的系数以及D^3折算到函数中） f=@(t,h,r,y,x,z)(-r*sin(t)。*(y-sin(t))-r*cos(t)。*(x-cos(t))) 。/ 。。。     ((x-r*cos(t))。^2+(y-r*sin(t))。^2+(z-h)。^2)。^(3/2) * u0*I^2/(4*pi)/D^3;   % 积分参数设置 tol = 1e-6; quadf = @quadl;trace = 1;  % 先使用triplequad计算内三重积分（积分限均为常数,h、r注意换算） Ithr = @(y,x,z) triplequad(f,0,2*pi,0,D*H,D*R1,D*R2,tol,quadf,y,x,z);   % 由于被积分函数必须能够接受向量输入并返回向量输出,所以外面的三重积分必须逐层进行Iy = @(x,z) quadf( @(Y)arrayfun(@(y)Ithr(y,x,z),Y),2*sqrt(R1^2-x),sqrt(R2^2-x。^2),tol,trace);Ix = @(z) quadf( @(X)arrayfun(@(x)Iy(x,z),X),R1,R2,tol,trace);Iz = quadf( @(Z)arrayfun(@(z)Ix(z),Z),0,D*H,tol,trace);注解及讨论1、对代码自身不想再多做解释,该说的注释里都说了。自己好好看看,实在看不懂再问吧（要更好地理解代码,需要好好熟悉一下匿名函数）。2、需要注意：重积分的计算时间不是根据积分的重数成倍增加,而是按照幂级数规律增加的,所以,计算六重积分所需的时间是比较长的（这段代码在我电脑上求解时间近8分钟）。3、按照当前的程序,求出的结果：Iz =  1。1554e-005由于y的积分下限看起来有点怪,像是sqrt(R1^2-x。^2)之误,但仔细分析应该没错。如果积分下限取sqrt(R1^2-x。^2),则计算结果为复数：Iz =  5。4730e-005 -6。7361e-005i原因是,对x的积分范围是R1-R2,也就意味着,x的取值范围最大可以是R2,而里面那一重积分的下限是sqrt(R1^2-x。^2),当x取值大于R1时就会出现复数。4、上面的代码侧重于实现方法,具体的数学表达式不保证绝对没问题,请自行仔细核实。 再问： 膜拜了。。。上面那个地方是我写错了，应该是sqrt(R1^2-x。^2)，问题问的确实也不大清楚，因为最终想得到的是Y（本式子的结果）与其中一些变量的关系（如R1、R2、H），所以想用图来表示，不知道可不可以实现 再答： 1、上面讨论过，y的积分下限如果是sqrt(R1^2-x。^2)，计算结果会出现复数，你确定是你想要的结果？如果不是，请认真核实公式（需要知道公式对应的实际问题是什么），然后再计算，以避免花费大量时间计算得到无意义的结果。 2、如果想改变R1、R2、H之类的变量计算相应的Y，把上面除第一行常数定义之外的代码用循环框起来计算即可，例如（代码只是说明大致的思路，未实际测试）：% 常数定义I=1; u0=4*pi*1e-007; 。。。。。。HH = 11:0。1:12;Y=zeros( size(HH) );for ii=1:length(HH)    H=HH(ii);    % 调用积分相关代码    。。。。。。    Y(ii)=Iz;end如果想同时改变两个变量，用两重循环即可。需要特别强调的是，单次计算六重积分的耗时已经很可观，而循环意味着多次计算，所以对需要的时间要有心理准备。 顺便说一下，上面程序中之所以设置trace=1，就是因为计算时间太长，看起来像是死机了一样，所以让迭代过程中输出一些信息，表明计算还在继续。如果不希望输出这些信息，设置trace=0即可（三个quadf可以分别设置）。 3、对于这个问题，我是花费了不少时间认真做了研究的，如果对你有帮助，希望能够采纳。