声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 21076|回复: 5

[综合讨论] 关于 solve 和 fsolve 的简单比较

[复制链接]
发表于 2009-6-8 19:33 | 显示全部楼层 |阅读模式

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

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

x
以matlab R2008a版本为例,各版本出错提示可能有所不同。有不对之处,欢迎指正。
1.solve 和 fsolve的基本含义
matlab给出的关于solve 和 fsolve的基本描述为:
solve——Symbolic solution of algebraic equations
fsolve——Solve system of nonlinear equations
可见solve用于解决代数方程(组)的符号(解析)解,而fsolve用来解决非线性方程(组)的数值解。
【在matlab里面solve命令主要是用来求解代数方程(即多项式)的解,但是也不是说其它方程一个也不能解,不过求解非代数方程的能力相当有限,通常只能给出很特殊的实数解。从计算机的编程实现角度讲,如今的任何算法都无法准确的给出任意非代数方程的所有解,但是我们有很多成熟的算法来实现求解在某点附近的解。matlab也不例外,它也只能给出任意非代数方程在某点附近的解,函数有两个:fzero和fsolve,具体用法请用help或doc命令查询吧。如果还是不行,你还可以将问题转化为非线性最优化问题,求解非线性最优化问题的最优解,可以用的命令有:fminbnd, fminsearch, fmincon等等。】(引自:http://blog.sina.com.cn/s/blog_4c4af5c101008w9f.html,作者:ggbondg)
下面举几个例子:
  1. 例1:>> solve('a*x-1')
  2. ans =
  3. 1/a
  4. 例2:>> solve('exp(x)+sin(x)-2')
  5. ans =
  6. .44867191635127271149118657202662
复制代码

注:对于solve结果的显示,有时看起来比较长,可用vpa进行精度控制,如:
>> vpa(solve('exp(x)+sin(x)-2'),3)
ans =
.449
  1. 例3:>> fsolve(@(x)exp(x)+sin(x)-2,0)
  2. Optimization terminated: first-order optimality is less than options.TolFun.
  3. ans =
  4.     0.4487
复制代码

2.关于solve 和 fsolve求解方程组时的书写规则
对于solve,方程可以直接书写,不需要运算符”.”;
对于fsolve,当未知量与未知量有乘除操作或未知量有开方、幂等操作时运算符”.”可写也可不写(记得好像必须写,试了试,发现不写也行)。下面举几个例子:
  1. 例4:>> solve('x+y.^2-1','x.^2-y-3')
  2. ??? Error using ==> solve at 77
  3. ' x+y.^2-1 ' is not a valid expression or equation.
  4. 例5:>> solve('x+y^2-1','x^2-y-3')
  5. ans =
  6.     x: [4x1 sym]
  7.     y: [4x1 sym]
复制代码
  1. 例6:function shiyan
  2. clc
  3. clear
  4. x0=[0,0];
  5. fsolve(@mf,x0)

  6. function F=mf(x)
  7. F=[x(1)+x(2)^2-1;
  8.    x(1)^2-x(2)-3];
  9. %%%%%Result%%%%%%
  10. Optimizer appears to be converging to a minimum that is not a root:
  11. Sum of squares of the function values is > sqrt(options.TolFun).
  12. Try again with a new starting point.
  13. ans =
  14.     1.6268   -0.1537
复制代码


例7:把例6中的mf函数,换成如下再试试。
function F=mf(x)

F=[x(1)+x(2).^2-1;
x(1).^2-x(2)-3];

例8:把例6的初值x0设为x0=[-2,2]; 运行结果为:

Optimization terminated: first-order optimality is less than options.TolFun.
ans =
-2.1875
1.7854
可见,用fsolve解非线性方程组,比较依赖处置的选择,因此建议用fsolve解方程时,能大致了解问题的求解区间,以便选择合适的初值。
3. 关于solve 和 fsolve求解时,参数为多数值的求解
问题来源:http://forum.vibunion.com/thread-82658-1-1.html
类似问题描述:k=(5.0e+4):(1e+3):(6e+4);h=1.6e-6;n1=2.2899;n0=1.5040;n2=1.000;解方程组:
p1=sqrt(k.^2.*n1.^2-b.^2;p2=sqrt(b.^2-k.^2.*n2.^2);p0=sqrt(b.^2-k.^2.*n0.^2);
p1*h-pi-atan(p0./p1)-atan(p2./p1)=0。(见:http://forum.vibunion.com/thread-83102-1-1.html)
解决方法:这是非线性方程组,不过我们可以先用solve试试:

  1. 例9:问题如前所述。考虑用solve是否能解。
  2. clear
  3. clc
  4. k=(5.0e+4):(1e+3):(6e+4);
  5. h=1.6e-6;
  6. n1=2.2899;
  7. n0=1.5040;
  8. n2=1.000;
  9. for i=1:length(k)
  10. y=solve(['p1=sqrt(',num2str(k(i)),'^2*',num2str(n1),'^2-b^2)'],...
  11.     ['p2=sqrt(b^2-',num2str(k(i)),'^2*',num2str(n2)','^2)'],...
  12.     ['p0=sqrt(b^2-',num2str(k(i)),'^2.*',num2str(n0),'^2)'],...
  13.     ['p1*',num2str(h),'-',num2str(pi),'-atan(p0/p1)-atan(p2/p1)=0']);
  14. end
复制代码

%%%%%Result%%%%%
> In solve at 140

In shiyan at 9
y=solve(['p1=sqrt(',num2str(k(i)),'^2*',num2str(n1),'^2-b^2)'],...

['p2=sqrt(b^2-',num2str(k(i)),'^2*',num2str(n2)','^2)'],...

['p0=sqrt(b^2-',num2str(k(i)),'^2.*',num2str(n0),'^2)'],...

['p1*',num2str(h),'-',num2str(pi),'-atan(p0/p1)-atan(p2/p1)=0']);
Warning: Warning, solutions may have been lost
Warning: Explicit solution could not be found.
>> whos y

Name  Size  Bytes  Class  Attributes
y     0x0    64     sym
由此说明利用solve并不能解决这个复杂的非线性方程组,考虑数值解法。
  1. 例10:问题如例9,考虑fsolve求解问题。
  2. % 主程序
  3. clear
  4. clc
  5. global k h n1 n2 n0
  6. k1=(5.0e+4):(1e+3):(6e+4);
  7. h=1.6e-6;
  8. n1=2.2899;
  9. n0=1.5040;
  10. n2=1.000;
  11. x0=[0 1 0 0];
  12. b=zeros(1,length(k1));
  13. for i=1:length(k1)
  14.     k=k1(i);
  15.     y=fsolve(@myfun,x0);
  16.     b(i)=y(4);
  17. end
  18. plot(k1,b);
  19. % 子程序
  20. function F=myfun(x)
  21. global k h n1 n2 n0
  22. % p0 p1 p2 b----x(1) x(2) x(3) x(4)
  23. F=[sqrt(k^2*n1^2-x(4)^2)-x(2);
  24.    sqrt(x(4)^2-k^2*n2^2)-x(3);
  25.    sqrt(x(4)^2-k^2*n0^2)-x(1);
  26.    x(2)*h-pi-atan(x(1)./x(2))-atan(x(3)./x(2))];
  27. %%%%% Result%%%%%%%%%
  28. y=fsolve(@myfun,x0);
  29. Optimizer appears to be converging to a minimum that is not a root:
  30. Sum of squares of the function values is > sqrt(options.TolFun).
  31. Try again with a new starting point.
复制代码


这说明所选初值并不合理。
  1. 例11 例10还可以这样改写:
  2. % 主程序
  3. clear
  4. clc
  5. % global k h n1 n2 n0
  6. k1=(5.0e+4):(1e+3):(6e+4);
  7. h=1.6e-6;
  8. n1=2.2899;
  9. n0=1.5040;
  10. n2=1.000;
  11. x0=[0 1 0 0];
  12. b=zeros(1,length(k1));
  13. for i=1:length(k1)
  14.     k=k1(i);
  15.     y=fsolve(@(x)myfun(x,k, h, n1, n2, n0),x0);
  16.     b(i)=y(4);
  17. end
  18. plot(k1,b);
  19. % 子程序
  20. function F=myfun(x,k, h, n1, n2, n0)
  21. % global k h n1 n2 n0
  22. % p0 p1 p2 b----x(1) x(2) x(3) x(4)
  23. F=[sqrt(k^2*n1^2-x(4)^2)-x(2);
  24.    sqrt(x(4)^2-k^2*n2^2)-x(3);
  25.    sqrt(x(4)^2-k^2*n0^2)-x(1);
  26.    x(2)*h-pi-atan(x(1)./x(2))-atan(x(3)./x(2))];
复制代码


例12 参见http://forum.vibunion.com/thread-82658-1-1.html此贴,看看大侠ChaChing,dingd(1stOpt解法)等的解法。
注:利用fsolve解数值解,初值的选择十分重要。而1stOpt则对初值的选择要求比较低,不妨一试。关于这一点,请参考如下文章:
作者:dingd 非线性方程组-1stOpt与fsolve的比较:
http://forum.vibunion.com/thread-48384-1-5.html

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2009-9-26 17:13 | 显示全部楼层
fsolve函数 求解非线性最小二乘解 是不是最好的方法??求解答
 楼主| 发表于 2009-9-26 22:43 | 显示全部楼层

回复 沙发 cammer534 的帖子

这个不一定。fsolve求解需要指定初始值。解的好坏应该和函数采用的算法有很大关系。

评分

1

查看全部评分

发表于 2011-3-30 23:23 | 显示全部楼层
讲的非常好,谢谢。
发表于 2011-4-1 09:57 | 显示全部楼层
什么软件都不是万能的,有时候也不能太相信!
发表于 2011-4-14 14:21 | 显示全部楼层
感谢楼主分享。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-29 10:06 , Processed in 0.058732 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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