求C语言求多重定积分的代码

变步长梯形求世野定积分代码如下：

#include

#include

#include

#define e 2.7182818

double f(double x)

{

return pow(e,x)*x*x//猜漏 函数x^2 * e^x

}

double Calc(double a,double b,double esp)//变步长搜兆喊梯形求[a,b]定积分,esp是精度

{

int done(0)

int n=1

double h,Tn,T2n,k,temp,x

h=b-a

Tn=h*(f(a)+f(b))/2.0

while(!done)

{

temp = 0

for(k=0k

问题

首先，说一下条件给的不清楚或可能有错的地方：

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、上面的代码侧重于实现方法，具体的数学表达式不保证绝对没问题，请自行仔细核实。

---