matlab在科学计算中的应用7 - 代数方程与最优化问题的求解

第七章代数方程与最优化问题的求解•代数方程的求解•无约束最优化问题的计算机求解•有约束最优化问题的计算机求解•整数规划问题的计算机求解7.1代数方程的求解7.1.1代数方程的图解法•一元方程的图解法例：ezplot('exp(-3*t)…*sin(4*t+2)+4*exp…(-0.5*t)*cos(2*t)-…0.5',[05])holdon,line([0,5],[0,0])%同时绘制横轴验证：symst;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)')%第一个方程曲线holdon%保护当前坐标系ezplot(‘x^2*…cos(x+y^2)+…y^2*exp(x+y)')•方程的图解法仅适用于一元、二元方程的求根问题。7.1.2多项式型方程的准解析解法•例：ezplot('x^2+y^2-1');holdon%绘制第一方程的曲线ezplot(‘0.75*x^3-y+0.9’)%绘制第二方程为关于x的6次多项式方程应有6对根。图解法只能显示求解方程的实根。•一般多项式方程的根可为实数，也可为复数。可用MATLAB符号工具箱中的solve()函数。最简调用格式：S=solve(eqn1,eqn2,…,eqnn)（返回一个结构题型变量S，如S.x表示方程的根。）直接得出根：（变量返回到MATLAB工作空间）[x,…]=solve(eqn1,eqn2,…,eqnn)同上，并指定变量[x,…]=solve(eqn1,eqn2,…,eqnn，’x,…’)•例：symsxy;[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]•验证[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的符号工具箱得出原始问题的高精度数值解，故称之为准解析解。•例：[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=271x(22),y(22),z(22)ans=-1.0869654762986136074917644096117ans=.37264668450644375527750811296627e-1ans=.89073290972562790151300874796949验证：err=[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];norm(double(eval(err)))eval可对字符型变量进行符号求解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•例：求解（含变量倒数）symsxy;[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=261err=[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•例：求解（带参数方程）symsabxy;[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*b-3*a^3)^(1/2))][1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))]y=[a+1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))+3][a+1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))+3]7.1.3一般非线性方程数值解非线性方程的标准形式为f(x)=0函数fzero格式x=fzero(fun,x0)%用fun定义表达式f(x)，x0为初始解。x=fzero(fun,x0,options)[x,fval]=fzero(…)%fval=f(x)[x,fval,exitflag]=fzero(…)[x,fval,exitflag,output]=fzero(…)说明该函数采用数值解求方程f(x)=0的根。例:求的根解：fun='x^3-2*x-5';z=fzero(fun,2)%初始估计值为2z=2.0946formatlongopt=optimset('Tolx',1.0e-8);y=fzero('x^3-2*x-5',2,opt)y=2.094551482117093250xx7.1.4一般非线性方程组数值解•格式：最简求解语句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)•例：自编函数：functionq=my2deq(p)q=[p(1)*p(1)+p(2)*p(2)-1;0.75*p(1)^3-p(2)+0.9];OPT=optimset;OPT.LargeScale='off';[x,Y,c,d]=fsolve('my2deq',[1;2],OPT)Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.x=0.35700.9341Y=1.0e-009*0.12150.0964c=1d=iterations:7funcCount:21algorithm:'trust-regiondogleg'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);%结果和上述完全一致，从略。Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.若改变初始值x0=[-1,0]T[x,Y,c,d]=fsolve(f,[-1,0]',OPT);x,Y,kk=d.funcCountOptimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.x=-0.98170.1904Y=1.0e-010*0.5619-0.4500kk=15初值改变有可能得出另外一组解。故初值的选择对解的影响很大，在某些初值下甚至无法搜索到方程的解。•例：用solve()函数求近似解析解symst;solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=.67374570500134756702960220427474％不允许手工选择初值，只能获得这样的一个解。可先用用图解法选取初值，再调用fsolve()函数数值计算formatlongy=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)NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius121.8634e-0095.16e-0051243.67694e-0193.61071e-0057.25e-0101Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.t=3.52026389294877f=-6.063776702980306e-010•重新设置相关的控制变量，提高精度。ff=optimset;ff.TolX=1e-16;ff.TolFun=1e-30;ff.Display='iter';[t,f]=fsolve(y,3.5203,ff)NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius12