K =
-1.5954 1.2040 -0.4201 -2.9046 fval = -5.0000 -3.0000 -1.0000 attainfactor = 8.4703e-022
本例中,优化器试图使对象值与目标值匹配,达到因子(取为8.4703e-022)显示目标与对象的匹配程度很高。
10.4 非线性方程的优化解
10.4.1 非线性方程的求解
非线性方程的求解问题可以看作是单变量的最小化问题,通过不断缩小搜索区间来逼近问题解的真值。单变量最小化问题的求解方法请参见前面的内容。在Matlab中,非线性方程求解所采用的算法是二分法、secant法和逆二次内插法的组合。
函数:fzero
功能:求单变量函数的零点。 格式:x = fzero(fun,x0)
x = fzero(fun,x0,options)
x = fzero(fun,x0,options,P1,P2,...) [x,fval] = fzero(...)
[x,fval,exitflag] = fzero(...)
[x,fval,exitflag,output] = fzero(...) 说明:
x = fzero(fun,x0) 如果x0为标量,函数试图找到x0 附近fun函数的零点。fzero函数返回的x值为fun函数改变符号处邻域内的点,或者是NaN(如果搜索失败)。这里,当函数发现Inf、NaN或复数时,搜索终止。若x0为一长度为2的向量,fzero函数假设x0为一区间,其中fun (x0 (1))的符号与fun (x0 (2))的符号相反。当该情况不为真时,发生错误。用此区间调用fzero函数可以保证fzero函数返回fun函数改变符号处附近的点。
x = fzero(fun,x0,options) 用options结构指定的参数进行最小化。
x = fzero(fun,x0,options,P1,P2,...) 提供其他变量如P1,P2等并传递给目标函数fun。如果没有选项设置,令options = [ ]。
[x,fval] = fzero(...) 返回解x处目标函数的值。
[x,fval,exitflag] = fzero(...) 返回exitflag参数,描述退出条件。
[x,fval,exitflag,output] = fzero(...) 返回包含优化信息的输出结构output。 各变量意义同前。 注意:
(1)调用fzero函数时,使用初值区间(二元素的x0)常常比用标量x0快;
(2)fzero命令给零点的定义是函数与x轴相交的点。函数与x轴接触但并没有穿过x轴的点不算作有效零点。如y = x^2函数曲线便在0处与x轴接触,但没有穿过x轴,所以没有发现零点。对于没有有效零点的函数,fzero函数将一直运行到发现Inf、NaN和复数值。
例10-18 通过计算sin函数在3附近的零点来计算?。 解:在Matlab中实现: >> x=fzero(@sin,3) x =
3.1416
例10-19 找到cos函数在1和2之间的零点。
解:在Matlab中实现: >> x=fzero(@cos,[1,2]) x =
1.5708
注意cos (1)和cos (2)的符号相反。 例10-20 找到下面函数的零点:
f(x)?x3?2x?5 解:寻找2附近的零点: >> x=fzero('x^3-2*x-5',2) x =
2.0946
因为此函数为一多项式,也可用: >> roots([1 0 -2 -5]) ans =
2.0946 -1.0473 + 1.1359i -1.0473 - 1.1359i
找到相同的实数零点和共轭复数解。
10.4.2 非线性方程组的求解
非线性方程组的数学模型为:
F (x) = 0 其中x为一向量,F (x)为一函数,返回向量值。
在Matlab中,用fsolve函数求解非线性方程组。其算法基于最小二乘法。 函数:fsolve
功能:求解非线性等式系统。 格式:x = fsolve(fun,x0)
x = fsolve(fun,x0,options)
x = fsolve(fun,x0,options,P1,P2, ... ) [x,fval] = fsolve(...)
[x,fval,exitflag] = fsolve(...)
[x,fval,exitflag,output] = fsolve(...)
[x,fval,exitflag,output,jacobian] = fsolve(...) 说明:
x = fsolve(fun,x0) 初值为x0,试图求解由fun函数描述的等式系统。 x = fsolve(fun,x0,options) 用options结构指定的参数进行最小化。
x = fsolve(fun,x0,options,P1,P2, ... ) 提供其他变量如P1,P2等并传递给目标函数fun。如果没有选项设置,令options = [ ]。
[x,fval] = fsolve(...) 返回解x处目标函数的值。
[x,fval,exitflag] = fsolve(...) 返回exitflag参数,描述退出条件。
[x,fval,exitflag,output] = fsolve(...) 返回包含优化信息的输出结构output。 [x,fval,exitflag,output,jacobian] = fsolve(...) 返回解x处的fun函数的Jacobian矩阵。
·fun变量
为目标函数,即需要最小化的目标函数。Fun函数需要输入向量参数x,返回x处的目标函数向量值F。可以将fun函数指定为命令行,如
x = fsolve(inline('sin(x.*x)'),x0)
同样,fun参数也可以是一个包含函数名的字符串。对应的函数可以是M文件、内部函数和MEX文件。若fun = ‘myfun’,则
x = fsolve(@myfun,x0) 其中M文件myfun .m必须是下面的形式:
function F = myfun(x)
F = ... % 计算x处的目标函数
若Jacobian矩阵可以算得,并且options. Jacobian设置为’on’,即 options = optimset (‘Jacobian’, ‘on’)
则fun函数必须在第二个输出变量中输出x处的Jacobian矩阵J。注意,如果fun函数只要求有一个输出变量(此时优化算法只要求目标函数值而不要求Jacobian值),则通过核对nargout参数可以避免计算Jacobian矩阵。
function [F, J] = myfun(x)
F = ... % x处的目标函数值 if nargout > 1 % 两个输出变量 J = ... % x处的Jacobian函数 end
若fun函数返回一个有m个元素的向量或矩阵,长度为n的x,则Jacobian矩阵J为m*n的矩阵,其中J (i, j)为F (i)对x (j)的偏导数。
其它参数意义同前。 注意:
(1)要求解的函数必须是连续的。当成功收敛时,fsolve函数只给出一个根。当函数收敛于非零点时,给出下面的信息:
Optimizer is stuck at a minimum that is not a root Try again with a new starting guess
此时,重新给定初值并运行fsolve函数进行试算。
(2)fsolve函数只对实数变量有效。当x为复数时,必须将它分为实数部分和虚数部分。
(3)大型化问题 现在,若在fun函数中提供了Jacobian矩阵,options参数不能与大型算法同时使用以比较解析Jacobian矩阵和有限差分Jacobian矩阵。可以通过将参数MaxIter设置为0来核对导数。然后用大型算法重新运行程序。
例10-21 求解下列方程组
2x1?x2?e?x1 ?x1?2x2?e?x2 解:我们试图求解下面的系统:
2x1?x2?e?x1?0
?x1?2x2?e?x2?0 初值为x0 = [-5, -5]
首先我们编写一个M文件Ex1021.m,计算x处的等式的值F: function F = Ex1021(x)
F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))]; 调用优化函数:
>> x0=[-5;-5]; %初值
>> options=optimset('Display','iter'); %输出显示的选项 >> [x,fval]=fsolve(@Ex1021,x0,options) %调用优化函数
Norm of First-order
Iteration Func-count f(x) step optimality CG-iterations
1 4 47071.2 1 2.29e+004 0 2 7 6527.47 1.45207 3.09e+003 1 3 10 918.372 1.49186 418 1 4 13 127.74 1.55326 57.3 1 5 16 14.9153 1.57591 8.26 1
6 19 0.779051 1.27662 1.14 1 7 22 0.00372453 0.484658 0.0683 1 8 25 9.21615e-008 0.0385552 0.000336 1 9 28 5.65915e-017 0.000193707 8.34e-009 1 表中各列的意义分别为:
Iteration 迭代次数;
Func-count 目标函数计算次数; f (x) 目标函数值; Norm of step 当前步长的范数;
First-order optimiltity 当前梯度的无限范数;
CG-Iterations 当前迭代中源于PCG的迭代次数。 Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun x =
0.5671 0.5671 fval =
1.0e-008 * -0.5319 -0.5319
说明优化过程成功终止,函数值的相对改变小于Options.TolFun。x1和x2的取值均为0.5671,两个函数的目标值均为1.0e-008 *-0.5319。
例10-22 求矩阵x,使其满足方程
?12? X*X*X???
34??初值为x = [1, 1; 1, 1]。
解:首先编写一待求等式的M文件Ex1022.m: function F = Ex1022(x) F = x*x*x-[1,2;3,4]; 然后调用优化过程:
>> x0=ones(2,2); %设初值
>> options=optimset('Display','off'); %取消显示 >> [x,Fval,exitflag]=fsolve(@Ex1022,x0,options) x =
-0.1291 0.8602 1.2903 1.1612 Fval =
1.0e-003 *
0.1541 -0.1163 0.0109 -0.0243 exitflag = 1
残差接近于0:
>> sum(sum(Fval.*Fval)) ans =
3.7973e-008
例10-23 求解方程 cosx+x = 0 解:在Matlab中实现: >> fsolve('cos(x)+x',0)
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun ans =
-0.7391
>> cos(ans) %检验解 ans =
0.7391