楼主: wuyou136

[我分享] 基于全局优化工具箱求解非线性方程组数值解(无需初值)

  [复制链接]

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-13 20:34:35 | 显示全部楼层
kastin 发表于 2017-3-13 15:22
这几天比较忙,没来得及回复,见谅。

一般我用Mathematica用于精度验证,因为它相当于matlab的符号函数版 ...

谢谢,你在论坛都回答都非常有质量,我也是很希望能得到你的建议。:)  

对于方程求解的精度问题,是否只需将求得的解代回到原方程就可以了,也就是函数返回的 fval,而不需要再与 mathematica 的解析解进行对比?因为 double 型最多精确到小数点后16位左右,所以这种数值方法最高是可以精确到小数点后16位的,比如对于这里的贝塞尔函数,fval 在10^-16量级,已经达到了 double 型数值计算的最高精度,也就是说 GlobalSolve 通过调节误差控制参数是可以达到16位精度的,不知道这样理解是否正确?


另外你举例的高震荡函数来对收敛速度和收敛精度进行评估我还是没有理解,这怎么对方程求解进行评估?是指例如与前100个解析根进行对比吗?

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-13 22:17:38 | 显示全部楼层
wuyou136 发表于 2017-3-13 20:34
谢谢,你在论坛都回答都非常有质量,我也是很希望能得到你的建议。  

对于方程求解的精度问题,是否只 ...

谢谢你这么看重~

之前我也曾想过如何求解方程在某个范围内的所有根,但是由于没有一般的全局搜索算法,所以作罢。没想到你恰好解决了这个问题,所以很欣喜,本打算自己来进行各种测试以便让你能进一步优化和提高该函数的功能,只是因为版本问题无法测试,所以只是将自己的一点想法(不一定对,因为没有实际测试)说明一下。

说到fsolve对于初值的敏感性我是深有体会,这种敏感性会导致一些振荡根出现遗漏,之前我在计算余切函数的有理级数的部分和的根的时候,出现这个问题:

右端用有限项(比如n=6)表示一个有理函数的和,用f(x)表示之,那么可以方程 f(x)=a 的根可以用反余切函数近似表示。但是当时为了考察这种近似程度与n的关系,所以数值求解f(x)=a的根。这是一个振荡曲线,可以试试你的函数是否遗漏,或是精度如何。

精度按照fval来计算也可,但这不太容易判断与精确解的近似位数多少。

关于收敛速度和精度的评估比较复杂,其实也没必要分析,因为这是用的优化搜索方法,内部算法对于不同问题最终选择的分支是不同的。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-14 17:09:38 | 显示全部楼层
kastin 发表于 2017-3-13 22:17
谢谢你这么看重~

之前我也曾想过如何求解方程在某个范围内的所有根,但是由于没有一般的全局搜索算法,所 ...

谢谢建议。我尝试计算了一下你给的这个函数,可以通过调节 MultiStart 的初值点数,来增强搜索能力找寻所有的根。

在 n 的值比较小的时候,找寻所有根比较容易。但是一旦 n 的值增加,要想不遗漏根的话,初值点数要给的很大才行,并且它漏掉的根比较有规律,都是最靠近原点的那个根,我觉得这也是 GlobalSolve 后续可以改进的地方。

下面附上我测试的一些数据,对于不同的 n 值,如果要求得所有的根,GlobalSolve 所耗费的时间:

  1. a = 1; n = 10; lb = -n-5; ub = n+5;
  2. syms x
  3. fun = sum(1/x+2*x./(x^2-(1:n).^2))-a;
  4. fun = matlabFunction(fun);
  5. fplot(fun,[lb,ub]); hold on
  6. %%
  7. x = linspace(lb,ub,2e3); y = fun(x);
  8. [~,locs] = findpeaks(-abs(y));
  9. h = plot(x(locs),-0.5,'bo'); %蓝色标记代表 findpeaks 找到的根
  10. %%
  11. tic
  12. [~,~,~,ms,mf] = GlobalSolve(fun,1,lb,ub,100);
  13. toc
  14. plot(ms,zeros(size(ms)),'ro') %红色标记代表 GlobalSolve 找到的根
复制代码


n = 10
时间已过 2.036699 秒。

n = 30
时间已过 5.845213 秒。

n = 70
时间已过 25.255925 秒。

上面都是以默认精度进行计算的,fval 在 10^-5 量级,我把计算精度调至最高,发现 fval 是在 10^-11 量级,相应的计算时间也会增加。

n = 20

n = 20

n = 50

n = 50

n = 100

n = 100

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-14 17:25:23 | 显示全部楼层
本帖最后由 kastin 于 2017-3-14 17:28 编辑
wuyou136 发表于 2017-3-14 17:09
谢谢建议。我尝试计算了一下你给的这个函数,可以通过调节 MultiStart 的初值点数,来增强搜索能力找寻所 ...

目前看来效果还是很不错的,不过时间消耗还是比较大的。

我有几个想法:
1. 是不是lb和ub相应缩小也会降低搜索时间?其实可以先通过估计,然后写一个合理的区间,这样可能会减少很多时间的消耗。比如上面的根如果a不是很大,区间可以设为[-1.5n,1.5n]
2. 如果在区间缩小(但保证靠近0那个点在区间内),看看是否还会遗漏,从而找出原因并解决。我怀疑靠近零点的根遗漏是因为容差原因,因为方程中出现1/x的项,所以这一项会很大从而掩盖了后面的小项的有效数字。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-14 19:10:21 | 显示全部楼层
kastin 发表于 2017-3-14 17:25
目前看来效果还是很不错的,不过时间消耗还是比较大的。

我有几个想法:

缩小 lb 和 ub 可以降低搜索时间,就相当于在单位区间内增加初值点数。我觉得你分析的很对,漏掉靠近原点的根应该是因为 1/x 这一项的原因。我用来计算零点的话,就不会出现漏掉靠近原点的根的情况。之前也是想到缩小搜索区间来重新计算靠近原点的根,刚尝试了下,缩小搜索区间的话,也没有漏掉靠近原点的根了。

要提升 GlobalSolve 搜索效率,很重要的一点就是估算 lb 和 ub,这个要根据实际情况来了,lb 和 ub 给的越精确,GlobalSolve 的搜索效率也就越高。

论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-14 21:59:32 | 显示全部楼层
wuyou136 发表于 2017-3-14 19:10
缩小 lb 和 ub 可以降低搜索时间,就相当于在单位区间内增加初值点数。我觉得你分析的很对,漏掉靠近原点 ...

既然这样,如果可以事先确定根的区间(包括人为给定区间),那么能有效提高效率和降低遗漏;对于复杂函数如果很难直接判断所有根的区间,不知道是否能够通过二分搜索根的区间来调整GlobalSolve的智能性。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-16 17:14:04 | 显示全部楼层
本帖最后由 wuyou136 于 2017-3-16 17:21 编辑
kastin 发表于 2017-3-14 21:59
既然这样,如果可以事先确定根的区间(包括人为给定区间),那么能有效提高效率和降低遗漏;对于复杂函数 ...

抱歉回复晚了。还是以前面的震荡曲线为例,现在想将区间从 [-30 30] 由算法缩小到 [-10 10],但是二分法我没想明白怎么做可以将这个区间判断出来。
如果要有什么办法,能够判断出来[-30 -10]这一段函数值与 -30 这个端点函数值符号都一致的话那就好了,可以缩小区间同时也可以应用到多元函数里面。

M.jpg



论坛优秀回答者

退役版主

2206 麦片

财富积分



46

主题

2115

帖子

303

最佳答案
  • 关注者: 160
发表于 2017-3-16 19:21:20 | 显示全部楼层
wuyou136 发表于 2017-3-16 17:14
抱歉回复晚了。还是以前面的震荡曲线为例,现在想将区间从 [-30 30] 由算法缩小到 [-10 10],但是二分法我 ...

我说的二分法就是类似函数fzero里面的算法,通过符号来判断根的区间(但如果有无穷间断点,或者切点是根,这个方法就会失效,需要进行特殊处理)。

主要思路是对粗糙的估计区间(或给定的区间,比如[-30,30])分成两个区间(对半分,或按照黄金分割比来分),然后分别在这两个区间内执行GlobalSolve,如此类似递归下去。

如果GlobalSolve对于没有根的区间能很快判断出结果,那么这种方法可以极大提高效率(排除大部分无根区间)——但如果GlobalSolve对无根区间的搜索时间和有根的区间差不多,这种方法就没有什么意义了。

一旦上述方法失效,唯一的办法就是你所说的判断某个区间的函数值符号是什么样的,能做到这一点的就是所谓的“区间算术”了,这是数学领域一个比较热门的话题。

目前有MATLAB工具箱INTLAB可实现这个功能,之前我在本论坛分享过这个工具箱(传送门),不过那是多年前的版本,现在官网上的应该是新版的,不过要收费了(50欧元起)。所以幸好当初我分享了,看我分享后有1万多的下载量就知道它帮助了多少人。

官网地址在下面,里面介绍了工具箱的特点和功能,例子演示以及各种应用,还有相关理论的论文(Further reading下面)。
  1. http://www.ti3.tu-harburg.de/rump/intlab/index.html
复制代码
实际上,有了区间算术功能,非线性方程所有根的求解就不是问题了——果然这个工具箱已经实现了非线性方程所有根的求解问题。对于单变量和躲变量方程求根也有涉及。自然,全局优化的问题也能得到部分解决——它还能验证全局最优问题和约束,参数识别(这个也是需要全局搜索,并且对初值很敏感)。

总的来说,这个工具箱还不错,只是我还没怎么用到它,不太熟悉,但我认得好东西,所以当初就下载分享了。你可以研究一下,应该有不少收获。
请点击"回复此楼",否则我将无法收到回帖提醒。
问题如果比较复杂或较难,请邮箱联系kimist@qq.com

论坛优秀回答者

权威

3630 麦片

财富积分



20

主题

3779

帖子

769

最佳答案
  • 关注者: 431
 楼主| 发表于 2017-3-16 19:51:16 | 显示全部楼层
kastin 发表于 2017-3-16 19:21
我说的二分法就是类似函数fzero里面的算法,通过符号来判断根的区间(但如果有无穷间断点,或者切点是根 ...

在你的答案里面总是能看到新的东西。我刚才测试了一下,GlobalSolve 对于无根区间的搜索,在同等大小区间的情况,比有根区间消耗的时间甚至要高3倍。二分法应该用不上了,也说明了区间优化很有必要。谢谢你提供的链接,我好好研究一下。:)

新手

51 麦片

财富积分


050


16

主题

105

帖子

3

最佳答案
  • 关注者: 1
发表于 2017-3-17 09:16:07 | 显示全部楼层
楼主这个函数怎么使啊,为什么我在用这个函数验证时,出现下面这样的错误呢?
QQ截图20170317092251.png
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表