octave-symbolic 1.1.0-3 / usr / share / octave / packages / symbolic-1.1.0 / symfsolve.m

## Copyright (C) 2003 Willem J. Atsma <watsma@users.sf.net>
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{x}, @var{inf}, @var{msg}] =} symfsolve (@dots{})
## Solve a set of symbolic equations using @command{fsolve}. There are a number of
## ways in which this function can be called.
##
## This solves for all free variables, initial values set to 0:
##
## @example
## symbols
## x=sym("x");   y=sym("y");
## f=x^2+3*x-1;  g=x*y-y^2+3;
## a = symfsolve(f,g);
## @end example
##
## This solves for x and y and sets the initial values to 1 and 5 respectively:
##
## @example
## a = symfsolve(f,g,x,1,y,5);
## a = symfsolve(f,g,@{x==1,y==5@});
## a = symfsolve(f,g,[1 5]);
## @end example
##
## In all the previous examples vector a holds the results: x=a(1), y=a(2).
## If initial conditions are specified with variables, the latter determine
## output order:
##
## @example
## a = symfsolve(f,g,@{y==1,x==2@});  # here y=a(1), x=a(2)
## @end example
##
## The system of equations to solve for can be given as separate arguments or
## as a single cell-array:
##
## @example
## a = symfsolve(@{f,g@},@{y==1,x==2@});  # here y=a(1), x=a(2)
## @end example
##
## If the variables are not specified explicitly with the initial conditions,
## they are placed in alphabetic order. The system of equations can be comma-
## separated or given in a cell-array. The return-values are those of
## fsolve; @var{x} holds the found roots.
## @end deftypefn
## @seealso{fsolve}

function [ x,inf,msg ] = symfsolve (varargin)

  ## separate variables and equations
  eqns = cell();
  vars = cell();

  if iscell(varargin{1})
    if !strcmp(typeinfo(varargin{1}{1}),"ex")
      error("First argument must be (a cell-array of) symbolic expressions.")
    endif
    eqns = varargin{1};
    arg_count = 1;
  else
    arg_count = 0;
    for i=1:nargin
      tmp = disp(varargin{i});
      if( iscell(varargin{i}) || ...
          all(isalnum(tmp) || tmp=="_" || tmp==",") || ...
          !strcmp(typeinfo(varargin{i}),"ex") )
        break;
      endif
      eqns{end+1} = varargin{i};
      arg_count = arg_count+1;
    endfor
  endif
  neqns = length(eqns);
  if neqns==0
    error("No equations specified.")
  endif

  ## make a list with all variables from equations
  tmp=eqns{1};
  for i=2:neqns
    tmp = tmp+eqns{i};
  endfor
  evars = findsymbols(tmp);
  nevars=length(evars);

  ## After the equations may follow initial values. The formats are:
  ##   [0 0.3 -3 ...]
  ##   x,0,y,0.3,z,-3,...
  ##   {x==0, y==0.3, z==-3 ...}
  ##   none - default of al zero initial values

  if arg_count==nargin
    vars = evars;
    nvars = nevars;
    X0 = zeros(nvars,1);
  elseif (nargin-arg_count)>1
    if mod(nargin-arg_count,2)
      error("Initial value symbol-value pairs don't match up.")
    endif
    for i=(arg_count+1):2:nargin
      tmp = disp(varargin{i});
      if all(isalnum(tmp) | tmp=="_" | tmp==",")
        vars{end+1} = varargin{i};
        X0((i-arg_count+1)/2)=varargin{i+1};
      else
        error("Error in symbol-value pair arguments.")
      endif
    endfor
    nvars = length(vars);
  else
    nvars = length(varargin{arg_count+1});
    if nvars!=nevars
      error("The number of initial conditions does not match the number of free variables.")
    endif
    if iscell(varargin{arg_count+1})
      ## cell-array of relations - this should work for a list of strings ("x==3") too.
      for i=1:nvars
        tmp = disp(varargin{arg_count+1}{i});
        vars{end+1} = sym (strtok (tmp, "=="));
        X0(i) = str2num(tmp((findstr(tmp,"==")+2):length(tmp)));
      endfor
    else
      ## straight numbers, match up with symbols in alphabetic order
      vars = evars;
      X0 = varargin{arg_count+1};
    endif
  endif

  ## X0 is now a vector, vars a list of variables.
  ## create temporary function:
  symfn = sprintf("function Y=symfn(X) ");
  for i=1:nvars
    symfn = [symfn sprintf("%s=X(%d); ",disp(vars{i}),i)];
  endfor
  for i=1:neqns
    symfn = [symfn sprintf("Y(%d)=%s; ",i,disp(eqns{i}))];
  endfor
  symfn = [symfn sprintf("endfunction")];

  eval(symfn);
  [x,inf,msg] = fsolve("symfn",X0);

endfunction

%!shared x,y,f,g
%! x = sym ("x");
%! y = sym ("y");
%! f = x ^ 2 + 3 * x - 1;
%! g = x * y - y ^ 2 + 3;
%!assert (symfsolve (f, g), [0.30278; -1.58727], 1e-5);
%!assert (symfsolve (f, g, x, 1, y, 5), [0.30278; 1.89004]', 1e-5);
%!assert (symfsolve (f, g, {x==1,y==5}), [0.30278; 1.89004]', 1e-5);
%!assert (symfsolve (f, g, [1 5]), [0.30278; 1.89004]', 1e-5);
%!assert (symfsolve ({f, g}, {y==1,x==2}), [1.89004; 0.30278]', 1e-5);