[Pkg-octave-commit] [SCM] octave-symbolic branch, master, updated. a718b53403d9f164f8c2a3df521385a01d28a11b
watsma
watsma at 416fae20-06d0-4450-9b69-c6c34d4b5f03
Mon Jan 3 03:40:37 UTC 2011
The following commit has been merged in the master branch:
commit 4c752c479a27209cce99a02c360f1c076e1c2431
Author: watsma <watsma at 416fae20-06d0-4450-9b69-c6c34d4b5f03>
Date: Mon May 5 21:43:57 2003 +0000
main/symbolic/numden.cc
Initial version - Oct file to get numerator and denominator of a
ratio of polynomials.
main/symbolic/findsymbols.cc
Initial version - Oct file to retrieve a sorted list of symbols in
expression.
main/symbolic/findsym.m
Initial version - M file to retrieve list of symbols, m****b-compatible.
main/symbolic/symlsolve.cc
Initial version - Oct file to call GiNaC::lsolve() - solve linear system
of equations.
main/symbolic/symfsolve.m
Initial version - M file generates temporary function file and applies
fsolve() to a (set of) symbolic equations.
main/symbolic/syminfo.cc
Initial version - Oct file to retrieve information about expressions -
calls GiNaC::info().
main/symbolic/Makefile
Added new targets.
git-svn-id: https://octave.svn.sourceforge.net/svnroot/octave/trunk/octave-forge/main/symbolic@923 416fae20-06d0-4450-9b69-c6c34d4b5f03
diff --git a/findsym.m b/findsym.m
new file mode 100644
index 0000000..11e7c4f
--- /dev/null
+++ b/findsym.m
@@ -0,0 +1,87 @@
+## Copyright (C) 2003 Willem J. Atsma
+##
+## This file is part of Octave.
+##
+## Octave 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, or (at your option) any later version.
+##
+## Octave 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 Octave; see the file COPYING. If not,
+## write to the Free Software Foundation, 59 Temple Place -
+## Suite 330, Boston, MA 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} @var{vars} = findsym (@var{f}, @var{n})
+## Find symbols in expression @var{f} and return them comma-separated in
+## string @var{vars}. The symbols are sorted in alphabetic order. If @var{n}
+## is specified, the @var{n} symbols closest to "x" are returned.
+##
+## Example:
+## @example
+## symbols
+## x=sym("x"); y=sym("y"); f=x^2+3*x*y-y^2;
+## vars = findsym (f);
+## vars2 = findsym (f,1);
+## @end example
+##
+## This is intended for m****b compatibility, calls findsymbols().
+## @end deftypefn
+## @seealso{findsymbols}
+
+## Author: Willem J. Atsma <watsma(at)users.sf.net>
+## Created: 18 April 2003
+
+function VARS = findsym(F,Nout)
+
+symlist = findsymbols(F);
+Nlist = length(symlist);
+if Nlist==0
+ warning("No symbols were found.")
+ VARS = "";
+ return
+endif
+
+if exist("Nout")!=1
+ VARS = disp(symlist{1});
+ for i=2:Nlist
+ VARS = [VARS "," disp(symlist{i})];
+ endfor
+ return
+else
+ # If Nout is specified, sort anew from x.
+ symstrings = disp(symlist{1});
+ for i=2:Nlist
+ symstrings = [symstrings ; disp(symlist{i})];
+ endfor
+
+ symasc = toascii(symstrings);
+
+ if Nlist<Nout
+ warning("Asked for %d, variables, only %d found.",Nout,Nlist);
+ Nout=Nlist;
+ endif
+ symasc(:,1) = abs(toascii("x")-symasc(:,1));
+
+ # Sort by creating an equivalent number for each entry
+ Nc = length(symasc(1,:));
+ powbase=zeros(Nc,1); powbase(Nc)=1;
+ for i=(Nc-1):-1:1
+ powbase(i) = powbase(i+1)*128;
+ endfor
+ [xs,I]=sort(symasc*powbase);
+
+ VARS = deblank(symstrings(I(1),:));
+
+ for i=2:Nout
+ VARS = [VARS "," deblank(symstrings(I(i),:))];
+ endfor
+
+endif
diff --git a/findsymbols.cc b/findsymbols.cc
new file mode 100644
index 0000000..9dc4209
--- /dev/null
+++ b/findsymbols.cc
@@ -0,0 +1,97 @@
+/*
+Copyright (C) 2003 Willem J. Atsma
+
+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, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+*/
+#include <octave/oct.h>
+
+#include <ginac/ginac.h>
+#include "ov-vpa.h"
+#include "ov-ex.h"
+#include "symbols.h"
+
+
+/* Travel down the hierarchical expression and insert new symbols in the
+ list so the endresult is sorted. */
+static void append_symbols(octave_value_list& symlist,const GiNaC::ex& expression)
+{
+ int i, j, n = expression.nops();
+
+ for(i=0;i<n;i++) {
+ if(GiNaC::is_exactly_a<GiNaC::symbol>(expression.op(i))) {
+ bool unique = true;
+ int insert_here = symlist.length();
+ GiNaC::ex ex_sym;
+ GiNaC::symbol sym, sym_new = GiNaC::ex_to<GiNaC::symbol>(expression.op(i));
+ std::string sym_name,sym_name_new = sym_new.get_name();
+ for(j=0;j<symlist.length();j++) {
+ /* have to convert back to compare: */
+ get_symbol(symlist(j),ex_sym);
+ sym = GiNaC::ex_to<GiNaC::symbol>(ex_sym);
+ if(sym==sym_new) {
+ unique = false;
+ break;
+ } else {
+ if(sym.get_name()>sym_name_new) {
+ insert_here = j;
+ break;
+ }
+ }
+ }
+ if(unique) {
+ octave_value_list tmp = symlist;
+ symlist.resize(symlist.length()+1);
+ symlist(insert_here) = octave_value(new octave_ex(expression.op(i)));
+ for(j=insert_here;j<tmp.length();j++)
+ symlist(j+1) = tmp(j);
+ }
+ } else append_symbols(symlist,expression.op(i));
+ }
+}
+
+
+DEFUN_DLD(findsymbols,args, ,
+"-*- texinfo -*-\n\
+ at deftypefn {Loadable Function} {@var{syms} =} findsymbols(@var{f})\n\
+\n\
+Returns the symbols in symbolic expression @var{f} in a list.\n\
+The list is sorted in alphabetical order.\n\
+ at end deftypefn")
+{
+ GiNaC::ex expression;
+ octave_value retval;
+ octave_value_list symlist;
+ int nargin = args.length();
+
+ if (nargin != 1) {
+ error("Need one argument.");
+ return retval;
+ }
+ try {
+ if (!get_expression (args(0), expression)) {
+ error("Argument must be a symbolic expression.");
+ return retval;
+ }
+ /* Add 1 to so this works for symbols too. */
+ append_symbols(symlist,expression+1);
+ retval = symlist;
+ } catch (std::exception &e) {
+ error (e.what ());
+ retval = octave_value ();
+ }
+
+ return retval;
+}
diff --git a/numden.cc b/numden.cc
new file mode 100644
index 0000000..f3e50f9
--- /dev/null
+++ b/numden.cc
@@ -0,0 +1,54 @@
+/*
+Copyright (C) 2003 Willem J. Atsma
+
+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, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+*/
+#include <octave/oct.h>
+
+#include <ginac/ginac.h>
+#include "ov-vpa.h"
+#include "ov-ex.h"
+#include "symbols.h"
+
+
+DEFUN_DLD(numden,args,nargout,
+"-*- texinfo -*-\n\
+ at deftypefn Loadable Function {[num,den] =} numden(@var{f})\n\
+\n\
+Return the numerator and denominator of symbolic expression @var{f}.\n\
+ at end deftypefn")
+{
+ GiNaC::ex expression, numden_list;
+ octave_value_list retval;
+ int nargin = args.length();
+ if (nargin != 1) {
+ print_usage ("numden");
+ return retval;
+ }
+ try {
+ if (!get_expression (args(0), expression)) {
+ error("Argument must be a symbolic expression.");
+ return retval;
+ }
+ numden_list = expression.numer_denom();
+ retval.append(new octave_ex(numden_list[0]));
+ retval.append(new octave_ex(numden_list[1]));
+ } catch(std::exception &e) {
+ error (e.what ());
+ retval = octave_value ();
+ }
+ return retval;
+}
diff --git a/symfsolve.m b/symfsolve.m
new file mode 100644
index 0000000..8d8a5b1
--- /dev/null
+++ b/symfsolve.m
@@ -0,0 +1,161 @@
+## Copyright (C) 2003 Willem J. Atsma
+##
+## 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, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[ @var{x}, at var{inf}, at var{msg} ] =} symfsolve (...)
+## Solve a set of symbolic equations using 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 list/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 list or cell-array. The return-values are those of
+## fsolve; @var{x} holds the found roots.
+## @end deftypefn
+## @seealso{fsolve}
+
+## Author: Willem J. Atsma <watsma(at)users.sf.net>
+##
+## 2003-04-22 Willem J. Atsma <watsma(at)users.sf.net>
+## * Initial revision
+
+function [ x,inf,msg ] = symfsolve (varargin)
+
+ #separate variables and equations
+ eqns = list;
+ vars = list;
+
+ if iscell(varargin{1}) | islist(varargin{1})
+ if !strcmp(typeinfo(varargin{1}{1}),"ex")
+ error("First argument must be (a cell-array/list 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}) | islist(varargin{i}) | ...
+ all(isalnum(tmp) | tmp=="_" | tmp==",") | ...
+ !strcmp(typeinfo(varargin{i}),"ex") )
+ break;
+ endif
+ eqns=append(eqns,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=append(vars,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}) | islist(varargin{arg_count+1})
+ # List/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 = append(vars,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
diff --git a/syminfo.cc b/syminfo.cc
new file mode 100644
index 0000000..ce90d9c
--- /dev/null
+++ b/syminfo.cc
@@ -0,0 +1,103 @@
+/*
+Copyright (C) 2003 Willem J. Atsma
+
+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, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+*/
+
+#include <octave/oct.h>
+#include <ginac/ginac.h>
+#include "ov-ex.h"
+#include "ov-relational.h"
+#include "ov-struct.h"
+#include "ov-bool.h"
+#include "symbols.h"
+
+
+inline static void
+add_info_field(Octave_map& info,const GiNaC::ex& exp,std::string field,int info_flag)
+{
+ bool flag = exp.info(info_flag);
+ info.assign(field,octave_value(flag));
+}
+
+static void get_info(Octave_map& oct_info, const GiNaC::ex& gobj)
+{
+ add_info_field(oct_info,gobj,"numeric",GiNaC::info_flags::numeric);
+ add_info_field(oct_info,gobj,"real",GiNaC::info_flags::real);
+ add_info_field(oct_info,gobj,"rational",GiNaC::info_flags::rational);
+ add_info_field(oct_info,gobj,"integer",GiNaC::info_flags::integer);
+ add_info_field(oct_info,gobj,"crational",GiNaC::info_flags::crational);
+ add_info_field(oct_info,gobj,"cinteger",GiNaC::info_flags::cinteger);
+ add_info_field(oct_info,gobj,"relation",GiNaC::info_flags::relation);
+ add_info_field(oct_info,gobj,"symbol",GiNaC::info_flags::symbol);
+ add_info_field(oct_info,gobj,"list",GiNaC::info_flags::list);
+ add_info_field(oct_info,gobj,"exprseq",GiNaC::info_flags::exprseq);
+ add_info_field(oct_info,gobj,"polynomial",GiNaC::info_flags::polynomial);
+ add_info_field(oct_info,gobj,"integer_polynomial",GiNaC::info_flags::integer_polynomial);
+ add_info_field(oct_info,gobj,"cinteger_polynomial",GiNaC::info_flags::cinteger_polynomial);
+ add_info_field(oct_info,gobj,"rational_polynomial",GiNaC::info_flags::rational_polynomial);
+ add_info_field(oct_info,gobj,"crational_polynomial",GiNaC::info_flags::crational_polynomial);
+ add_info_field(oct_info,gobj,"rational_function",GiNaC::info_flags::rational_function);
+ add_info_field(oct_info,gobj,"algebraic",GiNaC::info_flags::algebraic);
+ add_info_field(oct_info,gobj,"indexed",GiNaC::info_flags::indexed);
+ add_info_field(oct_info,gobj,"has_indices",GiNaC::info_flags::has_indices);
+ add_info_field(oct_info,gobj,"idx",GiNaC::info_flags::idx);
+}
+
+DEFUN_DLD(syminfo,args,nargout,
+"-*- texinfo -*-\n\
+ at deftypefn Loadable Function {@var{info} =} syminfo(@var{eqn})\n\
+\n\
+Returns information regarding the nature of symbolic expression @var{eqn} in\n\
+a structure. Uses the GiNaC::info() function.\n\
+ at end deftypefn")
+{
+ Octave_map oct_info;
+ octave_value retval;
+ int nargin = args.length();
+
+ if (nargin != 1) {
+ error("Wrong number of arguments.");
+ return retval = oct_info;
+ }
+
+ try {
+ GiNaC::ex expression;
+ if(!get_expression(args(0),expression)) {
+ GiNaC::numeric num;
+ if(get_numeric(args(0),num))
+ expression = GiNaC::ex(num);
+ else {
+ GiNaC::relational relation;
+ if(get_relation(args(0),relation)) {
+ expression = GiNaC::ex(relation);
+ std::cout << expression << std::endl;
+ } else {
+ if(!get_symbol(args(0),expression)) {
+ error("Expected a symbolic object as parameter.");
+ return retval = oct_info;
+ }
+ }
+ }
+ }
+ get_info(oct_info,expression);
+ retval = oct_info;
+ } catch (std::exception &e) {
+ error (e.what ());
+ retval = octave_value ();
+ }
+ return retval;
+}
diff --git a/symlsolve.cc b/symlsolve.cc
new file mode 100644
index 0000000..8280615
--- /dev/null
+++ b/symlsolve.cc
@@ -0,0 +1,94 @@
+/*
+Copyright (C) 2003 Willem J. Atsma
+
+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, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Based on differentiate.cc by Ben Sapp (2002).
+*/
+#include <octave/oct.h>
+
+#include <ginac/ginac.h>
+#include "ov-ex.h"
+#include "ov-relational.h"
+#include "symbols.h"
+
+DEFUN_DLD(symlsolve,args,nargout,
+"-*- texinfo -*-\n\
+ at deftypefn Loadable Function {@var{sols} =} symlsolve(@var{eqns}, at var{vars})\n\
+\n\
+Apply the GiNaC lsolve() method to the given linear system of equations and\n\
+variables. @var{sols}, @var{eqns} and @var{vars} must be symbolic expressions\n\
+or lists of these.\n\
+ at end deftypefn")
+{
+ GiNaC::lst eqns, vars;
+ GiNaC::relational relation;
+ GiNaC::ex expression, sols;
+ octave_value_list oct_sols;
+ octave_value retval;
+ int i, nargin = args.length();
+
+ if (nargin != 2) {
+ error("Need 2 arguments.");
+ return retval;
+ }
+
+ try {
+ if(args(0).is_list()) {
+ octave_value_list oct_eqn_list(args(0).list_value());
+ for(i=0;i<oct_eqn_list.length();i++) {
+ if(!get_relation(oct_eqn_list(i),relation)) {
+ if(!get_expression(oct_eqn_list(i),expression)) {
+ error("Equation %d is not a valid expression.",i+1);
+ return retval;
+ } else relation = (expression==0);
+ }
+ eqns.append(relation);
+ }
+ } else {
+ if(!get_relation(args(0),relation)) {
+ if(!get_expression(args(0),expression)) {
+ error("Equation is not a valid expression.");
+ return retval;
+ } else relation = (expression==0);
+ }
+ eqns.append(relation);
+ }
+
+ if(args(1).is_list()) {
+ octave_value_list oct_vars(args(1).list_value());
+ for(i=0;i<oct_vars.length();i++) {
+ if(!get_symbol(oct_vars(i),expression)) {
+ error("Variable %d is not a valid symbol.",i+1);
+ return retval;
+ }
+ vars.append(expression);
+ }
+ } else {
+ if(!get_symbol(args(1),expression)) {
+ error("Variable is not a valid symbol.");
+ return retval;
+ }
+ vars.append(expression);
+ }
+ sols = lsolve(eqns,vars);
+ for(i=0;i<(int)sols.nops();i++) oct_sols.append(new octave_ex(sols[i]));
+ retval = oct_sols;
+ } catch (std::exception &e) {
+ error (e.what ());
+ retval = octave_value ();
+ }
+ return retval;
+}
--
octave-symbolic
More information about the Pkg-octave-commit
mailing list