Matlab讨论区 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

dingd 发表于 2007-7-12 15:22

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

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=;
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 发表于 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求解的代码：function p=myfun1
p0 = ;
options=optimset('MaxFunEvals',20000,'MaxIter',2000);
= fsolve(@f,p0,options) ;

function F=f(p)
F=[p(1)+p(2)*(1-exp(-(p(3)*(0)^p(4))))-51.61;
p(1)+p(2)*(1-exp(-(p(3)*(9.78)^p(4))))-51.91;
p(1)+p(2)*(1-exp(-(p(3)*(30.68)^p(4))))-53.27;
p(1)+p(2)*(1-exp(-(p(3)*(59.7)^p(4))))-59.68;];
end
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 发表于 2007-7-18 15:31

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

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

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

rocwoods 发表于 2007-7-18 16:04

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

function p=myfun1
p0 = ;
options=optimset('MaxFunEvals',20000,'MaxIter',5000);
= fsolve(@f,p0,options) ;

function F=f(p)
F=[p(1)+p(2)*(1-exp(-(p(3)*(0)^p(4))))-51.61;
p(1)+p(2)*(1-exp(-(p(3)*(9.78)^p(4))))-51.91;
p(1)+p(2)*(1-exp(-(p(3)*(30.68)^p(4))))-53.27;
p(1)+p(2)*(1-exp(-(p(3)*(59.7)^p(4))))-59.68;];
end
end

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

dingd 发表于 2007-7-18 16:18

赞一个:victory:

xjzuo 发表于 2007-7-18 16:24

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

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

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

dingd 发表于 2007-7-19 11:24

rocwoods的Matlab代码运行一下，大概能以1/10的概率求得正解。代码中如果初值都用随机值，即把 “p0 = ; ”改为“p0 = ; ”，没有人为控制因素，成功的概率就更低些了。看样子，与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 发表于 2007-7-19 11:57

用如下代码：

function =myfun3
x0 = unifrnd(0,1,1,9);
options=optimset('MaxFunEvals',20000,'MaxIter',5000);
= fsolve(@f,x0,options) ;

function F=f(x)
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;
- 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;
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;
- 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;
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;
- 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;
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;
- 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;
x(1)*x(3)-x(2)*x(4)];
end
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.97396780200530211.64863306079609210.746159415319836
Columns 5 through 8
3.250825125902567 6.710547896595700-8.764002776048228 1.251334798196427
Column 9
-0.525131249969468

rocwoods 发表于 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

function =myfun3
x0 = unifrnd(-2,2,1,9);
options=optimset('MaxFunEvals',10000,'MaxIter',5000);
= fsolve(@f,x0,options) ;

function F=f(x)
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;
- 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;
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;
- 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;
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;
- 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;
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;
- 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;
x(1)*x(3)-x(2)*x(4)];
end
end
用while sum(x>=0)~=9
=myfun3;
end
15次找到该解，其中还得到其它的解。

dingd 发表于 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
...

rocwoods 发表于 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 编辑 ]

rocwoods 发表于 2007-7-19 15:16

X0=[];
X=[];
x=-1;
exitflag=0;
while sum(x>=0)~=9 | exitflag~=1
=myfun3;
if exitflag==1
X0=;
X=;
end
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 发表于 2007-7-19 15:41

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

mechanic05 发表于 2008-1-25 18:10

回复 2楼的帖子

你好！！
这个问题麻烦给看看．也许是我的初始值给的老不合适吧，老是求不出解．我是用ＭＡＴＬＡＢ解的．看能否用你的-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=
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
=fsolve(f,x0,h)
end

____________________________

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

[ 本帖最后由 sigma665 于 2008-1-25 20:58 编辑 ]

dingd 发表于 2008-1-26 00:01

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

页: [1] 2

声振论坛's Archiver

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

回复 2楼 的帖子

回复 2楼的帖子