gusucode.com > symbolic工具箱matlab源码程序 > symbolic/mfun.m
function y = mfun(fun,varargin) %MFUN Numeric evaluation of a special function. % MFUN will be removed in a future release. % Instead, use the appropriate special function (see MFUNLIST). % For example, use bernoulli(n) instead of mfun('bernoulli',n). % % MFUN('function',p1,p2,...,pk) numerically evaluates one of % the symbolic engine's special mathematical functions. The 'function' % input is the name of the function to evaluate and the p's are numeric % inputs to 'function'. The function MFUNLIST displays the allowed % function names and their inputs. % The last parameter specified may be a matrix. % % Example: % x = 0:0.1:5.0; % y = mfun('FresnelC',x) % % See also MFUNLIST, SYMENGINE. % Copyright 1993-2015 The MathWorks, Inc. warning(message('symbolic:mfun:FunctionToBeRemoved')); if isequal(lower(fun),'hypergeom') y = hypergeom(varargin{:}); return end oldDigits = digits(16); cleanupObj = onCleanup(@() digits(oldDigits)); [fun,varargin] = map2mupfun(fun,varargin); a = computeA(varargin{:}); [x,siz,nans] = computeX(varargin{:}); y = formatY(x,oldDigits); [r,st] = evalFun(y,a,fun); if st == 0 y = processResult(r,siz,nans); else y = processSingularity(r,siz,nans,fun,y,a); end function a = computeA(varargin) a = []; for k = 1:length(varargin)-1 t = varargin{k}; if ~isnumeric(t) t = str2num(t); %#ok if isempty(t) error(message('symbolic:mfun:errmsg1')) end end a = [a ',' char(sym(t))]; %#ok<AGROW> end if ~isempty(a), a = [a(2:end) ',']; end function [x,siz,nans] = computeX(varargin) x = varargin{length(varargin)}; if ~isnumeric(x) x = str2num(x); %#ok if isempty(x) error(message('symbolic:mfun:errmsg1')) end end siz = size(x); x = x(:).'; nans = isnan(x); x(nans) = 0; function y = formatY(x,d) % format arguments for integer and real x if all(imag(x) == 0) if all(x == fix(x)) form = ['%' int2str(d) '.0f,']; else form = ['%' int2str(d+6) '.' int2str(d) 'e,']; end y = sprintf(form,x); % format arguments for complex x else form = ['%' int2str(d+6) '.' int2str(d) 'e']; y = sprintf([form '#' form '*1i,'],[real(x); abs(imag(x))]); p = find(y == '#'); % add the correct signs for imaginary parts s = find(imag(x) >= 0); if any(s) y(p(s)) = char('+'*ones(1,length(s))); end s = find(imag(x) < 0); if ~isempty(s) && any(s) y(p(s)) = char('-'*ones(1,length(s))); end end % additional format for the MEX file y(length(y)) = []; y = ['[' y ']']; function [r,st] = evalFun(y,a,fun) [r,st] = evalin(symengine,['float(eval(map(' y ', s->' fun '(' a 's))))']); function y = processResult(r,siz,nans) r = double(r); if isempty(r) error(message('symbolic:mfun:errmsg4')) end r(nans) = NaN; y = reshape(r,siz); function y = processSingularity(r,siz,nans,fun,y,a) % singularities in r if ~isempty(strfind(r,'division by zero')) || ~isempty(strfind(r,'NaN')) || ~isempty(strfind(r,'Singularity')) y = [',' y(2:length(y)-1) ',']; % use commas as delimiters for ALL elements c = find(y == ','); u = NaN*ones(1,prod(siz)); for k = 2:length(c) r = y( c(k-1)+1 : c(k)-1 ); [r,st] = evalin(symengine,['float(' fun '(' a r '))']); if st == 0 u(k-1) = double(r); end end u(nans) = NaN; y = reshape(u,siz); % Overflow elseif ~isempty(strfind(r,'too large')) y = [',' y(2:length(y)-1) ',']; % use commas as delimiters for ALL elements c = find(y == ','); u = Inf*ones(1,prod(siz)); for k = 2:length(c) r = y( c(k-1)+1 : c(k)-1 ); [r,st] = evalin(symengine,['float(' fun '(' a r '))']); if st==0, u(k-1) = double(r); end end u(nans) = NaN; y = reshape(u,siz); else error(message('symbolic:mfun:errmsg2',r)) end % Change function names and argument order or normalization for MuPAD function [fun,args] = map2mupfun(fun,args) fun(1) = lower(fun(1)); switch fun case 'ellipticF' args{1} = asin(args{1}); args{2} = args{2}.^2; case {'ellipticK','ellipticCK','ellipticCE'} args{1} = args{1}.^2; case 'ellipticE' if length(args) == 1 args{1} = args{1}.^2; else args{1} = asin(args{1}); args{2} = args{2}.^2; end case 'ellipticPi' if length(args) == 3 args = {args{2},asin(args{1}),args{3}.^2}; else args{2} = args{2}.^2; end case 'ellipticCPi' args{2} = args{2}.^2; case {'psi','erfc'} if length(args) == 2 args = fliplr(args); end case 'ei' fun = 'Ei'; case 'li' fun = '((z)->Ei(ln(z)))'; case 'ci' fun = 'cosint'; case 'si' fun = 'sinint'; case 'chi' fun = 'Chi'; case 'shi' fun = 'Shi'; case 'w' fun = 'lambertW'; case 'gAMMA' if length(args) == 1 fun = 'gamma'; else fun = 'igamma'; end case 'harmonic' fun = '((z)->psi(z+1)+eulergamma)'; case 'lnGAMMA' fun = '((z)->ln(gamma(z)))'; case 'ssi' fun = '((z)->sinint(z)-PI/2)'; case 'zeta' if length(args) == 3 error(message('symbolic:mfun:ZetavNotSupported')); else args=fliplr(args); end % polynomials case 't' fun = 'orthpoly::chebyshev1'; case 'u' fun = 'orthpoly::chebyshev2'; case 'g' fun = 'orthpoly::gegenbauer'; case 'h' fun = 'orthpoly::hermite'; case 'p' if length(args)==2 fun = 'orthpoly::legendre'; else fun = 'orthpoly::jacobi'; end case 'l' fun = 'laguerreL'; if length(args)==2 args = [args(1) 0 args(2)]; end otherwise % assume matches MuPAD end