octave_packages/symbolic-1.1.0/symfsolve.m

   1 ## Copyright (C) 2003 Willem J. Atsma <watsma@users.sf.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{x}, @var{inf}, @var{msg}] =} symfsolve (@dots{})
  18 ## Solve a set of symbolic equations using @command{fsolve}. There are a number of
  19 ## ways in which this function can be called.
  20 ##
  21 ## This solves for all free variables, initial values set to 0:
  22 ##
  23 ## @example
  24 ## symbols
  25 ## x=sym("x");   y=sym("y");
  26 ## f=x^2+3*x-1;  g=x*y-y^2+3;
  27 ## a = symfsolve(f,g);
  28 ## @end example
  29 ##
  30 ## This solves for x and y and sets the initial values to 1 and 5 respectively:
  31 ##
  32 ## @example
  33 ## a = symfsolve(f,g,x,1,y,5);
  34 ## a = symfsolve(f,g,@{x==1,y==5@});
  35 ## a = symfsolve(f,g,[1 5]);
  36 ## @end example
  37 ##
  38 ## In all the previous examples vector a holds the results: x=a(1), y=a(2).
  39 ## If initial conditions are specified with variables, the latter determine
  40 ## output order:
  41 ##
  42 ## @example
  43 ## a = symfsolve(f,g,@{y==1,x==2@});  # here y=a(1), x=a(2)
  44 ## @end example
  45 ##
  46 ## The system of equations to solve for can be given as separate arguments or
  47 ## as a single cell-array:
  48 ##
  49 ## @example
  50 ## a = symfsolve(@{f,g@},@{y==1,x==2@});  # here y=a(1), x=a(2)
  51 ## @end example
  52 ##
  53 ## If the variables are not specified explicitly with the initial conditions,
  54 ## they are placed in alphabetic order. The system of equations can be comma-
  55 ## separated or given in a cell-array. The return-values are those of
  56 ## fsolve; @var{x} holds the found roots.
  57 ## @end deftypefn
  58 ## @seealso{fsolve}
  59
  60 function [ x,inf,msg ] = symfsolve (varargin)
  61
  62   ## separate variables and equations
  63   eqns = cell();
  64   vars = cell();
  65
  66   if iscell(varargin{1})
  67     if !strcmp(typeinfo(varargin{1}{1}),"ex")
  68       error("First argument must be (a cell-array of) symbolic expressions.")
  69     endif
  70     eqns = varargin{1};
  71     arg_count = 1;
  72   else
  73     arg_count = 0;
  74     for i=1:nargin
  75       tmp = disp(varargin{i});
  76       if( iscell(varargin{i}) || ...
  77           all(isalnum(tmp) || tmp=="_" || tmp==",") || ...
  78           !strcmp(typeinfo(varargin{i}),"ex") )
  79         break;
  80       endif
  81       eqns{end+1} = varargin{i};
  82       arg_count = arg_count+1;
  83     endfor
  84   endif
  85   neqns = length(eqns);
  86   if neqns==0
  87     error("No equations specified.")
  88   endif
  89
  90   ## make a list with all variables from equations
  91   tmp=eqns{1};
  92   for i=2:neqns
  93     tmp = tmp+eqns{i};
  94   endfor
  95   evars = findsymbols(tmp);
  96   nevars=length(evars);
  97
  98   ## After the equations may follow initial values. The formats are:
  99   ##   [0 0.3 -3 ...]
 100   ##   x,0,y,0.3,z,-3,...
 101   ##   {x==0, y==0.3, z==-3 ...}
 102   ##   none - default of al zero initial values
 103
 104   if arg_count==nargin
 105     vars = evars;
 106     nvars = nevars;
 107     X0 = zeros(nvars,1);
 108   elseif (nargin-arg_count)>1
 109     if mod(nargin-arg_count,2)
 110       error("Initial value symbol-value pairs don't match up.")
 111     endif
 112     for i=(arg_count+1):2:nargin
 113       tmp = disp(varargin{i});
 114       if all(isalnum(tmp) | tmp=="_" | tmp==",")
 115         vars{end+1} = varargin{i};
 116         X0((i-arg_count+1)/2)=varargin{i+1};
 117       else
 118         error("Error in symbol-value pair arguments.")
 119       endif
 120     endfor
 121     nvars = length(vars);
 122   else
 123     nvars = length(varargin{arg_count+1});
 124     if nvars!=nevars
 125       error("The number of initial conditions does not match the number of free variables.")
 126     endif
 127     if iscell(varargin{arg_count+1})
 128       ## cell-array of relations - this should work for a list of strings ("x==3") too.
 129       for i=1:nvars
 130         tmp = disp(varargin{arg_count+1}{i});
 131         vars{end+1} = {sym(strtok(tmp,"=="))};
 132         X0(i) = str2num(tmp((findstr(tmp,"==")+2):length(tmp)));
 133       endfor
 134     else
 135       ## straight numbers, match up with symbols in alphabetic order
 136       vars = evars;
 137       X0 = varargin{arg_count+1};
 138     endif
 139   endif
 140
 141   ## X0 is now a vector, vars a list of variables.
 142   ## create temporary function:
 143   symfn = sprintf("function Y=symfn(X) ");
 144   for i=1:nvars
 145     symfn = [symfn sprintf("%s=X(%d); ",disp(vars{i}),i)];
 146   endfor
 147   for i=1:neqns
 148     symfn = [symfn sprintf("Y(%d)=%s; ",i,disp(eqns{i}))];
 149   endfor
 150   symfn = [symfn sprintf("endfunction")];
 151
 152   eval(symfn);
 153   [x,inf,msg] = fsolve("symfn",X0);
 154
 155 endfunction
 156
 157 %!shared
 158 % x = sym ("x");
 159 % y = sym ("y");
 160 % f = x ^ 2 + 3 * x - 1;
 161 % g = x * y - y ^ 2 + 3;
 162 %!test
 163 % assert (symfsolve (f, g), [0.30278; -1.58727]', 1e-5);
 164 %!test
 165 % assert (symfsolve (f, g, x, 1, y, 5), [0.30278; 1.89004]', 1e-5);
 166 %!test
 167 % assert (symfsolve (f, g, {x==1,y==5}), [0.30278; 1.89004]', 1e-5);
 168 %!test
 169 % assert (symfsolve (f, g, [1 5]), [0.30278; 1.89004]', 1e-5);
 170 %!test
 171 % assert (symfsolve ({f, g}, {y==1,x==2}), [1.89004; 0.30278]', 1e-5);