目录
递推关系式的作图程序
x(1)=1;x(2)=10;
hold on
for t=3:10
x(t)=0.8*x(t-1)+1.5*x(t-2); %递推关系式，注意MATLAB里向量的下标不能从零开始。
end
plot(2:10,x(2:10))顶点覆盖近似算法
%只是这种近似算法结果很差
%首先输入关联矩阵F及节点个数n
F=[0 1 0 0 0 0 0;
1 0 1 0 0 0 0;
0 1 0 1 1 1 0;
0 0 1 0 0 1 0;
0 0 1 0 0 1 0;
0 0 1 1 1 0 1;
0 0 0 0 0 1 0];
n=7;
C=[];
l=0;
for i=1:n
for j=1:n
if F(i,j)~=0
if l==0
C=[i j];l=2;
else
p=0;q=0;
for a=1:l
if C(a)==i
p=1;
end
if C(a)==j
q=1;
end
end
if p==0
l=l+1;C(l)=i;
end
if q==0
l=l+1;C(l)=j;
end
F(i,:)=zeros(1,n);
F(:,j)=zeros(n,1);
end
end
end
end
disp(C);
方程求根
% 方程求根
% inv - 逆矩阵
% roots - 多项式的根
% fzero - 一元函数零点
% fsolve - 非线性方程组
% solve - 符号方程解
% *newton - 牛顿迭代法解非线性方程INV.m
%求逆矩阵
%用法 B=inv(A) 其中A为数值或符号方阵,B返回A的逆
%例如
% inv([1 2;3 4]) %数值
% syms a b c d;inv([[a,b;c,d]) %符号
%
%INV Matrix inverse.
% INV(X) is the inverse of the square matrix X.
% A warning message is printed if X is badly scaled or
% nearly singular.
%
% See also SLASH, PINV, COND, CONDEST, NNLS, LSCOV.
% Copyright (c) 1984-98 by The MathWorks, Inc.
% $Revision: 5.5 $ $Date: 1997/11/21 23:38:31 $
% Built-in function.
ROOTS.m
function r = roots(c)
%roots(p)可求得多项式p的所有复根.
%例 x^3+2x^2-3=0
% 求解
% roots([1 2 0 -3])
%
%ROOTS Find polynomial roots.
% ROOTS(C) computes the roots of the polynomial whose coefficients
% are the elements of the vector C. If C has N+1 components,
% the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1).
%
% See also POLY, RESIDUE, FZERO.
% J.N. Little 3-17-86
% Copyright (c) 1984-98 by The MathWorks, Inc.
% $Revision: 5.6 $ $Date: 1997/11/21 23:41:05 $
% ROOTS finds the eigenvalues of the associated companion matrix.
if size(c,1)>1 & size(c,2)>1
error('Must be a vector.')
end
c = c(:).';
n = size(c,2);
r = zeros(0,1);
inz = find(c);
if isempty(inz),
% All elements are zero
return
end
% Strip leading zeros and throw away.
% Strip trailing zeros, but remember them as roots at zero.
nnz = length(inz);
c = c(inz(1):inz(nnz));
r = zeros(n-inz(nnz),1);
% Polynomial roots via a companion matrix
n = length(c);
if n > 1
a = diag(ones(1,n-2),-1);
a(1,:) = -c(2:n) ./ c(1);
r = [r;eig(a)];
end
FSOLVE.m
function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolve(FUN,x,options,varargin)
%优化工具箱函数,可求多元非线性方程组的实根. 用法与fzero类似。
%例 先写一个M函数rooteg4fun.m
% function y=rooteg4fun(x)
% y(1)=4*x(1)-x(2)+exp(x(1))/10-1;
% y(2)=-x(1)+4*x(2)+x(1).^2/8;
% 使用
% [x,f,h]=fsolve('rooteg4fun',[0,0]) %初值x(1)=0,x(2)=0
% x返回解向量，f返回误差向量，h>0表明算法收敛
% 注意：方程变量必须拼成一个向量变量，即用x(1),x(2),...
%
%FSOLVE Solves nonlinear equations by a least squares method.
%
% FSOLVE solves equations of the form:
%
% F(X)=0 where F and X may be vectors or matrices.
%
% X=FSOLVE(FUN,X0) starts at the matrix X0 and tries to solve the
% equations described in FUN. FUN is usually an M-file which returns
% an evaluation of the equations for a particular value of X: F=FUN(X).
%
% X=FSOLVE(FUN,X0,OPTIONS) minimizes with the default optimization
% parameters replaced by values in the structure OPTIONS, an argument
% created with the OPTIMSET function. See OPTIMSET for details. Used
% options are Display, TolX, TolFun, DerivativeCheck, Diagnostics, Jacobian,
% JacobPattern, LineSearchType, LevenbergMarquardt, MaxFunEvals, MaxIter,
% DiffMinChange and DiffMaxChange, LargeScale, MaxPCGIter, PrecondBandWidth,
% TolPCG, TypicalX. Use the Jacobian option to specify that FUN may be called
% with two output arguments where the second, J, is the Jacobian matrix:
% [F,J] = feval(FUN,X). If FUN returns a vector (matrix) of m components when
% X has length n, then J is an m-by-n matrix where J(i,j) is the partial
% derivative of F(i) with respect to x(j). (Note that the Jacobian J is the
% transpose of the gradient of F.)
%
% X=FSOLVE(FUN,X0,OPTIONS,P1,P2,...) passes the problem-dependent
% parameters P1,P2,... directly to the function FUN: FUN(X,P1,P2,...).
% Pass an empty matrix for OPTIONS to use the default values.
%
% [X,FVAL]=FSOLVE(FUN,X0,...) returns the value of the objective function
% at X.
%
% [X,FVAL,EXITFLAG]=FSOLVE(FUN,X0,...) returns a string EXITFLAG that
% describes the exit condition of FSOLVE.
% If EXITFLAG is:
% > 0 then FSOLVE converged to a solution X.
% 0 then the maximum number of function evaluations was reached.
% < 0 then FSOLVE did not converge to a solution.
%
% [X,FVAL,EXITFLAG,OUTPUT]=FSOLVE(FUN,X0,...) returns a structure OUTPUT
% with the number of iterations taken in OUTPUT.iterations, the number of
% function evaluations in OUTPUT.funcCount, the algorithm used in OUTPUT.algorithm,
% the number of CG iterations (if used) in OUTPUT.cgiterations, and the first-order
% optimality (if used) in OUTPUT.firstorderopt.
%
% [X,FVAL,EXITFLAG,OUTPUT,JACOB]=FSOLVE(FUN,X0,...) returns the
% Jacobian of FUN at X.
% Copyright (c) 1990-98 by The MathWorks, Inc.
% $Revision: 1.26 $ $Date: 1998/10/22 19:28:31 $
% Andy Grace 7-9-90.
% Grandfathered FSOLVE call for Optimization Toolbox versions prior to 2.0:
% [X,OPTIONS]=FSOLVE(FUN,X0,OPTIONS,GRADFUN,P1,P2,...)
%
% ------------Initialization----------------
defaultopt = optimset('display','final','LargeScale','on', ...
'TolX',1e-6,'TolFun',1e-6,'DerivativeCheck','off',...
'Jacobian','off','MaxFunEvals','100*numberOfVariables',...
'Diagnostics','off',...
'DiffMaxChange',1e-1,'DiffMinChange',1e-8,...
'PrecondBandWidth',0,'TypicalX','ones(numberOfVariables,1)','MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
'TolPCG',0.1,'MaxIter',400,'JacobPattern',[], ...
'LineSearchType','quadcubic','LevenbergMarq','off');
% If just 'defaults' passed in, return the default options in X
if nargin==1 & nargout <= 1 & isequal(FUN,'defaults')
x = defaultopt;
return
end
if nargin < 2, error('FSOLVE requires two input arguments');end
if nargin < 3, options=[]; end
% These are added so that we can have the same code as in lsqnonlin which
% actually has upper and lower bounds.
LB = []; UB = [];
%[x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolve(FUNin,x,options,varargin)
% Note: don't send varargin in as a comma separated list!!
numargin = nargin; numargout = nargout;
[calltype, GRADFUN, varargin] = parse_call(FUN,options,numargin,numargout,varargin);
if isequal(calltype,'new') % fsolve version 2.*
xstart=x(:);
numberOfVariables=length(xstart);
large = 'large-scale';
medium = 'medium-scale';
l = []; u = [];
options = optimset(defaultopt,options);
switch optimget(options,'display')
case {'off','none'}
verbosity = 0;
case 'iter'
verbosity = 2;
case 'final'
verbosity = 1;
case 'testing'
verbosity = Inf;
otherwise
verbosity = 1;
end
diagnostics = isequal(optimget(options,'diagnostics','off'),'on');
gradflag = strcmp(optimget(options,'Jacobian'),'on');
line_search = strcmp(optimget(options,'largescale','off'),'off'); % 0 means trust-region, 1 means line-search
% Convert to inline function as needed
if ~isempty(FUN) % will detect empty string, empty matrix, empty cell array
[funfcn, msg] = fprefcnchk(FUN,'fsolve',length(varargin),gradflag);
else
errmsg = sprintf('%s\n%s', ...
'FUN must be a function name, valid string expression, or inline object;', ...
' or, FUN may be a cell array that contains these type of objects.');
error(errmsg)
end
x(:) = xstart;
switch funfcn{1}
case 'fun'
fuser = feval(funfcn{3},x,varargin{:});
f = fuser(:);
nfun=length(f);
JAC = zeros(nfun,numberOfVariables);
case 'fungrad'
[fuser,JAC] = feval(funfcn{3},x,varargin{:});
f = fuser(:);
nfun=length(f);
case 'fun_then_grad'
fuser = feval(funfcn{3},x,varargin{:});
f = fuser(:);
JAC = feval(funfcn{4},x,varargin{:});
nfun=length(f);
otherwise
error('Undefined calltype in FSOLVE');
end
% check size of JAC
[Jrows, Jcols]=size(JAC);
if Jrows~=nfun | Jcols ~=numberOfVariables
errstr = sprintf('%s\n%s%d%s%d\n',...
'User-defined Jacobian is not the correct size:',...
' the Jacobian matrix should be ',nfun,'-by-',numberOfVariables);
error(errstr);
end
YDATA = []; caller = 'fsolve';
% trustregion and enough equations (as many as variables)
if ~line_search & nfun >= numberOfVariables
OUTPUT.algorithm = large;
% trust region and not enough equations -- switch to line_search
elseif ~line_search & nfun < numberOfVariables
warnstr = sprintf('%s\n%s\n', ...
'Large-scale method requires at least as many equations as variables; ',...
' switching to line-search method instead.');
warning(warnstr);
OUTPUT.algorithm = medium;
% line search and no bounds
elseif line_search & isempty(l) & isempty(u)
OUTPUT.algorithm = medium;
% line search and bounds and enough equations, switch to trust region
elseif line_search & (~isempty(LB) | ~isempty(UB)) & nfun >= numberOfVariables
warnstr = sprintf('%s\n%s\n', ...
'Line-search method does not handle bound constraints; ',...
' switching to trust-region method instead.');
warning(warnstr);
OUTPUT.algorithm = large;
% can't handle this one:
elseif line_search & (~isempty(LB) | ~isempty(UB)) & nfun < numberOfVariables
errstr = sprintf('%s\n%s\n%s\n', ...
'Line-search method does not handle bound constraints ',...
' and trust-region method requires at least as many equations as variables; ',...
' aborting.');
error(errstr);
end
if diagnostics > 0
% Do diagnostics on information so far
constflag = 0; gradconstflag = 0; non_eq=0;non_ineq=0;lin_eq=0;lin_ineq=0;
confcn{1}=[];c=[];ceq=[];cGRAD=[];ceqGRAD=[];
hessflag = 0; HESS=[];
msg = diagnose('fsolve',OUTPUT,gradflag,hessflag,constflag,gradconstflag,...
line_search,options,xstart,non_eq,...
non_ineq,lin_eq,lin_ineq,LB,UB,funfcn,confcn,f,JAC,HESS,c,ceq,cGRAD,ceqGRAD);
end
% Execute algorithm
if isequal(OUTPUT.algorithm, large)
if ~gradflag
Jstr = optimget(options,'JacobPattern',[]);
if isempty(Jstr)
% Put this code separate as it might generate OUT OF MEMORY error
Jstr = sparse(ones(Jrows,Jcols));
end
else
Jstr = [];
end
l = []; u = []; computeLambda = 0;
[x,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT]=...
snls(funfcn,x,l,u,verbosity,options,f,JAC,YDATA,caller,Jstr,computeLambda,varargin{:});
else
[x,FVAL,JACOB,EXITFLAG,OUTPUT] = ...
nlsq(funfcn,x,verbosity,options,f,JAC,YDATA,caller,varargin{:});
end
Resnorm = FVAL'*FVAL; % assumes FVAL still a vector
if Resnorm > 10*optimget(options,'TolFun',1e-4) & verbosity>0
if verbosity > 0
disp('Optimizer is stuck at a minimum that is not a root')
disp('Try again with a new starting guess')
end
EXITFLAG = -1;
end
% Reset FVAL to shape of the user-function output, fuser
FVAL = reshape(FVAL,size(fuser));
% end FSOLVE 2.*
else % version 1.5 FSOLVE
if length(options)<5;
options(5)=0;
end
% Switch methods making Gauss Newton the default method.
if options(5)==0; options(5)=1; else options(5)=0; end
% Convert to inline function as needed.
if ~isempty(FUN)
[funfcn, msg] = fcnchk(FUN,length(varargin));
if ~isempty(msg)
error(msg);
end
else
error('FUN must be a function name or valid expression.')
end
if ~isempty(GRADFUN)
[gradfcn, msg] = fcnchk(GRADFUN,length(varargin));
if ~isempty(msg)
error(msg);
end
else
gradfcn = [];
end
[x,options] = nlsqold(funfcn,x,options,gradfcn,varargin{:});
if options(8)>10*options(3) & options(1)>0
disp('Optimizer is stuck at a minimum that is not a root')
disp('Try again with a new starting guess')
end
% Set the second output argument FVAL to be options as in old calling syntax
FVAL = options;
% end fsolve version 1.5.*
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [calltype, GRADFUN, otherargs] = parse_call(FUN,options,numargin,numargout,otherargs)
% PARSE_CALL Determine which calling syntax is being used: the FSOLVE prior to 2.0, or
% in version 2.0 or later of the Toolbox.
% old call: [X,OPTIONS]=FSOLVE(FUN,X0,OPTIONS,GRADFUN,varargin)
% new call: [X,FVAL,EXITFLAG,OUTPUT,JACOB]=FSOLVE(FUN,X0,OPTIONS,varargin)
if numargout > 2 % [X,FVAL,EXITFLAG,...]=FSOLVE (...)
calltype = 'new';
GRADFUN = [];
elseif isa(FUN,'cell') % FUN == {...}
calltype = 'new';
GRADFUN = [];
elseif ~isempty(options) & isa(options,'double') % OPTIONS == scalar or and array
calltype = 'old';
if length(otherargs) > 0
GRADFUN = otherargs{1};
otherargs = otherargs(2:end);
else
GRADFUN = [];
end
elseif isa(options,'struct') % OPTIONS has fields
calltype = 'new';
GRADFUN = [];
else % Ambiguous
warnstr = sprintf('%s\n%s\n%s\n',...
'Cannot determine from calling sequence whether to use new (2.0 or later) FSOLVE ', ...
'function or grandfathered FSOLVE function. Assuming new syntax; if call was grandfathered', ...
'FSOLVE syntax, this may give unexpected results.');
warning(warnstr)
calltype = 'new';
GRADFUN = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allfcns,msg] = fprefcnchk(funstr,caller,lenVarIn,gradflag)
%PREFCNCHK Pre- and post-process function expression for FUNCHK.
% [ALLFCNS,MSG] = PREFUNCHK(FUNSTR,CALLER,lenVarIn,GRADFLAG) takes
% the (nonempty) expression FUNSTR from CALLER with LenVarIn extra arguments,
% parses it according to what CALLER is, then returns a string or inline
% object in ALLFCNS. If an error occurs, this message is put in MSG.
%
% ALLFCNS is a cell array:
% ALLFCNS{1} contains a flag
% that says if the objective and gradients are together in one function
% (calltype=='fungrad') or in two functions (calltype='fun_then_grad')
% or there is no gradient (calltype=='fun'), etc.
% ALLFCNS{2} contains the string CALLER.
% ALLFCNS{3} contains the objective function
% ALLFCNS{4} contains the gradient function (transpose of Jacobian).
%
% NOTE: we assume FUNSTR is nonempty.
% Initialize
msg='';
allfcns = {};
funfcn = [];
gradfcn = [];
if gradflag
calltype = 'fungrad';
else
calltype = 'fun';
end
% {fun}
if isa(funstr, 'cell') & length(funstr)==1
% take the cellarray apart: we know it is nonempty
if gradflag
calltype = 'fungrad';0
end
[funfcn, msg] = fcnchk(funstr{1},lenVarIn);
if ~isempty(msg)
error(msg);
end
% {fun,[]}
elseif isa(funstr, 'cell') & length(funstr)==2 & isempty(funstr{2})
if gradflag
calltype = 'fungrad';
end
[funfcn, msg] = fcnchk(funstr{1},lenVarIn);
if ~isempty(msg)
error(msg);
end
% {fun, grad}
elseif isa(funstr, 'cell') & length(funstr)==2 % and ~isempty(funstr{2})
[funfcn, msg] = fcnchk(funstr{1},lenVarIn);
if ~isempty(msg)
error(msg);
end
[gradfcn, msg] = fcnchk(funstr{2},lenVarIn);
if ~isempty(msg)
error(msg);
end
calltype = 'fun_then_grad';
if ~gradflag
warnstr = ...
sprintf('%s\n%s\n%s\n','Jacobian function provided but OPTIONS.Jacobian=''off'';', ...
' ignoring Jacobian function and using finite-differencing.', ...
' Rerun with OPTIONS.Jacobian=''on'' to use Jacobian function.');
warning(warnstr);
calltype = 'fun';
end
elseif ~isa(funstr, 'cell') %Not a cell; is a string expression, function name string or inline object
[funfcn, msg] = fcnchk(funstr,lenVarIn);
if ~isempty(msg)
error(msg);
end
if gradflag % gradient and function in one function/M-file
gradfcn = funfcn; % Do this so graderr will print the correct name
end
else
errmsg = sprintf('%s\n%s', ...
'FUN must be a function name, valid string expression, or inline object;', ...
' or, FUN may be a cell array that contains these type of objects.');
error(errmsg)
end
allfcns{1} = calltype;
allfcns{2} = caller;
allfcns{3} = funfcn;
allfcns{4} = gradfcn;
allfcns{5}=[];
SOLVE.m
function varargout = solve(varargin)
%求各种类型方程组的解析解,
%当找不到解析解时,将自动寻求原点附近的一个近似解, 并可达任意精度.
%用法
% solve('方程1','方程2',...,'方程N','变量1','变量2',...,'变量M')
%SOLVE Symbolic solution of algebraic equations.
% SOLVE('eqn1','eqn2',...,'eqnN')
% SOLVE('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
% SOLVE('eqn1','eqn2',...,'eqnN','var1','var2',...'varN')
%
% The eqns are symbolic expressions or strings specifying equations. The
% vars are symbolic variables or strings specifying the unknown variables.
% SOLVE seeks zeros of the expressions or solutions of the equations.
% If not specified, the unknowns in the system are determined by FINDSYM.
% If no analytical solution is found and the number of equations equals
% the number of dependent variables, a numeric solution is attempted.
%
% Three different types of output are possible. For one equation and one
% output, the resulting solution is returned, with multiple solutions to
% a nonlinear equation in a symbolic vector. For several equations and
% an equal number of outputs, the results are sorted in lexicographic
% order and assigned to the outputs. For several equations and a single
% output, a structure containing the solutions is returned.
%
% Examples:
%
% solve('p*sin(x) = r') chooses 'x' as the unknown and returns
%
% ans =
% asin(r/p)
%
% [x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0') returns
%
% x =
% [ 1]
% [ 3]
%
% y =
% [ 1]
% [ -3/2]
%
% S = solve('x^2*y^2 - 2*x - 1 = 0','x^2 - y^2 - 1 = 0') returns
% the solutions in a structure.
%
% S =
% x: [8x1 sym]
% y: [8x1 sym]
%
% [u,v] = solve('a*u^2 + v^2 = 0','u - v = 1') regards 'a' as a
% parameter and solves the two equations for u and v.
%
% S = solve('a*u^2 + v^2','u - v = 1','a,u') regards 'v' as a
% parameter, solves the two equations, and returns S.a and S.u.
%
% [a,u,v] = solve('a*u^2 + v^2','u - v = 1','a^2 - 5*a + 6') solves
% the three equations for a, u and v.
%
% [x,y] = solve('sin(x+y)-exp(x)*y = 0','x^2-y = 2') cannot find
% an analytic solution, so returns a numeric solution.
%
% See also DSOLVE.
% Copyright (c) 1993-98 by The MathWorks, Inc.
% $Revision: 1.13 $ $Date: 1997/11/29 01:06:35 $
% Set the Explicit solution (for equations of order 4 or less)
% option on in the Maple Workspace.
maple('_EnvExplicit := true;');
% Collect input arguments together in either equation or variable lists.
eqns = [];
vars = [];
for k = 1:nargin
v = varargin{k};
if all(isletter(v) | ('0'<=v & v<='9') | v == '_' | v == ',')
vars = [vars ',' v];
else
[t,stat] = maple(v);
if stat
error(['''' v ''' is not a valid expression or equation.'])
end
eqns = [eqns ',' v];
end
end
if isempty(eqns)
warning('List of equations is empty.')
varargout = cell(1,nargout);
varargout{1} = sym([]);
return
else
eqns(1) = '{'; eqns(end+1) = '}';
end
neqns = sum(commas(eqns)) + 1;
% Determine default variables and sort variable list.
if isempty(vars)
vars = ['[' findsym(sym(eqns),neqns) ']'];
else
vars(1) = '['; vars(end+1) = ']';
end
v = vars;
[vars,stat] = maple('sort',v);
if stat
error(['''' v ''' is not a valid variable list.'])
end
vars(1) = '{'; vars(end) = '}';
nvars = sum(commas(vars)) + 1;
if neqns ~= nvars
warning([int2str(neqns) ' equations in ' int2str(nvars) ' variables.'])
end
% Seek analytic solution.
[R,stat] = maple('solve',eqns,vars);
% If no analytic solution, seek numerical solution.
if (isempty(R) & (nvars == neqns))
[R,stat] = maple('fsolve',eqns,vars);
end
if stat
error(R)
end
% If still no solution, give up.
if isempty(R) | findstr(R,'fsolve')
warning('Explicit solution could not be found.');
varargout = cell(1,nargout);
varargout{1} = sym([]);
return
end
% Expand any RootOf.
while ~isempty(findstr(R,'RootOf'))
k = min(findstr(R,'RootOf'));
p = findstr(R,'{'); p = max(p(p<k));
q = findstr(R,'}'); q = min(q(q>k));
s = R(p:q);
t = s(min(findstr(s,'RootOf'))+6:end);
e = min(find(cumsum((t=='(')-(t==')'))==0));
if isempty(find(commas(t(2:e-1))))
% RootOf with one argument, possibly an imbedded RootOf.
[s,stat] = maple('allvalues',s,'dependent');
if stat
error(s)
end
else
% RootOf with a good estimate as the second argument.
s = maple('evalf',s);
end
R = [R(1:p-1) s R(q+1:end)];
end
% Parse the result.
if nvars == 1 & nargout <= 1
% One variable and at most one output.
% Return a single scalar or vector sym.
S = sym([]);
c = find(commas(R) | R == '}');
for p = find(R == '=')
q = min(c(c>p));
t = trim(R(p+1:q-1)); % The solution (xk)
S = [S; sym(t)];
end
varargout{1} = S;
else
% Several variables.
% Create a skeleton structure.
c = [1 find(commas(vars)) length(vars)];
S = [];
for j = 1:nvars
v = trim(vars(c(j)+1:c(j+1)-1));
S = setfield(S,v,[]);
end
% Complete the structure.
c = [1 find(commas(R) | R == '{' | R == '}') length(R)];
for p = find(R == '=')
q = max(c(c<p));
v = trim(R(q+1:p-1)); % The variable (x)
q = min(c(c>p));
t = trim(R(p+1:q-1)); % The solution (xk)
S = setfield(S,v,[getfield(S,v); sym(t)]);
end
if nargout <= 1
% At most one output, return the structure.
varargout{1} = S;
elseif nargout == nvars
% Same number of outputs as variables.
% Match results in lexicographic order to outputs.
v = fieldnames(S);
for j = 1:nvars
varargout{j} = getfield(S,v{j});
end
else
error([int2str(nvars) ' variables does not match ' ...
int2str(nargout) ' outputs.'])
end
end
%-------------------------
function s = trim(s);
%TRIM TRIM(s) deletes any leading or trailing blanks.
while s(1) == ' ', s(1) = []; end
while s(end) == ' ', s(end) = []; end
%-------------------------
function c = commas(s)
%COMMAS COMMAS(s) is true for commas not inside parentheses.
p = cumsum((s == '(') - (s == ')'));
c = (s == ',') & (p == 0);
NEWTON.m
function x=newton(f_name,x0,tol)
% Purpose: Solve a nonlinear equation by Newton iteration
% Synopsis: x=newton('f_name',x0)
% x=newton('f_name',x0,tol)
%
% x: return a root of f(x)=0 near x0
% f_name:name of the function f(x) that defines the equation
% x0: initial guess
% tol: tolerence(Default:1e-6)
% L.J.Hu 1-8-1998
h=0.0001;M=500;
if nargin<3, tol=1e-6;end
x=x0;xb=x-999;
n=0;
while abs(x-xb)>tol
xb=x;
if n>M break;end
y=feval(f_name, x);
fprintf(' n=%3.0f, x=%12.5e, y=%12.5e, \n',n,x,y)
y_driv=(feval(f_name,x+h)-y)/h;
x=xb-y/y_driv;
n=n+1;
end
fprintf(' n=%3.0f, x=%12.5e, y=%12.5e, ',n,x,y)
if n>500,
fprintf('\n');disp('Warning: iterations exceeds the limitation, probable devergent');
end哈密尔顿回路
TSP模拟退火
accept.m
function y=accept(t,df)
y=(df<0)|(((df/t)<88)&(exp(-df/t)>rand(1,1)));annealing.m
%function R=annealing(N,L,s,t,dt,C,R)
%N为问题规模，即节点个数；L可取较大值，如500、1000；
%s取1、2等；t为初始温度，参考范围为0.5--2；
%dt为衰减因子，一般不小于0.9;
%C为边权矩阵，应是一个强连通图的边权矩阵
%R为初始路径，结果路径也存放在R中
%L、s、t、dt应通过多次试验来确定，以获得优化的结果
%参考《非数值并行算法--模拟退火算法》科学出版社
function R=annealing(N,L,s,t,dt,C,R)
s0=0;
while 1
a=0;
for k=1:L
[r,df]=calculate(R,C,N);
if accept(t,df)
R=r;a=1;
disp(cost_sum(R,C,N));
end
end
t=t*dt
if a==0
s0=s0+1;
else
s0=0;
end
if s0==s
break;
end
end
calculate.m
function [r,df]=calculate(R,C,N)
judge=rand(1,1);
if judge<0.5
r=exchange2(R);
df=cost_sum(r,C,N)-cost_sum(R,C,N);
else
r=exchange3(R);
df=cost_sum(r,C,N)-cost_sum(R,C,N);
end
cost_sum.m
function y=cost_sum(x,C,N)
y=0;
for i=1:(N-1)
y=y+C(x(i),x(i+1));
end
y=y+C(x(N),x(1));
exchange2.m
function r=exchange2(R)
N=length(R);
I=1+fix(unifrnd(0,N));
J=1+fix(unifrnd(0,N-1));
if I==J
J=J+1;
end
if J<I
I=I+J;
J=I-J;
I=I-J;
end
r=R;
if J-I~=1&J-I~=N-1
for p=1:(J-I)
r(I+p)=R(J-p+1);
end
end
if J-I==1
r(I)=R(J);
r(J)=R(I);
end
if J-I==N-1
for p=1:(N-2)
r(p+2)=R(N+1-p)
end
endexchange3.m
function r=exchange3(R)
N=length(R);
K=fix(unifrnd(0,N));
J=fix(unifrnd(0,N-1));
I=fix(unifrnd(0,N-2));
if I==J
J=J+1;
end
if I==K
K=K+2;
end
if J==K
K=K+1;
end
if I>J
I=I+J;
J=I-J;
I=I-J;
end
if I>K
I=I+K;
K=I-K;
I=I-K;
end
if J>K
J=J+K;
K=J-K;
J=J-K;
end
r=R;
if J-I~=1&K-J~=1&K-I~=N-1
for q=1:(J-I)
r(I+q)=R(J+1-q);
end
for q=1:(K-J)
r(J+q)=R(K+1-q);
end
end
if J-I==1&K-J==1
r(K)=R(J);r(J)=R(K);
end
if J-I==1&K-J~=1&K-I~=N-1
for q=1:(K-J)
r(I+q)=R(I+1+q);
end
r(K)=R(J);
end
if K-J==1&J-I~=1&K~=N
for q=1:(J-I)
r(I+1+q)=R(I+q);
end
r(I+1)=R(K);
end
if I==1&J==2&K==N
for q=1:(N-2)
r(1+q)=R(2+q);
end
r(N)=R(2);
end
if I==1&J==(N-1)&K==N
for q=1:(N-2)
r(q)=R(1+q);
end
r(N-1)=R(1);
end
if J-I~=1&K-I==N-1
for q=1:(J-1)
r(q)=R(1+q);
end
r(J)=R(1);
end
if J==(N-1)&K==N&J-I~=1
r(J+1)=R(N);
for q=1:(N-J-1)
r(J+1+q)=R(J+q);
end
end
三边交换简单算法
bianquan.m
N=13;
for i=1:N
for j=1:N
C(i,j)=inf;
end
end
for i=1:N
C(i,i)=0;
end
C(1,2)=6.0;C(1,13)=12.9;
C(2,3)=5.9;C(2,4)=10.3;
C(3,4)=12.2;C(3,5)=17.6;
C(4,13)=8.8;C(4,7)=7.4;
C(4,5)=11.5;
C(5,2)=17.6;C(5,6)=8.2;
C(6,9)=14.9;C(6,7)=20.3;
C(7,9)=19.0;C(7,8)=7.3;
C(8,9)=8.1;C(8,13)=9.2;
C(9,10)=10.3;
C(10,11)=7.7;
C(11,12)=7.2;
C(12,13)=7.9;
for i=1:N
for j=1:N
if C(i,j) < inf
C(j,i)=C(i,j);
end
end
end
for i=1:N
C(i,i)=0;
end
R=[4 7 6 5 3 2 1 13 12 11 10 9 8];cost_sum.m
function y=cost_sum(x,C,N)
y=0;
for i=1:(N-1)
y=y+C(x(i),x(i+1));
end
y=y+C(x(N),x(1));
jiaohuan3.m
R=1:N;
n=0;
for I=1:(N-2)
for J=(I+1):(N-1)
for K=(J+1):N
n=n+1;
Z(n,:)=[I J K];
end
end
end
for m=1:(N*(N-1)*(N-2)/6)
I=Z(m,1);J=Z(m,2);K=Z(m,3);
r=R;
if J-I~=1&K-J~=1&K-I~=N-1
for q=1:(J-I)
r(I+q)=R(J+1-q);
end
for q=1:(K-J)
r(J+q)=R(K+1-q);
end
end
if J-I==1&K-J==1
r(K)=R(J);r(J)=R(K);
end
if J-I==1&K-J~=1&K-I~=N-1
for q=1:(K-J)
r(I+q)=R(I+1+q);
end
r(K)=R(J);
end
if K-J==1&J-I~=1&K~=N
for q=1:(J-I)
r(I+1+q)=R(I+q);
end
r(I+1)=R(K);
end
if I==1&J==2&K==N
for q=1:(N-2)
r(1+q)=R(2+q);
end
r(N)=R(2);
end
if I==1&J==(N-1)&K==N
for q=1:(N-2)
r(q)=R(1+q);
end
r(N-1)=R(1);
end
if J-I~=1&K-I==N-1
for q=1:(J-1)
r(q)=R(1+q);
end
r(J)=R(1);
end
if J==(N-1)&K==N&J-I~=1
r(J+1)=R(N);
for q=1:(N-J-1)
r(J+1+q)=R(J+q);
end
end
if cost_sum(r,C,N)<cost_sum(R,C,N)
R=r
end
end
fprintf('总长为%f\n',cost_sum(R,C,N))提供一种求解最优哈密尔顿的算法---三边交换调整法，要求在运行jiaohuan3（三交换法）之前，给定邻接矩阵C和节点个数N，结果路径存放于R中。
bianquan.m文件给出了一个参数实例，可在命令窗口中输入bianquan,得到邻接矩阵C和节点个数N以及一个任意给出的路径R,，回车后再输入jiaohuan3,得到了最优解。
由于没有经过大量的实验，又是近似算法，对于网络比较复杂的情况，可以尝试多运行几次jiaohuan3,看是否能到进一步的优化结果。
画等温线
dengwen.m
clear;close all;
fphn=fopen('hunan.txt','r');
hnb=fgetl(fphn);
hnmap=fscanf(fphn,'%f %f',[2,59]); % It has 59 rows now.湖南省界经纬度
fclose(fphn);
hnmap=hnmap';
xa=hnmap(:,[1]);
ya=hnmap(:,[2]);
fp=fopen('LATLON57.txt','r');
LL57=fscanf(fp,'%d %f %f',[3,97]); % It has 97 rows now.湖南省97县名称号码，经纬度
fclose(fp);
LL57=LL57';
x=LL57(:,[3])/10;
y=LL57(:,[2])/10;
fpy=fopen('etw00100.txt','r');
ymd57=fscanf(fpy,'%d',[3,1]);%实在不懂这句是什么意思，并且后面也没有用到这个变量
yu97=fscanf(fpy,'%d %f %f',[3,97]); % It has 97 rows now.湖南省97县温度
fclose(fpy);
yu97=yu97';
z=yu97(:,[2]);%湖南省97县温度
hold on;
plot(xa,ya,'.','markersize',5,'color','red');%湖南省界
plot(x,y,'.','markersize',6);%湖南省97县位置
[xi,yi]=meshgrid(linspace(min(x),max(x),25),linspace(min(y),max(y),25));
zi=griddata(x,y,z,xi,yi,'cubic');%湖南省97县温度立方插值
hold on;[c,h]=contour(xi,yi,zi,'b-');%温度等值线，注意是平面的等值线
clabel(c,h);hold off;%三个txt数据文件的意义及其数据格式举例
%97县名称号码、纬度、经度
%LATLON57.txt：
57554 294 1101
57562 296 1114
57574 294 1124
57584 294 1130
57682 287 1136
57662 289 1116
57671 289 1124
57687 282 1128
57649 283 1097
57669 284 1112
57655 258 1104
57745 275 1096
57774 275 1122
57766 273 1114
57776 273 1128
57872 269 1126
57853 267 1106
57845 262 1098
57866 262 1116
57972 258 1130
57965 255 1116
57544 295 1095
57643 290 1098
57642 287 1096
57640 286 1095
57646 286 1100
57657 283 1102
57740 280 1096
57558 291 1105
57564 294 1111
57565 297 1118
57566 295 1117
57577 294 1122
57661 289 1115
57663 289 1120
57674 286 1124
57666 285 1121
57575 295 1126
57680 288 1131
57585 295 1135
57673 287 1129
57679 282 1130
57678 283 1126
57688 282 1126
57773 279 1128
57771 279 1125
57772 278 1125
57780 279 1132
57781 277 1135
57779 270 1134
57882 268 1136
57886 265 1138
57746 276 1100
57658 280 1102
57743 279 1098
57752 279 1106
57744 274 1092
57754 274 1102
57842 269 1097
57841 266 1097
57763 277 1120
57761 278 1113
57760 277 1114
57762 277 1117
57860 270 1113
57768 273 1115
57769 272 1117
57758 271 1116
57767 271 1120
57846 266 1102
57857 264 1103
57851 265 1109
57871 270 1124
57777 272 1129
57778 271 1130
57875 269 1125
57870 268 1121
57874 264 1124
57876 264 1129
57868 266 1119
57867 264 1113
57962 260 1117
57971 259 1122
57966 256 1120
57975 254 1122
57969 253 1113
59063 252 1115
57881 267 1133
57887 261 1131
57889 261 1140
57981 260 1134
57973 258 1127
57974 256 1124
57978 253 1126
57976 254 1129
57985 256 1137
57865 264 1116
%边界经纬度
%hunan.txt：
59
113.5 29.8
113.7 29.5
113.7 29.1
113.8 29.1
114.1 28.8
114.1 28.5
114.3 28.3
113.6 27.5
113.7 27.3
113.9 27.3
113.9 26.7
114.1 26.6
113.9 26.2
114.2 26.2
114.2 26.1
114.0 26.0
114.0 25.6
113.8 25.4
113.3 25.5
112.9 25.4
112.9 25.3
113.0 25.3
113.0 25.0
112.8 25.0
112.7 25.2
112.2 25.3
112.2 24.8
111.5 24.6
111.4 25.2
111.3 25.2
111.2 25.0
111.1 25.0
111.1 25.2
111.4 25.6
111.3 26.3
110.8 26.3
110.4 26.0
110.2 26.1
110.2 26.3
110.0 25.9
109.8 25.9
109.8 26.0
109.4 26.3
109.6 27.1
109.0 27.0
108.9 27.1
109.4 27.6
109.3 27.8
109.2 29.1
109.2 29.5
109.9 29.8
110.3 29.8
110.5 29.7
110.7 29.8
110.6 30.1
111.2 30.0
111.9 29.8
112.3 29.5
112.8 29.8
%97县温度最高值、最低值
%etw00100.txt：
2001 8 1
57554 22.6 12.0
57562 24.3 10.6
57574 24.7 10.0
57584 26.2 19.2
57682 25.2 4.2
57662 24.7 9.8
57671 24.4 0.4
57687 25.0 16.5
57649 22.2 46.2
57669 21.7 2.7
57655 22.7 23.8
57745 23.3 2.8
57774 24.9 0.0
57766 25.5 0
57776 17.8 0.2
57872 25.7 0.1
57853 24.7 1.1
57845 23.5 1.1
57866 25.0 5.1
57972 24.3 37.7
57965 24.6 22.0
57544 22.1 59.5
57643 22.5 43.6
57642 23.6 49.8
57640 23.2 33.4
57646 23.7 27.9
57657 23.1 13.7
57740 22.0 13.5
57558 22.7 20.5
57564 24.0 22.0
57565 24.8 4.2
57566 0 16.1
57577 0 0.9
57661 23.3 6.7
57663 23.9 0.9
57674 25.3 0.0
57666 23.9 0.0
57575 0 8.4
57680 0 0.2
57585 23.6 13.7
57673 25.0 0.3
57679 0 0
57678 24.4 0.1
57688 26.8 0.0
57773 25.2 2.1
57771 24.8 2.9
57772 25.4 1.3
57780 25.1 1.3
57781 25.8 0.0
57779 25.4 0.0
57882 25.8 0.1
57886 23.2 53.6
57746 23.4 6.5
57658 23.1 14.5
57743 0 8.2
57752 22.3 21.6
57744 24.0 0.8
57754 24.0 0.9
57842 25.1 0.2
57841 24.2 2.1
57763 25.2 0.0
57761 0 0.0
57760 23.6 0.0
57762 24.3 0.6
57860 24.5 5.0
57768 0 0.0
57769 24.7 0.0
57758 25.2 0.2
57767 25.4 0.3
57846 24.4 0.0
57857 23.8 0.1
57851 24.6 9.3
57871 25.5 0.9
57777 25.4 0.0
57778 0 0.0
57875 0 0.0
57870 25.0 15.5
57874 25.2 0.4
57876 24.4 74.2
57868 26.5 7.0
57867 25.0 42.4
57962 25.1 56.0
57971 24.3 16.1
57966 24.2 57.6
57975 24.0 9.2
57969 23.4 39.0
59063 23.6 18.2
57881 25.1 9.4
57887 23.9 44.2
57889 22.0 17.4
57981 23.9 0
57973 23.6 84.6
57974 24.2 37.0
57978 23.9 0
57976 23.7 23.2
57985 22.0 21.0
57865 0 57.4
4589
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
