声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 18396|回复: 23

[综合讨论] 非线性方程组-1stOpt与fsolve的比较!

[复制链接]
发表于 2007-7-12 15:22 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
huright 在帖子http://forum.vibunion.com/forum/thread-47883-1-1.html介绍了Matlab求非线性方程组fsolve的应用和实例。不过个人认为非线性方程组,最强最好用的当属目前的1stOpt了,基本无需编程,也不用猜初值,求解能力非MatLab的fsolve可比。下面第一例是上面介绍的实例6,系数是变化的;第二例是道看上去并不复杂的的非线性方程组,有兴趣的试试,看fsolve能否求出来,化多少时间猜初值?。
例1:
Constant H=0.32, Pc0=0.23, w=0.18;
LoopConstant k=[0:0.1:1];
Function Pc0+H*(1+1.5*(x1/W-1)-0.5*(x1/w-1)^3)-x2;
            x1-k*x2^0.5;

例2:四维非线性方程组
p0+p1*(1-exp(-(p2*(0)^p3)))=51.61;
p0+p1*(1-exp(-(p2*(9.78)^p3)))=51.91;
p0+p1*(1-exp(-(p2*(30.68)^p3)))=53.27;
p0+p1*(1-exp(-(p2*(59.7)^p3)))=59.68;

有用Matlab解出第二例答案的吗?

[ 本帖最后由 ChaChing 于 2009-9-6 11:09 编辑 ]
34.jpg

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2007-7-18 14:45 | 显示全部楼层
虽然找合适的初值是个问题,但是楼主也不要把MATLAB的fsolve说的一无是处。楼主这篇帖子明显带有贬低其它软件的倾向了。
bainhome曾经在下面帖子里说过一些话,基本上我的观点和他一致。http://forum.vibunion.com/thread-17903-1-1.html
根据MATLAB的提示,试验了一会找到了一个合理的初值。以下是
MATLAB解出来的结果。
p =
  51.609999999999999
  -0.465222188130621
  -0.053730715604489
   0.976132921502588
楼主可以试验一下上面的结果和1stOpt的结果相比差多少。
MATLAB求解的代码:
  1. function p=myfun1
  2.    p0 = [51.61;-1;0;1];
  3.    options=optimset('MaxFunEvals',20000,'MaxIter',2000);
  4.     [p,fval] = fsolve(@f,p0,options) ;
  5.    
  6. function F=f(p)
  7. F=[p(1)+p(2)*(1-exp(-(p(3)*(0)^p(4))))-51.61;
  8. p(1)+p(2)*(1-exp(-(p(3)*(9.78)^p(4))))-51.91;
  9. p(1)+p(2)*(1-exp(-(p(3)*(30.68)^p(4))))-53.27;
  10. p(1)+p(2)*(1-exp(-(p(3)*(59.7)^p(4))))-59.68;];
  11. end
  12. end
复制代码
推荐两个地址:
http://www.simwe.com/forum/thread-787570-1-1.html
http://www.simwe.com/forum/thread-788981-1-1.html

[ 本帖最后由 rocwoods 于 2007-7-18 14:59 编辑 ]

评分

1

查看全部评分

 楼主| 发表于 2007-7-18 15:31 | 显示全部楼层
初值猜的很精准啊,佩服!不过对于不少人而言,“功力”是不容易达到这个水准的。

比较一下,没什么坏处。Matlab是数学里的老大,Maple之类的也不能与之争锋,何况其它?

其实还有一组解,再试下。(如果是参照1stOpt的结果来定初值,就没必要试了)
发表于 2007-7-18 16:04 | 显示全部楼层
p =
  51.609999999999999
  -0.134520891677992
  -0.241198737477775
   0.693446385900985
第二组值,说实在的,不知道是我运气好还是怎么着,还是受1stOpt初值随机给的启发。我在fsolve中也随机给初值。代码如下:

  1. function p=myfun1
  2.    p0 = [51.61;rand;unifrnd(-1,1);rand];
  3.    options=optimset('MaxFunEvals',20000,'MaxIter',5000);
  4.     [p,fval] = fsolve(@f,p0,options) ;
  5.    
  6. function F=f(p)
  7. F=[p(1)+p(2)*(1-exp(-(p(3)*(0)^p(4))))-51.61;
  8. p(1)+p(2)*(1-exp(-(p(3)*(9.78)^p(4))))-51.91;
  9. p(1)+p(2)*(1-exp(-(p(3)*(30.68)^p(4))))-53.27;
  10. p(1)+p(2)*(1-exp(-(p(3)*(59.7)^p(4))))-59.68;];
  11. end
  12. end
复制代码

试了两次就得到了上面的解。
其实,只要有了自己的一套算法语言,无所谓谁好谁坏。只有算法才是灵魂。本人喜爱MATLAB,同时我也很喜欢1stOpt这样优秀的国产软件。之所以喜欢,并不是因为它们的名字,而是它们是得力的工具。能够帮我解决我遇到的问题。同时能够提供给我们一个平台去实现自己的想法。废话也不多说了,有了算法,无所谓哪个软件好,哪个软件坏!

评分

1

查看全部评分

 楼主| 发表于 2007-7-18 16:18 | 显示全部楼层
赞一个:victory:
发表于 2007-7-18 16:24 | 显示全部楼层
我也还是持一贯的观点: 充分利用各个软件的优点才是最重要的. 不宜贬低某个软件,当然前提是都比较强.
拟合优化是1stOpt的强项,当然我们就应该利用它;但就象我们不能因为它绘图、矩阵计算不强,无法解微分方程等等,就把1stOpt贬低得一无是处一样。

猜测一下:看dingd 的回帖,似乎本身是很想知道Matlab是如何处理该问题的...

[ 本帖最后由 eight 于 2007-7-19 14:54 编辑 ]
 楼主| 发表于 2007-7-19 11:24 | 显示全部楼层
rocwoods的Matlab代码运行一下,大概能以1/10的概率求得正解。代码中如果初值都用随机值,即把 “p0 = [51.61;rand;unifrnd(-1,1);rand]; ”改为“p0 = [rand;rand;rand;rand]; ”,没有人为控制因素,成功的概率就更低些了。看样子,与1stOpt思路不一样。

1stOpt安装包里附带的一个方程组例子,也有点意思,大家可试试:x1至x9,共9个参数,要求必须大于等于0。老版的1stOpt成功率也不高,新版,参数设置好了,近100%成功率。

23.3037*x2+(1-x1*x2)*x3*(exp(x5*( 0.485- 0.0052095*x7- 0.0285132*x8))-1)- 28.5132;
- 28.5132*x1+(1-x1*x2)*x4*(exp(x6*( 0.116- 0.0052095*x7+ 0.0233037*x9))-1)+ 23.3037;
101.779*x2+(1-x1*x2)*x3*(exp(x5*( 0.752- 0.0100677*x7- 0.1118467*x8))-1)- 111.8467;
- 111.8467*x1+(1-x1*x2)*x4*(exp(x6*(- 0.502- 0.0100677*x7+ 0.101779*x9))-1)+ 101.779;
111.461*x2+(1-x1*x2)*x3*(exp(x5*( 0.869- 0.0229274*x7- 0.1343884*x8))-1)- 134.3884;
- 134.3884*x1+(1-x1*x2)*x4*(exp(x6*( 0.166- 0.0229274*x7+ 0.111461*x9))-1)+ 111.461;
191.267*x2+(1-x1*x2)*x3*(exp(x5*( 0.982- 0.0202153*x7-0.2114823*x8))-1)-211.4823;
- 211.4823*x1+(1-x1*x2)*x4*(exp(x6*(- 0.473- 0.0202153*x7+ 0.191267*x9))-1)+ 191.267;
x1*x3-x2*x4;
发表于 2007-7-19 11:57 | 显示全部楼层
用如下代码:

  1. function [x,x0]=myfun3
  2.    x0 = unifrnd(0,1,1,9);
  3.    options=optimset('MaxFunEvals',20000,'MaxIter',5000);
  4.     [x,fval] = fsolve(@f,x0,options) ;
  5.    
  6. function F=f(x)
  7. F=[23.3037*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.485- 0.0052095*x(7)- 0.0285132*x(8)))-1)- 28.5132;
  8. - 28.5132*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*( 0.116- 0.0052095*x(7)+ 0.0233037*x(9)))-1)+ 23.3037;
  9. 101.779*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.752- 0.0100677*x(7)- 0.1118467*x(8)))-1)- 111.8467;
  10. - 111.8467*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*(- 0.502- 0.0100677*x(7)+ 0.101779*x(9)))-1)+ 101.779;
  11. 111.461*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.869- 0.0229274*x(7)- 0.1343884*x(8)))-1)- 134.3884;
  12. - 134.3884*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*( 0.166- 0.0229274*x(7)+ 0.111461*x(9)))-1)+ 111.461;
  13. 191.267*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.982- 0.0202153*x(7)-0.2114823*x(8)))-1)-211.4823;
  14. - 211.4823*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*(- 0.473- 0.0202153*x(7)+ 0.191267*x(9)))-1)+ 191.267;
  15. x(1)*x(3)-x(2)*x(4)];
  16. end
  17. end
复制代码

实验了两次得到一组解:不过参数有大于0有小于0。有时间试试怎么得到全大于0的。
由于是随机生成迭代初始值,故x0=0.2217        0.9159        0.7538        0.0229        0.3783        0.7437        0.5186        0.5524        0.279
求出来的结果
x =
  Columns 1 through 4
   0.898509997792166   0.973967802005302  11.648633060796092  10.746159415319836
  Columns 5 through 8
   3.250825125902567   6.710547896595700  -8.764002776048228   1.251334798196427
  Column 9
  -0.525131249969468
发表于 2007-7-19 12:52 | 显示全部楼层
全部大于等于0的解
x =
  Columns 1 through 4
   0.899999952616886   0.449987471977047   1.000006482444359   2.000068541604411
  Columns 5 through 8
   7.999971440536423   7.999692684210950   5.000031275992702   0.999987723454627
  Column 9
   2.000052483507424

  1. function [x,x0]=myfun3
  2.    x0 = unifrnd(-2,2,1,9);
  3.    options=optimset('MaxFunEvals',10000,'MaxIter',5000);
  4.     [x,fval] = fsolve(@f,x0,options) ;
  5.    
  6. function F=f(x)
  7. F=[23.3037*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.485- 0.0052095*x(7)- 0.0285132*x(8)))-1)- 28.5132;
  8. - 28.5132*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*( 0.116- 0.0052095*x(7)+ 0.0233037*x(9)))-1)+ 23.3037;
  9. 101.779*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.752- 0.0100677*x(7)- 0.1118467*x(8)))-1)- 111.8467;
  10. - 111.8467*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*(- 0.502- 0.0100677*x(7)+ 0.101779*x(9)))-1)+ 101.779;
  11. 111.461*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.869- 0.0229274*x(7)- 0.1343884*x(8)))-1)- 134.3884;
  12. - 134.3884*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*( 0.166- 0.0229274*x(7)+ 0.111461*x(9)))-1)+ 111.461;
  13. 191.267*x(2)+(1-x(1)*x(2))*x(3)*(exp(x(5)*( 0.982- 0.0202153*x(7)-0.2114823*x(8)))-1)-211.4823;
  14. - 211.4823*x(1)+(1-x(1)*x(2))*x(4)*(exp(x(6)*(- 0.473- 0.0202153*x(7)+ 0.191267*x(9)))-1)+ 191.267;
  15. x(1)*x(3)-x(2)*x(4)];
  16. end
  17. end
  18. 用while sum(x>=0)~=9
  19. [x,x0]=myfun3;
  20. end
  21. 15次找到该解,其中还得到其它的解。
复制代码
 楼主| 发表于 2007-7-19 13:58 | 显示全部楼层
运行上面的Matlab代码,一直在重新赋值、重新计算,也没有马上结束的意思。用的是Matlab6.5,是否有影响?

Optimizer appears to be converging to a point which is not a root.
Norm of relative change in X is less than max(options.TolX^2,eps) but
sum-of-squares of function values is greater than or equal to sqrt(options.TolFun)
Try again with a new starting guess.
Maximum number of function evaluations reached:
increase options.MaxFunEvals
Optimization terminated successfully:
First-order optimality is less than options.TolFun.
Maximum number of function evaluations reached:
increase options.MaxFunEvals
Maximum number of function evaluations reached:
increase options.MaxFunEvals
...
发表于 2007-7-19 14:19 | 显示全部楼层
奇怪!我用的是当前最新的版本:2007a(7.4)程序中求解的方程组是放在nested function 结构里的,而nested function 结构是7.0以后的版本才支持的,6.5居然不报错,比较疑惑。
此外循环中凡是出现“Optimization terminated successfully:
First-order optimality is less than options.TolFun”地方说明那次循环已经找到解了,只不过不满足解的分量全部大于0的条件。dingd试试初值x0=[
  
   0.700364053278940  -0.405912701790691   0.983501160015575   1.785846699149356

   0.846205556765350   1.564006336325064  -0.268554818225819  -1.559900753328203

  -0.742964895856200]时候的结果。这个是我退出循环时候的初值。

[ 本帖最后由 rocwoods 于 2007-7-19 14:22 编辑 ]

评分

1

查看全部评分

发表于 2007-7-19 15:16 | 显示全部楼层

  1. X0=[];
  2. X=[];
  3. x=-1;
  4. exitflag=0;
  5. while sum(x>=0)~=9 | exitflag~=1
  6. [x,x0,fval,exitflag]=myfun3;
  7. if exitflag==1
  8. X0=[X0;x0];
  9. X=[X;x];
  10. end
  11. end
复制代码

用以上代码寻找解,以及全部大于0的解。找到一组解就记录下来,同时直到找到全部分量大于0的解退出循环。我这里24次循环找到了4组(其实是3组,两组重复了)解。
下面给出了解以及相应的初值:
解X
X =
  Columns 1 through 6
    0.8985    0.9740   11.6486   10.7462    3.2508    6.7105
    0.8232   -0.5533    0.6719   -0.9997    8.8545   -2.7651
    0.8985    0.9740   11.6486   10.7462    3.2508    6.7105
    0.9000    0.4500    1.0000    2.0001    8.0000    7.9997
  Columns 7 through 9
   -8.7640    1.2513   -0.5251
    6.0466    0.9759   -1.7085
   -8.7640    1.2513   -0.5251
    5.0000    1.0000    2.0001

相应的初值X0
X0 =
  Columns 1 through 6
   -0.0511   -0.3876   -0.7803   -0.4013    0.3746   -0.9944
    0.2531    0.0271    1.3996   -1.7883    1.9266   -1.8581
   -0.1861    1.7561   -1.2493   -1.2420    0.2723    0.3607
   -1.4724    1.0453    0.2343    1.3138    1.4768    1.7574
  Columns 7 through 9
    0.2907   -0.6741    0.8691
   -0.7516   -0.8574   -1.0143
    1.1024   -0.9404    0.7201
   -1.0865    1.3680   -0.4536
发表于 2007-7-19 15:41 | 显示全部楼层
呵呵,学习永远要向你们这样,我要继续努力
发表于 2008-1-25 18:10 | 显示全部楼层

回复 2楼 的帖子

你好!!
这个问题麻烦给看看.也许是我的初始值给的老不合适吧,老是求不出解.我是用MATLAB解的.看能否用你的-1stOpt帮忙计算一下.见笑了,有-1stOpt这样一个专门的软件吗?****************,并麻烦推荐一本相关的介绍书籍,我也学习下.非常感谢!!:@)
要求的是一个频率方程,即下面矩阵B里面的X值.
clc;clear
syms s1 s2   x
s1=sqrt(-2*pi^2+sqrt(4*pi^4+x^2));
s2=sqrt(2*pi^2+sqrt(4*pi^4+x^2));
B=[0 1 0 1;s1 0 s2 0;sin(s1) cos(s1) sinh(s2) cosh(s2);s1*cos(s1) -s1*sin(s1) s2*cosh(s2) s2*sinh(s2)]
f2=det(sym(B))
ff=simple(f2)
h=optimset;h.Display='off';h.Tolx=1e-16;h.TolFun=1e-16;
f=inline('ff');
for x0=0:1:1000
[x,y,flag,out]=fsolve(f,x0,h)
end

____________________________

本版不讨论软件获得问题,所以帮你重新编辑了下

[ 本帖最后由 sigma665 于 2008-1-25 20:58 编辑 ]
 楼主| 发表于 2008-1-26 00:01 | 显示全部楼层
最好把问题写成标准的公式形式,好懂就容易回答。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 16:21 , Processed in 0.074117 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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