数值优化—三种复杂函数数值积分方法实例演示
过冷水前段时间在学习过程中遇到了一道难题，需要计算
在0.1~1 区间上的值，初步看该方程的积分项比较复杂不易给出原函数。用MATLAB也无法直接求出原函数。自然而然就想该函数如何在不求积分项原函数的情况下计算出积分项的具体值。在抓耳挠腮之际想起了公众号的一篇推文：蒙特卡洛法应用。可以直接求函数指定区间的面积，相当于求积分。蒙特卡洛算法求面积示意图如下：
在该思路的启发下过冷书立刻实践给出了对应的代码，求得函数解。原代码如下：
clear all
warning off
feature jit off
X=linspace(0.12,1,50);
n=100000; % 随机生成n个点
for i=1:length(X);
x=(1/X(i)).*rand(1,n);y0=5.*rand(1,n);
y=(x.^4.*exp(x))./(exp(x)-1).^2;
k=0;
points=find(y0<=y);
k=length(points);
S(i)=(1/X(i))*5*(k/n);
end
Y1=100*X.^3.*S;
X=0.7;x=(1/X)*rand(1,n);y0=2*rand(1,n);
y=(x.^4.*exp(x))./(exp(x)-1).^2;
points=find(y0<=y)；
x1=linspace(0,2,50);y1=(x1.^4.*exp(x1))./(exp(x1)-1).^2
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
% 创建 plot
plot(x1,y1,'LineWidth',1.5,'Color',[1 0 0]);
plot(x(points),y0(points),'MarkerFaceColor',[1 0 1],'MarkerEdgeColor','none','MarkerSize',3,'Marker','o','LineStyle','none','Color',[0 0 0]);
% 创建 rectangle
rectangle('Parent',axes1,'Position',[0 0 1/X 2],'EdgeColor',[0 0.447058826684952 0.74117648601532]);
% 创建 text
text('Parent',axes1,'FontSize',16,'FontName','Helvetica-Narrow','Interpreter','latex','String','（0.7,1）','Position',[1.4 2]);
text('Parent',axes1,'FontWeight','bold','FontSize',16,'FontName','Times New Roman','Interpreter','latex','String','$$S_{int}=P.S_{rectangle}$$','Position',[1.6 1]);
% 创建 label
xlabel('x');ylabel('y');
% 创建 title
title('蒙特卡洛求复杂函数积分数值','FontWeight','bold');
% 设置坐标轴属性
xlim(axes1,[0 2]);
ylim(axes1,[0 3]);
set(axes1,'FontSize',12,'FontWeight','bold','LineWidth',1.5);之前同伴在给大家分享蒙特卡洛算法计算面积的时候，过冷书吐槽算法算复杂积分精度不好，涉及到循环函数麻烦，可以用泰勒公式将复杂函数转换为多项式的形式，多项式原函数很容易求。既然牛已经吹出去了，现行下有实际问题，我们不妨用多项式拟合替代原函数试试看，证明方法论是否可行。历经多次失败&偶然成功，使得我关于用多项式替代复杂函数得到以下结论。
(1) 多项式替代复杂函数局限性很大，很多复杂函数基本上做不到使得多项式在较宽自变量区间吻合复杂函数；
(2) 多项式表达式只是可以求原函数，级数多了以后表达式也是很麻烦的；
(3）多项式计算积分的时候涉及到y负轴求积分计算结果存在负值，这是用程序自动计算积分的弊端，需要自行调整程序。
说了这么多多项式替代复杂函数的不好，我们只认为方法论是可行的，在我们研究的在0.1~1 区间上的是可以用多项式替代复杂函数。让我们来看看具体存在什么问题。首先我们得给出吻合复杂函数的多项表达式
由图可知多项式和原函数的符合程度较好，我们可以认为用多项式替换原函数是可行的。然后求多项式的原函数正常计算即可。代码如下：
warning off
feature jit off
X=linspace(0.12,1,50);x=1./X;y=(x.^4.*exp(x))./(exp(x)-1).^2;
p=polyfit(x,y,5);
f= p(1)*x.^5 + p(2)*x.^4 + p(3)*x.^3 + p(4)*x.^2 + p(5)*x + p(6);
syms x
f= p(1)*x.^5 + p(2)*x.^4 + p(3)*x.^3 + p(4)*x.^2 + p(5)*x + p(6);
F=int(f);
X=linspace(0.12,1,50);Y2=100*X.^3.*(feval(inline(F),1./X)-feval(inline(F),0));%feval:符号函数直接求解
figure1 = figure;
% 创建 axes
axes1 = axes('Parent',figure1);
hold(axes1,'on');
% 使用 plot 的矩阵输入创建多行
plot1 = plot([x' x'],[y' f']);
set(plot1(1),'DisplayName','$$\frac{x^4.exp(x)}{(exp(x)-1)^2}$$','MarkerFaceColor',[1 0 0],'Marker','o','LineStyle','none','Color',[1 0 0]);
set(plot1(2),'DisplayName','$$p_1*x.^5 + p_2*x.^4 + p_3*x.^3 + p_4*x.^2 + p_5*x + p_6$$','LineWidth',1.5,'Color',[0 0 1]);
% 创建 label
xlabel('$$x$$','FontAngle','italic','Interpreter','latex');
ylabel('$$f(x)$$','FontAngle','italic','Interpreter','latex');
% 创建 title
title('多项式拟合复杂函数');
% 设置坐标轴属性
xlim(axes1,[0.699454538886681 8.69945453888668]);
ylim(axes1,[0.732361963190186 5.23236196319018]);
set(axes1,'FontSize',14,'FontWeight','bold');
% 创建 legend
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.319778619498812 0.210874316368273 0.339180486457926 0.123570128076172],'Interpreter','latex','FontSize',14);两种方法的计算结果都不能作为标准参考来验证计算结果是否好，可以采用MATLAB符号运算法求得函数值，符号计算法较为简单，代码如下：
warning off
feature jit off
syms x
X=linspace(0.12,1,50);
for i=1:length(X);
f=(x.^4.*exp(x))./(exp(x)-1).^2;
y(i)=double((int(f,0,1/X(i))));%double:符号转数值
end
Y3=100*X.^3.*y;
%%三种方法方程值比较代码
Y=[Y1',Y2',Y3']
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot1 = plot([x',x',x'],Y,'LineWidth',1.5,'Parent',axes1);
set(plot1(1),'DisplayName','蒙特卡洛算法','Color',[1 0 0]);
set(plot1(2),'DisplayName','多项式替代','Color',[0.929411768913269 0.694117665290833 0.125490203499794]);
set(plot1(3),'DisplayName','符号积分','Color',[0 0 1]);
% 创建 label
xlabel('$$x$$','FontAngle','italic','FontName','Times New Roman', 'Interpreter','latex');
ylabel('$$Y$$','FontWeight','bold','FontAngle','italic','FontName','Times New Roman','Interpreter','latex');
% 创建 title
title('不同方法求解析值比较');
% 设置坐标轴属性
xlim(axes1,[0 10]);
ylim(axes1,[1.65644171779141 46.6564417177914]);
set(axes1,'FontSize',14,'FontWeight','bold','LineWidth',1.5);
% 创建 legend
legend(axes1,'show');三种方法得到的函数值比较如图
根据图像分析可得如下结论：
（1）：三种方法计算的函数值大致走势一致，三种方法互证可行性，自变量较大是三者一致性较好，自变量减小时，差别明显；
（2）蒙特卡洛算法和符号算法整体吻合程度较高，在精度要求不是非常高的计算中可以用蒙特卡洛方法思路解决问题。蒙特卡洛算法在自变量较小时存在数值明显的波动，概率法求值很容易出现波动，这说明我们概率矩阵有可能设置不合理，或者取点次数太小，概率不稳定，关于改进方法在此不详述；
（3）多项式替换法在自变量较较小时误差较大，但是我们绘制多项式和原函数的比较图时，根本无法看出两者差别，也说明了多项式替换可用但是存在较大误差，使用时要进行多方分析，实际若x:0~2 偏差会更加明显。
过冷水也是在实践过程中发现各种问题，觉得有点意思就把自己的问题拿出来和大家共同分享，促进过冷水和读者之间相互交流以及知识增长。过冷水本期和大家分享的知识点就这些。在实际应用中三种方法的可行性都比较高，能够解决复杂函数积分的问题，实际在解决数学问题中方法是很多的，蒙特卡洛算法、多项式应用较广，感兴趣的可深入研究。
如需转载，请在公众号中回复“转载”获取授权，未经授权擅自搬运抄袭的，必将追究其责任！