第六章 代数方程的解
? 代数方程的图解法
? 多项式型方程的准解析解法
? 一般非线性方程数值解
6.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])
% 同时绘制横轴
验证,
>> syms x 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('y^2 *…
cos(y+x^2) +…
x^2*exp(x+y)')
? 方程的图解法
仅适用于一元,
二元方程的求根
问题。
6.1.2 多项式型方程的准解析解法
? 例,
>> ezplot('x^2+y^2-1'); hold on % 绘制第一方程的曲线
>> ezplot(‘0.75*x^3-y+0.9’) % 绘制第二方程
为关于 x的 6次多项式方程
应有 6对根。图解法只能
显示求解方程的实根。
? 一般多项式方程的根可为实数，也可为复数。
可用 MATLAB符号工具箱中的 solve( )函数。
最简调用格式,
S=solve(eqn1,eqn2,…,eqn n)
（ 返回一个结构题型变量 S，如 S.x表示方程的根。）
直接得出根,（变量返回到 MATLAB工作空间）
[x,…]=solve(eqn 1,eqn2,…,eqn n)
同上，并指定变量
[x,…]=solve(eqn 1,eqn2,…,eqn n，’ x,…’)
? 例,
>> 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]
? 验证
>> [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') ;
>> x(22),y(22),z(22)
ans =
-1.0869654762986136074917644096117
ans =
,37264668450644375527750811296627e-1
ans =
,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)))
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
? 例：求解
（含变量倒数）
>> 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
? 例：求解
（带参数方程）
>> 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*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]
6.1.3 一般非线性方程数值解
? 格式,
最简求解语句
x=fsolve(Fun,x0)
一般求解语句
[x,f,flag,out]=fsolve(Fun,x0,opt,p1,p2,…)
若返回的 flag 大于 0，则表示求解成功，否则求
解出现问题。
获得默认的常用变量
opt=optimset;
如：用这两种方法修改参数（解误差控制量）
opt.Tolx=1e-10; 或 set(opt,‘Tolx’,1e-10)
? 例,
自编函数,
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; OPT.LargeScale='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
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,
若改变初始值 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
初值改变有可能得出另外一组解。故初值的选择
对解的影响很大，在某些初值下甚至无法搜索到方程
的解。
? 例,
用 solve( )函数求近似解析解
>> syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*
cos(2*t) - 0.5)
ans =
.67374570500134756702960220427474
％不允许手工选择初值，只能获得这样的一个解。
可先用用图解法选取初值，再调用 f solve( )函数
数值计算
>> 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
? 重新设置相关的控制变量，提高精度。
>> ff=optimset; ff.TolX=1e-16; ff.TolFun=1e-30;
>> 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
3 6 0 5.07218e-010 0 1
Optimization terminated successfully,
First-order optimality is less than options.TolFun,
t =
3.52026389244155
f =
0
? 代数方程的图解法
? 多项式型方程的准解析解法
? 一般非线性方程数值解
6.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])
% 同时绘制横轴
验证,
>> syms x 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('y^2 *…
cos(y+x^2) +…
x^2*exp(x+y)')
? 方程的图解法
仅适用于一元,
二元方程的求根
问题。
6.1.2 多项式型方程的准解析解法
? 例,
>> ezplot('x^2+y^2-1'); hold on % 绘制第一方程的曲线
>> ezplot(‘0.75*x^3-y+0.9’) % 绘制第二方程
为关于 x的 6次多项式方程
应有 6对根。图解法只能
显示求解方程的实根。
? 一般多项式方程的根可为实数，也可为复数。
可用 MATLAB符号工具箱中的 solve( )函数。
最简调用格式,
S=solve(eqn1,eqn2,…,eqn n)
（ 返回一个结构题型变量 S，如 S.x表示方程的根。）
直接得出根,（变量返回到 MATLAB工作空间）
[x,…]=solve(eqn 1,eqn2,…,eqn n)
同上，并指定变量
[x,…]=solve(eqn 1,eqn2,…,eqn n，’ x,…’)
? 例,
>> 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]
? 验证
>> [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') ;
>> x(22),y(22),z(22)
ans =
-1.0869654762986136074917644096117
ans =
,37264668450644375527750811296627e-1
ans =
,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)))
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
? 例：求解
（含变量倒数）
>> 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
? 例：求解
（带参数方程）
>> 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*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]
6.1.3 一般非线性方程数值解
? 格式,
最简求解语句
x=fsolve(Fun,x0)
一般求解语句
[x,f,flag,out]=fsolve(Fun,x0,opt,p1,p2,…)
若返回的 flag 大于 0，则表示求解成功，否则求
解出现问题。
获得默认的常用变量
opt=optimset;
如：用这两种方法修改参数（解误差控制量）
opt.Tolx=1e-10; 或 set(opt,‘Tolx’,1e-10)
? 例,
自编函数,
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; OPT.LargeScale='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
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,
若改变初始值 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
初值改变有可能得出另外一组解。故初值的选择
对解的影响很大，在某些初值下甚至无法搜索到方程
的解。
? 例,
用 solve( )函数求近似解析解
>> syms t; solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*
cos(2*t) - 0.5)
ans =
.67374570500134756702960220427474
％不允许手工选择初值，只能获得这样的一个解。
可先用用图解法选取初值，再调用 f solve( )函数
数值计算
>> 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
? 重新设置相关的控制变量，提高精度。
>> ff=optimset; ff.TolX=1e-16; ff.TolFun=1e-30;
>> 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
3 6 0 5.07218e-010 0 1
Optimization terminated successfully,
First-order optimality is less than options.TolFun,
t =
3.52026389244155
f =
0
