数值优化—复杂函数重积分计算方法实例演示

发布于
过冷水前段时间做了一篇数值优化—三种复杂函

既然要的是数值解，为何还使用符号解？能坐车进城，就决不骑摩托车。复杂函数用数值积分函数quad(f(x),xmin,xmax)完美求解，perfect!不巧的是疑难杂症都让过冷水碰上了,在原问题的基础上需要解决这么个问题:

之前求解的f(x)是数值解，没有函数表达式，不能再次使用int函数，更不可能用int(int(f(x)))求解，这是不合理的。私以为可以坐公交车到城里KTV去蹦迪，谁料新型冠状病使得客车停运了，这时摩托车就实用了。这种情况下就需要使用函数拟合方法、蒙特卡洛算法解决二次或多次积分的问题。使用两种方法需要注意的问题有：

（1）：若是图像变化趋势复杂没有合适的拟合函数，在三重积分以上情况下就只能使用蒙特卡洛算法；

（2）：能够使用函数近似替代法尽量使用函数替代法，蒙特卡洛算法涉及到的随机取点求积分值的语句运算量较大，耗时；

（3）：本文中蒙特卡洛算法和之前方法不一样之处在于面积上下限不一样，需要灵活调整。判断函数表达式和之前的语句不一样。具体调整方法代码体现在多重积分计算时要将概率稳定性和运行时间综合考虑选取实验次数。

蒙特卡洛代码如下：

代码语言：javascript
复制
clear all
warning off
feature jit off%程序加速；
X=linspace(0.12,1,50);
n=5000;
f=@(x)((x.^4.*exp(x))./(exp(x)-1).^2);%定义原始积分项函数；
for j=1:length(X);
    xx=min(X)+(X(j)-min(X)).*rand(1,n);%设置随机点X轴范围；
    %arrayfun：将积分函数作用于每个变量中，输出一组值；
    %quad：求函数数值解；
    %  y=arrayfun(@(x)(quad(f,0.12,x)),xx)：由于xx很大，该语句运行耗时较长；
    y=arrayfun(@(x)(quad(f,0.12,x)),xx);
    y1=100*xx.^-3.*y;
    y2=min(y)+(40-min(y)).*rand(1,n);%设置随机点Y轴范围；
    points=find(y2<=y1);
    k=length(points);
    Y1(j)=(X(j)-min(X))*(40-min(y))*(k/n)
end

再来看一下函数近似拟合替换程序：

代码语言：javascript
复制
syms x
X=linspace(0.12,1,50);
f=@(x)((x.^4.*exp(x))./(exp(x)-1).^2);%定义原始积分项函数；
y=arrayfun(@(x)(quad(f,0.12,x)),X);
y1=100*X.^-3.*y;
p=polyfit(X,y1,5);%多项式拟合；
f1= p(1)*x.^5 + p(2)*x.^4 + p(3)*x.^3 + p(4)*x.^2 + p(5)*x + p(6);
F=int(f1);%符合计算所求原函数；
Y2=feval(inline(F),X)-feval(inline(F),min(X));%feval:符号函数直接求解；

二过冷水就开心的骑着小摩托准备进城了，走在路上碰到好友，他说天这么冷骑什么车进城啊！我开奔驰送你进城。有奔驰不坐是傻子，于是过冷水开心的准备坐奔驰。

Matalab除了不能生孩子外，不能做的事还有很多，当然能做的事也很多。Matalab自带函数能力很齐全，我们常用的只是一角，很多你所能想到的麻烦问题都有特定的函数去解决类似问题，你只需要熟练运用即可。可以把Y(x)看成是一个二重积分用integral2解决我们的难题。

代码如下：

代码语言：javascript
复制
X=linspace(0.12,1,50);
f=@(x,t)(100.*x.^-3.*t.^4.*exp(t))./(exp(t)-1).^2;
tmax=@(c)c;%定义为一个变量X，表示第二项积分变量；
for i=1:length(X)
  Y3(i)=integral2(f,min(X),X(i),min(X),tmax) ;
end

哟，真香！别看代码这么简单，过冷水失败了好多次才运行正常，这是什么原因呢？我们下面这个函数代码：

代码语言：javascript
复制
X=linspace(0.12,1,50);
f=@(t,x)(100.*x.^-3.*t.^4.*exp(t))./(exp(t)-1).^2;
tmax=@(c)c;%定义为一个变量X，表示第二项积分变量；
for i=1:length(X)
  Y33(i)=integral2(f,min(X),X(i),min(X),tmax) ;
end

disp(Y3(1:5))

disp(Y33(1:5))

0    0.1117    0.3814    0.7473    1.1751

0    0.1409    0.5902    1.3877    2.5731

程序差别在定义函数时@(x,t)、@(t,x)两者的运行结果差别很大。Why?现在简单给讲一下integral2函数

二重积分涉及到积分先后顺序和积分限的问题在此不讲，该函数的积分限和积分先后顺序是固定好的，无须讨论。我们来看引起前述差别的原因是：

举个具体案例f(x,t)=t;

函数二重积分代码如下：

代码语言：javascript
复制
clear all

f1=@(x,t)t;

f2=@(t,x)t;

tmax=@(c)c;

z=3;

y1=integral2(f1,0,z,0,tmax);

y2=integral2(f2,0,z,0,tmax);

disp(y1);

disp(y2);

4.5000

9.0000

二重积分函数的调用注意项讲解完毕，回到正题。现在有三种方法求得我们想要的结果。

代码：

代码语言：javascript
复制
figure1 = figure;

axes1 = axes('Parent',figure1);

hold(axes1,'on');

%Y1值计算较耗时，若嫌耗时可以Y2换Y1使程序正常运行；

plot1 = plot(X,[Y1' Y2' Y3'],'Parent',axes1);

set(plot1(1),'DisplayName','蒙特卡洛算法','LineWidth',1.5);

set(plot1(2),'DisplayName','函数拟合替代');

set(plot1(3),'DisplayName','

interal2

','LineWidth',1);

xlabel('

x$$&#39;,&#39;Interpreter&#39;,&#39;latex&#39;);
ylabel(&#39;

f(x)$$','Interpreter','latex');

title('不同方法计算结果比较');

set(axes1,'FontSize',16,'LineWidth',1.5);

legend1 = legend(axes1,'show');

set(legend1,'Position',[0.671774724027605 0.779508549971586 0.124450949196585 0.136111107385821],'Interpreter','latex');

由图可知三种方法的差别很小，整体来看计算结果都较好，蒙特卡洛算法存在波动，说明我们实验次数取得不够大。函数拟合替换和二重积分函数基本一致很难区分，暗示在多重积分计算中函数拟合有较好的潜在的应用价值。

过冷水在学习过程中很明显感觉到用Matalab解决一个问题时需要的知识很杂。在用函数近似替换和蒙特卡洛算法时，和数值积分是交叉使用的，为了使得程序运行简单还混用其它方法，来减少程序运行时间，说明在用Matalab解决一个复杂的问题时，会涉及到各方面知识，综合性很高，Matalab基础知识要广博才能较好的解决一个问题。读者跟随Matalab爱好者平台我们共同学习，日日精进。

如需转载，请在公众号中回复“转载”获取授权，未经授权擅自搬运抄袭的，必将追究其责任！

#int#perfect