多项式线性方程-准解析解方法
简单方程求解
- 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]
- 匿名函数
- f=@(x)[xxx;xxx];
- Inline函数 （不推荐）
- f=inline(‘[xxx;xxx]’,’x’)
- M-函数
-
求解
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=@(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=@myout; x0=[2 1]; x=fminunc(f,x0,ff)
>> problem.solver='fminunc'; problem.options=optimset; problem.objective=@(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=@(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=@(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];
