Matlab-代数方程与最优问题

Matlab

Posted by Kingtous on January 18, 2019

多项式线性方程-准解析解方法

简单方程求解

  • solve
  • vpasolve
1
2
x+y=35
2x+4y=94
1
syms x y; [x0 y0]=solve(x+y==35,2*x+4*y==94)
1
2
3
[x1 y1]=vpasolve(x^2+y^2-1==0,...
						75*x^3/100-y+9/10==0)
<<<<<<< HEAD
1
2
3
4
5
syms x y z;
F=[equation1,equation2,equation3];
[x0 y0 z0]=vpasolve(F,[x,y,z]),size(x0)
%检验
norm(subs(F,{x,y,z},{x0,y0,z0}))
1
2
3
4
5
%含有参数的方程,需要声明要求解的变量
syms a b x y
[x1 y1]=solve(x^2+a*x^2...,[x,y])
%表明要求解的是x,y
%更高次带有参数的方程没有解析解
  • 优点
    • 求出解析解或高精度数值解
    • 可以指定初值
    • 可以同时求出方程的实数解和复数解
  • 缺点
    • 适用于多项式型方程
    • 相比数值解方法速度慢

非线性方程-数值解方法

1
x=fsolve(fun,x0)
1
[x,f,flag,out]=fsolve(fun,x0,opt)
1
[x,f,flag,out]=fsolve(fun,x0,opt,p1,p2,...)

其中f是误差向量,flag为负值则表示求解有问题,out存储中间信息

  • 变换标准形式

    • y=f(x)=0
  • 描述等式
    • M-函数
      • 建立m文件
      • function y=myfun(x)
      • y=[xxx;xxx]
    • 匿名函数
    • Inline函数 (不推荐)
      • f=inline(‘[xxx;xxx]’,’x’)
  • 求解

    1
    2
    3
    4
    5
    
    OPT=optimset;%控制精度
    %opt.TolX=1e-10
    %或者 set(opt,'TolX',1e-10)
    %还有 MaxIter MaxFcnEvals
    OPT.LargeScale='off'; [x,Y,c,d]=fsolve(f,[1;2],OPT)
    
    • 初值选实数,那就是实数根
    • 初值选复数,则都有可能
  • 检验解的正确性

多解矩阵方程

  • 随机选择初值,调用fsolve求解
  • 若得出的解已经被记录,则放弃这个解
  • 如果一段时间没有得出新的解,则终止程序
  • 可以用ctrl+c随时终止程序

无约束最优化问题

  • x=fminunc(fun,x0)
  • x=fminsearch(fun,x0)

将x,y换为x1,x2

1
2
3
4
5
>> f[email protected](x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); x0=[0; 0]; [x,b,c,d]=fminsearch(f,x0),
>> [x,b,c,d]=fminunc(f,[0;.0])
>> [x,y]=meshgrid(-3:.1:3, -2:.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); contour(x,y,z,30); ff=optimset; ff.OutputFcn[email protected]myout; x0=[2 1]; x=fminunc(f,x0,ff)

>> problem.solver='fminunc'; problem.options=optimset; problem.objective[email protected](x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); problem.x0=[2; 1]; [x,b,c,d]=fminunc(problem)

全局最优解

1
2
3
4
5
6
7
>> ezsurf('20+x1^2+x2^2-10*(cos(pi*x1)+cos(pi*x2))') 
>> view(0,90), shading flat 
>> f[email protected](x)20+(x(1)/30-1)^2+(x(2)/20-1)^2-10*(cos(pi*(x(1)/30-1))+cos(pi*(x(2)/20-1))); F=[]; tic, for i=1:100, [x,f0]=fminunc_global(f,-100,100,2,50); F=[F,f0]; end, toc
>> [x,y]=meshgrid(0.5:0.01:1.5); z=100*(y.^2-x).^2+(1-x).^2;contour3(x,y,z,100), zlim([0,310])
>> f[email protected](x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;ff=optimset; ff.TolX=1e-10; ff.TolFun=1e-20; x=fminunc(f,[0;0],ff)
>> syms x1 x2; f=100*(x2-x1^2)^2+(1-x1)^2; J=jacobian(f,[x1,x2])
>> ff.GradObj='on'; x=fminunc(@c6fun3,[0;0],ff)
1
2
3
4
5
6
7
8
9
%[x,f_min]=fminunc_global(fun,a,b,n,N)
% a b 搜索区间ab 决策变量n 尝试次数N
function [x,f0]=finunc_global(f,a,b,n,N,varargin)
k0=0;f0=Inf; if strcmp(class(f),'struct'),k0=1;end
for i=1:N ,x0=a+(b-a)*rand(n,1)
	if k0==1, f.x0=x0; [x1 f1 key]=fminunc(f)
    else, [x1 f1 key]=fminunc(f,x0,varargin{:});end
    if key>0 & f1<f0, x=x1;f0=f1;end
end
1
2
3
4
function [y,Gy]=c6fun3(x)
y=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
Gy=[-400*(x(2)-x(1)^2)*x(1)-2+2*x(1);
         200*x(2)-200*x(1)^2];

有约束最优化问题