一般区域n重积分MATLAB计算方法[有代码]

coohit

这里说的n理论上讲不限大小，但是鉴于n重积分的特点，如果积分重数高，同时积分区域范围较大，数值积分计算量十分巨大，如果各重积分原函数都能显式表示，采用符号积分是第一选择。当只能数值积分的时候，很多时候4重5重积分计算时间都较长。所以虽说不限制n的大小，但是能实际应用中n一般不会很大。
    下面介绍下这个函数：f = nIntegrate(fun,low,up)
      f，函数返回值是n重积分积分结果。fun是被积函数字符串形式，不同的变量依次以x1,x2,...xn表示，（注意，必须以x1,x2,...xn这种形式，其余像y1,y2,...yn或是其他表示方法都不行），low和up都是长度为n的cell数组，low存储从外到内各重积分的积分下限函数，up存储从外到内各重积分的积分上限函数，（都是字符串形式）。举例说明：
计算如下积分：x1*x2*x3，x1从1到2，x2从x1到2*x1，x3从x1*x2到2*x1*x2，构造输入变量
fun = 'x1*x2*x3';
up = {'2','2*x1','2*x1*x2'};
low = {'1','x1','x1*x2'};
然后按f = nIntegrate(fun,low,up)调用即可。
需要说明的是这次没打算上传源m文件，而是上传的加密后的p文件，但是不影响函数使用。该函数可运行7.0以后的MATLAB版本，其中又以2009a为界设置了不同的求解方法。如果用的是2009a,那么速度明显会比09a以前的版本快不少。譬如下面附件里比较复杂的问题：在我的电脑（CPU 酷睿2 T8100 2G内存）上用2009a算得到结果要560多秒，而用08a要5000多秒。需要说明的是由于函数是调用现成的MATLAB函数来编写n重积分，所以函数运行效率也不是很高。

• fun5 = 'sin(x1*exp(x2*sqrt(x3)))+x4^x5'
• up5 = {'1','exp(x1)','x1+sin(x2)','x1+x3','2*x4'}
• low5 = {'0','exp(x1)/2','x1/2+sin(x2)/2','x1/2+x3/2','x4'}
• tic;f = nIntegrate(fun5,low5,up5),toc
• fun5 =
• sin(x1*exp(x2*sqrt(x3)))+x4^x5
• up5 =
•     '1'    'exp(x1)'    'x1+sin(x2)'    'x1+x3'    '2*x4'
• low5 =
•   Columns 1 through 4
•     '0'    'exp(x1)/2'    'x1/2+sin(x2)/2'    'x1/2+x3/2'
•   Column 5
•     'x4'
• f =
•     5.7981
• Elapsed time is 564.504780 seconds.
  

nIntegrate.p.rar (2.21 KB, 下载次数: 0)

2016-7-14 10:26 上传

点击文件名下载附件

X-man

Osmond

1＃的积分用Forcal求解：
高斯算法（通常耗时较长）：

• mvar:
• f(x1,x2,x3,x4,x5)=sin{x1*exp[x2*sqrt(x3)]}+x4^x5;
• yx(j,x1,x2,x3,x4,x5,y0,y1)=
• {
•     which{
•         j==0:[y0=0, y1=1],
•         j==1:[y0=exp(x1)/2, y1=exp(x1)],
•         j==2:[y0=[x1+sin(x2)]/2, y1=x1+sin(x2)],
•         j==3:[y0=(x1+x3)/2, y1=x1+x3],
•         j==4:[y0=x4, y1=2*x4]
•     }
• };
• t0=sys::clock();
• XSLSF::igaus[HFor("f"),HFor("yx"),10,10,10,10,10];  //一般，区间分隔数10,10,10,10,10越大越精确，但耗时越长
• [sys::clock()-t0]/1000;
• 
  

[color=rgb(51, 102, 153) !important]复制代码
结果：
5.798135338576838
耗时113.875秒

用连分式法混合了IMSL算法（精度取1e-6）：

• mvar:
• Fx4x5(x4,x5)=sin{x1*exp[x2*sqrt(x3)]}+x4^x5;
• x4x5(x4,y0,y1)=
• {
•     y0=x4,
•     y1=2*x4
• };
• Fx2x3(_x2,_x3)= x2=_x2, x3=_x3, XSLSF::pqg2[HFor("Fx4x5"),(x1+x3)/2,(x1+x3),HFor("x4x5"),1e-6];
• Gx2(x2)=[x1+sin(x2)]/2;
• Hx2(x2)=x1+sin(x2);
• Fx1(_x1)= x1=_x1, IMSL::TWODQ[HFor("Fx2x3"),exp(x1)/2,exp(x1),HFor("Gx2"),HFor("Hx2"),0,1e-6,2,0];
• t0=sys::clock();
• IMSL::QDAGS[HFor("Fx1"),0,1,0,1e-6,0];
• [sys::clock()-t0]/1000;
• 
  

[color=rgb(51, 102, 153) !important]复制代码
结果：
5.798135342841069

Osmond

8＃的积分用Forcal求解：

• mvar: a=1, omiga=2, c1=3, c2=4;
• Ft(tao)=exp(-tao*tao);
• F(t)=cos(omiga*t)*{c1*exp(-t*t)-c2*XSLSF::fpqg[HFor("Ft"),-a,t,1e-6]};
• IMSL::QDAGS[HFor("F"),-a,a,0,1e-6,0];
  

结果：-0.1417543457105889


完全用IMSL函数，且相对精度取1e-6：

• mvar: a=1, omiga=2, c1=3, c2=4;
• Ft(tao)=exp(-tao*tao);
• F(t)=cos(omiga*t)*{c1*exp(-t*t)-c2*IMSL::QDAGS[HFor("Ft"),-a,t,0,1e-6,0]};
• IMSL::QDAGS[HFor("F"),-a,a,0,1e-6,0];
  

结果：
-0.1417543562201364

Osmond

完全使用IMSL的Forcal解法
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
用单重积分构造5重积分（相对精度取1e-1）：

• mvar:
• Fx5(x5)=sin{x1*exp[x2*sqrt(x3)]}+x4^x5;
• Fx4(_x4)= x4=_x4, IMSL::QDAGS[HFor("Fx5"),x4,2*x4,0,1e-1,0];
• Fx3(_x3)= x3=_x3, IMSL::QDAGS[HFor("Fx4"),(x1+x3)/2,(x1+x3),0,1e-1,0];
• Fx2(_x2)= x2=_x2, IMSL::QDAGS[HFor("Fx3"),[x1+sin(x2)]/2,x1+sin(x2),0,1e-1,0];
• Fx1(_x1)= x1=_x1, IMSL::QDAGS[HFor("Fx2"),exp(x1)/2,exp(x1),0,1e-1,0];
• t0=sys::clock();
• IMSL::QDAGS[HFor("Fx1"),0,1,0,1e-1,0];
• [sys::clock()-t0]/1000;
  

结果：
5.798135465850306
时间2.953秒

如果相对精度取1e-2或更小，计算时间延长，但结果越来越不准确，这个还不知道为什么？看来用单重构造多重不是很好的方法。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
用单重＋二重＋二重积分构造5重积分（相对精度取1e-3）：

• mvar:
• Fx4x5(x4,x5)=sin{x1*exp[x2*sqrt(x3)]}+x4^x5;
• Gx4(x4)=x4;
• Hx4(x4)=2*x4;
• Fx2x3(_x2,_x3)= x2=_x2, x3=_x3, IMSL::TWODQ[HFor("Fx4x5"),(x1+x3)/2,(x1+x3),HFor("Gx4"),HFor("Hx4"),0,1e-3,2,0];
• Gx2(x2)=[x1+sin(x2)]/2;
• Hx2(x2)=x1+sin(x2);
• Fx1(_x1)= x1=_x1, IMSL::TWODQ[HFor("Fx2x3"),exp(x1)/2,exp(x1),HFor("Gx2"),HFor("Hx2"),0,1e-3,2,0];
• t0=sys::clock();
• IMSL::QDAGS[HFor("Fx1"),0,1,0,1e-3,0];
• [sys::clock()-t0]/1000;
  

结果：
5.798135365426686
耗时6.141秒。

相对精度取1e-2或1e-3，计算结果是一样的，时间也几乎相等；但相对精度取1e-4或更小时，计算时间延长，且结果越来越不准确，这个还不知道为什么？
用单重积分构造5重积分，似乎不如用单重＋二重＋二重积分构造5重积分精确。