--- /dev/null
+## 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 = sym ("x");
+% y = sym ("y");
+% f = x ^ 2 + 3 * x - 1;
+% g = x * y - y ^ 2 + 3;
+%!test
+% assert (symfsolve (f, g), [0.30278; -1.58727]', 1e-5);
+%!test
+% assert (symfsolve (f, g, x, 1, y, 5), [0.30278; 1.89004]', 1e-5);
+%!test
+% assert (symfsolve (f, g, {x==1,y==5}), [0.30278; 1.89004]', 1e-5);
+%!test
+% assert (symfsolve (f, g, [1 5]), [0.30278; 1.89004]', 1e-5);
+%!test
+% assert (symfsolve ({f, g}, {y==1,x==2}), [1.89004; 0.30278]', 1e-5);
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff