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

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

I=1;u0=4*pi*1e-007; d=0.15;h=11.4;  a1=7.5;a2=14

问题
首先,说一下条件给的不清楚或可能有错的地方：
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、对于这个问题，我是花费了不少时间认真做了研究的，如果对你有帮助，希望能够采纳。