专栏首页matlab爱好者数值优化—三种复杂函数数值积分方法实例演示

数值优化—三种复杂函数数值积分方法实例演示

过冷水前段时间在学习过程中遇到了一道难题,需要计算

在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 偏差会更加明显。

过冷水也是在实践过程中发现各种问题,觉得有点意思就把自己的问题拿出来和大家共同分享,促进过冷水和读者之间相互交流以及知识增长。过冷水本期和大家分享的知识点就这些。在实际应用中三种方法的可行性都比较高,能够解决复杂函数积分的问题,实际在解决数学问题中方法是很多的,蒙特卡洛算法、多项式应用较广,感兴趣的可深入研究。

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

本文分享自微信公众号 - matlab爱好者(matlabaihaozhe)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-02-09

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 统计分布讲解

    随机现象中,变量的取值是不确定的,称之为随机变量。描述随机变量取值概率的函数称为概率分布。对于随机变量,通常主要关心它的两个主要数字特征:数学期望用于描述随机变...

    matlab爱好者
  • 非线性方程组求解迭代算法&图像寻初始值讲解

    前段时间过冷水在学习中遇到了一个解非线性方程组的问题,遇到非线性方程组的的问题过冷水果断一如既往、毫不犹豫的 fsolve()、feval()函数走起,直到有人...

    matlab爱好者
  • 二维热导方程Matlab数值解案例

    本次和大家一起来看看如何根据一维热传导有限差分法的思想求解二维热传导方程进行求解。

    matlab爱好者
  • QQ Mac版 Touch Bar功能设计

    腾讯ISUX
  • 陆奇最新演讲:如何成为一个优秀的工程师

    大数据文摘
  • 百度陆奇:如何成为一个优秀的工程师

    7月11日,陆奇出席百度内部Engineering Leadership Talk。作为计算机科学博士及优秀的管理者,他提出的五点要求,对每一位工程师都适用。

    华章科技
  • 陆奇最新演讲:如何成为一个优秀的工程师

    一位工程师,如何才能称得上优秀?除了写得一手好Code,什么样的工作态度和方法才是一个优秀工程师的必备?

    疯狂的技术宅
  • 复制状态与变量记录表 | performance_schema全方位介绍

    不知不觉中,performance_schema系列快要接近尾声了,今天将带领大家一起踏上系列第六篇的征程(全系共7个篇章),在这一期里,我们将为大家全面讲解p...

    沃趣科技
  • 如何使用Grouper2来查找活动目录组策略中的漏洞

    Grouper2是一款针对AD组策略安全的渗透测试工具,该工具采用C#开发,在Grouper2的帮助下,渗透测试人员可以轻松在活动目录组策略中查找到安全相关的错...

    FB客服
  • QtCreator-启用/关闭FakeVim模式

    Qt君

扫码关注云+社区

领取腾讯云代金券