优化积分方法在Matlab中的实践



早前在上数学分析课的时候，研究了matlab上的不同积分方法的速度以及误差率。

分别是：

一般积分法

示性函数法

蒙特卡洛法

并在这一过程中发掘出了Matlab新带的Mupad展现的强大画图功能（之前Matlab的画图效果就不吐嘲了）

一并记录如下。

以下是习题12.3部分例题基于matlab2014b的解答编写

(鉴于matlab2014b加入了带有便捷功能的mupad模块，额外试用mupad进行处理.)
(默认数值显示方式 format rat）

求下面曲面所围区域的体积

例1，z=x^{2}+y^{2},abs{x}=1,abs{y}=1,z=0.

(1)绘制积分区域

在命令窗口打开mupad窗口

>>mupad

在mupad命令窗口绘制积分区域

>>plot(plot::Implicit3d(x^2 + y^2 =z,x=-2..2, y=-2..2, z=0..2),     plot::Implicit3d(abs(x)=1, x=-2..2, y=-2..2, z=0..2),     plot::Implicit3d(abs(y)=1, x=-2..2, y=-2..2, z=0..2),     plot::Implicit3d(z=0, x=-2..2, y=-2..2, z=0..2))

（2）求解积分

在mupad窗口输入

4*int(int(int(1, z=0..x^2+y^2), y=0..1), x=0..1)

得

8/3

下面分别用3种matlab命令计算其积分，并试看其误差

解法1 一般累次积分

>> ticsyms x y z4*int(int(int(1,z,0,x^2+y^2),y,0,1),x,0,1)toc ans = 8/3 时间已过 0.059050 秒.

计算用时0.06秒，误差为0.

解法2 试性函数

>> tic f1=@(x,y,z)(x.^2+y.^2>=z);v1=triplequad(f1,0,1,0,1,0,sqrt(2));tocv1时间已过 1.053799 秒。

>> vpa((2/3-v1)/(2/3),4)  ans = 0.01493

计算用时1秒，误差约百分之一

解法3 蒙特卡罗方法

>>  tic X=unifrnd(0,1,1,10000000); Y=unifrnd(0,1,1,10000000); Z=unifrnd(0,sqrt(2),1,10000000); B=mean(X.^2+Y.^2>=Z); v=sqrt(2)*B; toc时间已过 0.636313 秒。>> vpa((2/3-v)/(2/3),4) ans = 0.01534

计算用时0.6秒，误差约百分之一

例2，z=x*y,x+y+z=1,z=0.

(1)绘制积分区域

在命令窗口打开mupad窗口

>>mupad

在mupad命令窗口绘制区域

>>plot(plot::Implicit3d(z=x*y,x=-2..2, y=-2..2, z=0..2),     plot::Implicit3d(x+y+z=1, x=-2..2, y=-2..2, z=0..2),     plot::Implicit3d(z=0, x=-2..2, y=-2..2, z=0..2))

(局部放大)

这里注意，Mupad界面提供了图像颜色，网格，透明度等直接设置的按钮。

（2）求解积分

在mupad窗口输入

int(int(int(1, z=0..x*y), y=0..(1-x)/(1+x)), x=0..1)+int(int(int(1, z=0..1-x-y), y=(1-x)/(1+x)..1-x), x=0..1)

得

ln(16) - ln(64) + 17/12

下面分别用3种matlab命令计算其积分，并试看其误差.

解法1 一般累次积分

>> ticsyms x y zint(int(int(1, z,0,x*y), y,0,(1-x)/(1+x)), x,0,1)+int(int(int(1, z,0,1-x-y), y,(1-x)/(1+x),1-x), x,0,1);toc ans = log(16) - 11/4 时间已过 0.090197 秒。

计算用时0.09，误差为O.

解法2 示性函数

>> tic f1=@(x,y,z)((z)<=(x*y)).*(x+y+z)<=1.*(z)>=0;v1=triplequad(f1,0,1,0,1,0,0.5);toc时间已过 2.549438 秒。>> vpa((v1-(log(16) - 11/4)),4) ans = 0.0073

误差值仅千分之七，但是相对误差

>> vpa((v1-(log(16) - 11/4))/(log(16) - 11/4),4) ans = 0.3254

综上
计算用时2.5秒，误差百分之三十.

解法3 蒙特卡罗方法

>>  tic X=unifrnd(0,1,1,10000000); Y=unifrnd(0,1,1,10000000); Z=unifrnd(0,1,1,10000000); B=mean((Z<=X.*Y).*(X+Y+Z<1)); toc时间已过 0.745913 秒。>> vpa((B-(log(16) - 11/4)),4) ans = 0.0078>> vpa((B-(log(16) - 11/4))/(log(16) - 11/4),4) ans = 0.3434

 同理可见，计算用时0.7秒，误差百分之三十。

例3，x^{2}+z^{2}=1,abs{1}+abs{1}=1.

(1)绘制积分区域

在命令窗口打开mupad窗口

>>mupad

在mupad命令窗口绘制区域

plot(plot::Implicit3d(x^2 + z^2 - 1,x=-2..2, y=-2..2, z=-2..2),     plot::Implicit3d(abs(x)+abs(y)=1, x=-2..2, y=-2..2, z=-2..2))

（2）求解积分

在mupad窗口输入

4*int(int(int(1, z=0..sqrt(1-x^2)), x=y-1..1-y), y=0..1)

得

2*PI - 8/3

下面分别用3种matlab命令计算其积分，并试看其误差.

解法1 一般累次积分

>> ticsyms x y z4*int(int(int(1, z,0,sqrt(1-x^2)), x,y-1,1-y), y,0,1);toc时间已过 0.052944 秒。>> ans ans = 2*pi - 8/3

计算用时0.05，误差为O.

解法2 示性函数

>> tic f1=@(x,y,z) (x.^2+z.^2<=1).*(abs(x)+abs(y)<=1);v1=triplequad(f1,-1,1,-1,1,-1,1);toc时间已过 10.021040 秒。下面检测其误差>> v=2*pi-8/3;vpa((v1-v)/v,4) ans = -6.252e-6

可见，计算用时10秒，误差十万分之六

解法3 蒙特卡洛方法

>> tic X=unifrnd(0,1,1,10000000); Y=unifrnd(0,1,1,10000000); Z=unifrnd(0,1,1,10000000); B=mean((X.^2+Z.^2<=1).*(abs(X)+abs(Y)<=1)); v1=8*B; toc时间已过 0.748515 秒。下面检验其误差>> v=2*pi-8/3;vpa(((v1-v)/v),6)  ans = 0.0000831074

可见，计算用时0.74秒，误差万分之八

综上所述，更一般的，蒙特卡洛方法的积分在更高维有更加显著的效果，所谓蒙特卡洛方法，通俗意义上讲，就是随机撒点，根据大数定律，逼近真实值，在其它算法中有更优化的体现，关于此，希望我能有进一步研究和体会。