非线性方程组-1stOpt与fsolve的比较！

dingd

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 编辑 ]

rocwoods

虽然找合适的初值是个问题，但是楼主也不要把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 编辑 ]

dingd

初值猜的很精准啊，佩服！不过对于不少人而言，“功力”是不容易达到这个水准的。

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

其实还有一组解，再试下。（如果是参照1stOpt的结果来定初值，就没必要试了）

rocwoods

p =
  51.609999999999999
  -0.134520891677992
  -0.241198737477775
   0.693446385900985
第二组值，说实在的，不知道是我运气好还是怎么着，还是受1stOpt初值随机给的启发。我在fsolve中也随机给初值。代码如下：

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

复制代码

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

dingd

赞一个:victory:

xjzuo

我也还是持一贯的观点: 充分利用各个软件的优点才是最重要的. 不宜贬低某个软件,当然前提是都比较强.
拟合优化是1stOpt的强项,当然我们就应该利用它；但就象我们不能因为它绘图、矩阵计算不强，无法解微分方程等等，就把1stOpt贬低得一无是处一样。

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

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

dingd

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;

rocwoods

用如下代码：

1. 
   
2. function [x,x0]=myfun3
   
3.    x0 = unifrnd(0,1,1,9);
   
4.    options=optimset('MaxFunEvals',20000,'MaxIter',5000);
   
5.     [x,fval] = fsolve(@f,x0,options) ;
   
6.    
   
7. function F=f(x)
   
8. 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;
   
9. - 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;
   
10. 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;
    
11. - 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;
    
12. 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;
    
13. - 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;
    
14. 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;
    
15. - 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;
    
16. x(1)*x(3)-x(2)*x(4)];
    
17. end
    
18. 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

rocwoods

全部大于等于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. 
   
2. function [x,x0]=myfun3
   
3.    x0 = unifrnd(-2,2,1,9);
   
4.    options=optimset('MaxFunEvals',10000,'MaxIter',5000);
   
5.     [x,fval] = fsolve(@f,x0,options) ;
   
6.    
   
7. function F=f(x)
   
8. 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;
   
9. - 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;
   
10. 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;
    
11. - 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;
    
12. 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;
    
13. - 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;
    
14. 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;
    
15. - 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;
    
16. x(1)*x(3)-x(2)*x(4)];
    
17. end
    
18. end
    
19. 用while sum(x>=0)~=9
    
20. [x,x0]=myfun3;
    
21. end
    
22. 15次找到该解，其中还得到其它的解。
    

复制代码

dingd

运行上面的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
...

rocwoods

奇怪！我用的是当前最新的版本：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 编辑 ]

rocwoods

1. 
   
2. X0=[];
   
3. X=[];
   
4. x=-1;
   
5. exitflag=0;
   
6. while sum(x>=0)~=9 | exitflag~=1
   
7. [x,x0,fval,exitflag]=myfun3;
   
8. if exitflag==1
   
9. X0=[X0;x0];
   
10. X=[X;x];
    
11. end
    
12. 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

caizi2008

呵呵,学习永远要向你们这样,我要继续努力

mechanic05

你好！！
这个问题麻烦给看看．也许是我的初始值给的老不合适吧，老是求不出解．我是用ＭＡＴＬＡＢ解的．看能否用你的-1stOpt帮忙计算一下．见笑了，有-1stOpt这样一个专门的软件吗？****************，并麻烦推荐一本相关的介绍书籍，我也学习下．非常感谢！！:@)
要求的是一个频率方程，即下面矩阵Ｂ里面的Ｘ值．
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 编辑 ]

dingd

最好把问题写成标准的公式形式，好懂就容易回答。