第7章 matlab代数方程与最优化问题的求解

• 验证
>> [eval('x.^2+y.^2-1') eval('75*x.^3/100-y+9/10')] ans = [ 0, 0] [ 0, 0] [ 0, 0] [ -.1e-31, 0] [ .5e-30+.1e-30*i, 0] [ .5e-30-.1e-30*i, 0] 由于方程阶次可能太高，不存在解析解。然而， 可利用MATLAB的符号工具箱得出原始问题的高精 度数值解，故称之为准解析解。
c= 1 d= iterations: 7 funcCount: 21 algorithm: 'trust-region dogleg' firstorderopt: 1.3061e-010 ％解回代的精度 调用inline( )函数： >> f=inline('[p(1)*p(1)+p(2)*p(2)-1; 0.75*p(1)^3-p(2)+0.9]','p'); >> [x,Y] = fsolve(f,[1; 2],OPT); % 结果和上述完全一致，从略。 Optimization terminated successfully: First-order optimality is less than options.TolFun.
• 一般多项式方程的根可为实数，也可为复数。 可用MATLAB符号工具箱中的solve( )函数。 最简调用格式： S=solve(eqn1,eqn2,…,eqnn) （返回一个结构题型变量S，如S.x表示方程的根。）
直接得出根： （变量返回到MATLAB工作空间） [x,…]=solve(eqn1,eqn2,…,eqnn) 同上，并指定变量 [x,…]=solve(eqn1,eqn2,…,eqnn，’x,…’)
验证： >> syms t ; t=3.52028; vpa(exp(-3*t)*sin(4*t+2)+… 4*exp(-0.5*t)*cos(2*t)-0.5) ans = -.19256654148425145223200161126442e-4
• 二元方程的图解法 例：
>> ezplot('x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)') % 第一个方程曲线 >> hold on % 保护当前坐标系 >> ezplot(‘x^2 *… cos(x+y^2) +… y^2*exp(x+y)')
• 例：求解 （含变量倒数） >> syms x y; >> [x,y]=solve('x^2/2+x+3/2+2/y+5/(2*y^2)+3/x^3=0',... 'y/2+3/(2*x)+1/x^4+5*y^4','x,y'); >> size(x) ans = 26 1 >> err=[x.^2/2+x+3/2+2./y+5./(2*y.^2)+3./x.^3,y/2+3./ (2*x)+1./x.^4+5*y.^4]; ％验证 >> norm(double(eval(err))) ans = 8.9625e-030
可先用用图解法选取初值，再调用fsolve( )函数 数值计算
>> format long >> y=inline('exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5','t'); >> ff=optimset; ff.Display='iter'; [t,f]=fsolve(y,3.5203,ff) Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 1 2 1.8634e-009 5.16e-005 1 2 4 3.67694e-019 3.61071e-005 7.25e-010 1 Optimization terminated successfully: First-order optimality is less than options.TolFun. t= 3.52026389294877 f= -6.063776702980306e-010
• 例： >> syms x y; >> [x,y]=solve('x^2+y^2-1=0','75*x^3/100-y+9/10=0')
x= [ -.98170264842676789676449828873194] [ -.55395176056834560077984413882735-.35471976465080793456863789934944*i] [ -.55395176056834560077984413882735+.35471976465080793456863789934944*i] [ .35696997189122287798839037801365] [ .86631809883611811016789809418650-1.2153712664671427801318378544391*i] [ .86631809883611811016789809418650+1.2153712664671427801318378544391*i] y= [ .19042035099187730240977756415289] [ .92933830226674362852985276677202-.21143822185895923615623381762210*i] [ .92933830226674362852985276677202+.21143822185895923615623381762210*i] [ .93411585960628007548796029415446] [ -1.4916064075658223174787216959259-.70588200721402267753918827138837*i] [ -1.4916064075658223174787216959259+.70588200721402267753918827138837*i]
• 方程的图解法 仅适用于一元、 二元方程的求根 问题。
7.1.2 多项式型方程的准解析解法
• 例：
>> ezplot('x^2+y^2-1'); hold on % 绘制第一方程的曲线 >> ezplot(‘0.75*x^3-y+0.9’) % 绘制第二方程 为关于x的6次多项式方程 应有6对根。图解法只能 显示求解方程的实根。
若改变初始值 x0=[-1,0]T
>> [x,Y,c,d]=fsolve(f,[-1,0]',OPT); x, Y, kk=d.funcCount Optimization terminated successfully: First-order optimality is less than options.TolFun. x= -0.9817 0.1904 Y= 1.0e-010 * 0.5619 -0.4500 kk = 15
7.1.3 一般非线性方程数值解
• 格式： 最简求解语句 x=fsolve(Fun, x0) 一般求解语句 [x, f, flag, out]=fsolve(Fun, x0,opt, p1, p2,…) 若返回的flag 大于0，则表示求解成功，否则求 解出现问题， opt 求解控制参数，结构体数据。 获得默认的常用变量 opt=optimset; 用这两种方法修改参数（解误差控制量） opt.Tolx=1e-10; 或 set(opt.‘Tolx’, 1e-10)
初值改变有可能得出另外一组解。故初值的选择 对解的影响很大，在某些初值下甚至无法搜索到方程 的解。
• 例：
用solve( )函数求近似解析解
>> syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)* cos(2*t) - 0.5) ans = .67374570500134756702960220427474 ％不允许手工选择初值，只能获得这样的一个解。
第七章代数方程与最优化问题的求解
• 代数方程的求解 • 无约束最优化问题的计算机求解 • 有约束最优化问题的计算机求解 • 整数规划问题的计算机求解
7.1代数方程的求解
7.1.1 代数方程的图解法
• 一元方程的图解法 例：
>> ezplot('exp(-3*t)… *sin(4*t+2)+4*exp… (-0.5*t)*cos(2*t)-… 0.5',[0 5]) >> hold on, >> line([0,5],[0,0]) % 同时绘制横轴
• 重新设置相关的控制变量，提高精度。
>> ff=optimset; ff.TolX=1e-16; ff.TolFun=1e-30; >> ff.Display='iter'; [t,f]=fsolve(y,3.5203,ff)
First-order Trust-region Iteration Func-count f(x) optimality radius 1 2 1.8634e-009 5.16e-005 1 2 4 3.67694e-019 3.61071e-005 7.25e-010 1 3 6 0 5.07218e-010 0 1 Optimization terminated successfully: First-order optimality is less than options.TolFun. t= 3.52026389244155err=[x+3*y.^3+2*z.^2-1/2, x.^2+3*y+z.^3-2, x.^3+2*z+2*y.^22/4]; >> norm(double(eval(err))) ans = 1.4998e-027
• 多项式乘积形式也可，如把第三个方程替换一下。
>> [x,y,z]=solve('x+3*y^3+2*z^2=1/2','x^2+3*y+z^3=2','x^3+ 2*z*y^2=2/4'); >> err=[x+3*y.^3+2*z.^2-1/2, x.^2+3*y+z.^3-2, x.^3+2*z.*y.^2-2/4]; >> norm(double(eval(err))) % 将解代入求误差 ans = 5.4882e-028
• 例：
自编函数： function q = my2deq(p) q=[p(1)*p(1)+p(2)*p(2)-1; 0.75*p(1)^3-p(2)+0.9]; >> OPT=optimset; rgeScale='off'; >> [x,Y,c,d] = fsolve('my2deq',[1; 2],OPT) Optimization terminated successfully: First-order optimality is less than options.TolFun. x= 0.3570 0.9341 Y= 1.0e-009 * 0.1215 0.0964
• 例：求解 （带参数方程） >> syms a b x y; >> [x,y]=solve('x^2+a*x^2+6*b+3*y^2=0','y=a+(x+3)','x,y') x= [ 1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b3*a^3)^(1/2))] [ 1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b3*a^3)^(1/2))] y= [ a+1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b3*a^3)^(1/2))+3] [ a+1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b3*a^3)^(1/2))+3]
• 例：
>> [x,y,z]=solve('x+3*y^3+2*z^2=1/2', 'x^2+3*y+z^3 =2 ' ,'x^3+2*z+2*y^2=2/4') ; >> size(x) ans = 27 1 >> x(22),y(22),z(22) ans = -1.0869654762986136074917644096117 ans = .37264668450644375527750811296627e-1 ans = .89073290972562790151300874796949