Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

Diff of /trunk/octave-forge/main/symbolic/inst/sym2poly.m [r8777] .. [r8778] Maximize Restore

  Switch to unified view

a/trunk/octave-forge/main/symbolic/inst/sym2poly.m b/trunk/octave-forge/main/symbolic/inst/sym2poly.m
1
## Copyright (C) 2003 Willem J. Atsma
1
## Copyright (C) 2003 Willem J. Atsma <watsma@users.sf.net>
2
##
2
##
3
## This program is free software; you can redistribute it and/or
3
## This program is free software; you can redistribute it and/or
4
## modify it under the terms of the GNU General Public
4
## modify it under the terms of the GNU General Public
5
## License as published by the Free Software Foundation;
5
## License as published by the Free Software Foundation;
6
## either version 2, or (at your option) any later version.
6
## either version 2, or (at your option) any later version.
...
...
14
## You should have received a copy of the GNU General Public
14
## You should have received a copy of the GNU General Public
15
## License along with this software; see the file COPYING.  If not,
15
## License along with this software; see the file COPYING.  If not,
16
## see <http://www.gnu.org/licenses/>.
16
## see <http://www.gnu.org/licenses/>.
17
17
18
## -*- texinfo -*-
18
## -*- texinfo -*-
19
## @deftypefn {Function File} @var{c} = sym2poly (@var{p}, @var{x})
19
## @deftypefn {Function File} {@var{c} =} sym2poly (@var{p}, @var{x})
20
## Returns the coefficients of the symbolic polynomial expression @var{p}
20
## Returns the coefficients of the symbolic polynomial expression @var{p}
21
## as a vector. If there is only one free variable in @var{p} the
21
## as a vector. If there is only one free variable in @var{p} the
22
## coefficient vector @var{c} is a plain numeric vector. If there is more
22
## coefficient vector @var{c} is a plain numeric vector. If there is more
23
## than one free variable in @var{p}, a second argument @var{x} specifies the
23
## than one free variable in @var{p}, a second argument @var{x} specifies the
24
## free variable and the function returns a cell vector of symbolic expressions.
24
## free variable and the function returns a cell vector of symbolic expressions.
25
## The coefficients correspond to decreasing exponent of the free variable.
25
## The coefficients correspond to decreasing exponent of the free variable.
26
##
26
##
27
## Example:
27
## Example:
28
## @example
28
## @example
29
## symbols
29
## symbols
30
## x = sym("x"); y=sym("y");
30
## x = sym("x");
31
## y = sym("y");
31
## c = sym2poly (x^2+3*x-4);    # c = [1,3,-4]
32
## c = sym2poly (x^2+3*x-4);    # c = [1,3,-4]
32
## c = sym2poly (x^2+y*x,x);    # c = @{2,y,0@}
33
## c = sym2poly (x^2+y*x,x);    # c = @{2,y,0@}
33
## @end example
34
## @end example
34
##
35
##
35
## If @var{p} is not a polynomial the result has no warranty.
36
## If @var{p} is not a polynomial the result has no warranty.
36
##
37
##
37
## @seealso{poly2sym,polyval,roots}
38
## @seealso{poly2sym,polyval,roots}
38
## @end deftypefn
39
## @end deftypefn
39
40
40
## Author: Willem J. Atsma <watsma(at)users.sf.net>
41
## Created: 18 April 2003
41
## Created: 18 April 2003
42
## Changed: 25 April 2003
42
## Changed: 25 April 2003
43
##    Removed the use of differentiate to get to coefficients - round-off
43
##    Removed the use of differentiate to get to coefficients - round-off
44
##     errors cause problems. Now using newly created sumterms().
44
##     errors cause problems. Now using newly created sumterms().
45
## Changed: 6 May 2003
45
## Changed: 6 May 2003
46
##    Removed the attempt to use ldegree(), degree() and coeff() - results
46
##    Removed the attempt to use ldegree(), degree() and coeff() - results
47
##     with these are inconsistent.
47
##     with these are inconsistent.
48
48
49
function c = sym2poly(p,x)
49
function c = sym2poly(p,x)
50
50
51
BADPOLY_COEFF_LIMIT = 500;
51
  BADPOLY_COEFF_LIMIT = 500;
52
52
53
if is_vpa(p)
53
  if is_vpa(p)
54
  # polynomial is one vpa number
54
    ## polynomial is one vpa number
55
  c = to_double(p);
55
    c = to_double(p);
56
  if length(c)!=1
56
    if length(c)!=1
57
      error("Argument is not a polynomial.");
57
      error("Argument is not a polynomial.");
58
  endif
58
    endif
59
  return
59
    return
60
endif
60
  endif
61
61
62
if !is_ex(p)
62
  if !is_ex(p)
63
  error("Argument has to be a symbolic expression.")
63
    error("Argument has to be a symbolic expression.")
64
endif
64
  endif
65
65
66
pvars = findsymbols(p);
66
  pvars = findsymbols(p);
67
if isempty(pvars)
67
  if isempty(pvars)
68
  # It is possible that we get an expression without any symbols.
68
    ## It is possible that we get an expression without any symbols.
69
  c = to_double(p);
69
    c = to_double(p);
70
  return;
70
    return;
71
endif
71
  endif
72
nvars = length(pvars); 
72
  nvars = length(pvars); 
73
73
74
if nvars>1 && exist("x")!=1
74
  if nvars>1 && exist("x")!=1
75
  error("Symbolic expression has more than 1 free variable; no variable specified.")
75
    error("Symbolic expression has more than 1 free variable; no variable specified.")
76
elseif exist("x")!=1
76
  elseif exist("x")!=1
77
  x = pvars{1};
77
    x = pvars{1};
78
endif
78
  endif
79
79
80
p = expand(p);
80
  p = expand(p);
81
81
82
## GiNaC has commands to access coefficients directly, but in octave this often
82
  ## GiNaC has commands to access coefficients directly, but in octave this often
83
## does not work, because for example x^2 typed in octave results in a 
83
  ## does not work, because for example x^2 typed in octave results in a 
84
## non-integer power in GiNaC: x^2.0 .
84
  ## non-integer power in GiNaC: x^2.0 .
85
85
86
[num,den] = numden(p);
86
  [num,den] = numden(p);
87
tmp = findsymbols(den);
87
  tmp = findsymbols(den);
88
for i=1:length(tmp)
88
  for i=1:length(tmp)
89
  if tmp{i}==x
89
    if tmp{i}==x
90
      error("Symbolic expression is a ratio of polynomials.")
90
      error("Symbolic expression is a ratio of polynomials.")
91
  endif
91
    endif
92
endfor
92
  endfor
93
93
94
p = expand(p);
94
  p = expand(p);
95
p_terms = sumterms(p);
95
  p_terms = sumterms(p);
96
# if this is well behaved, I can find the coefficients by dividing with x
96
  ## if this is well behaved, I can find the coefficients by dividing with x
97
c_ex = cell;
97
  c_ex = cell;
98
for i=1:length(p_terms)
98
  for i=1:length(p_terms)
99
  tmp = p_terms{i};
99
    tmp = p_terms{i};
100
  for j=1:BADPOLY_COEFF_LIMIT
100
    for j=1:BADPOLY_COEFF_LIMIT
101
      if disp(differentiate(tmp,x))=="0"
101
      if disp(differentiate(tmp,x))=="0"
102
          break;
102
        break;
103
      endif
103
      endif
104
      tmp = tmp/x;
104
      tmp = tmp/x;
105
  endfor
105
    endfor
106
  if j==BADPOLY_COEFF_LIMIT
106
    if j==BADPOLY_COEFF_LIMIT
107
      printf("Please examine your code or adjust this function.\n");
107
      printf("Please examine your code or adjust this function.\n");
108
      printf("This error may occur because the passed expression is not a polynomial.\n");
108
      printf("This error may occur because the passed expression is not a polynomial.\n");
109
      error("Reached the set limit (%d) for the number of coefficients.",BADPOLY_COEFF_LIMIT)
109
      error("Reached the set limit (%d) for the number of coefficients.",BADPOLY_COEFF_LIMIT)
110
  endif
110
    endif
111
  if (length(c_ex)<j) || isempty(c_ex{j})
111
    if (length(c_ex)<j) || isempty(c_ex{j})
112
      c_ex(j)=tmp;
112
      c_ex(j)=tmp;
113
  else
113
    else
114
      c_ex(j) = c_ex{j}+tmp;
114
      c_ex(j) = c_ex{j}+tmp;
115
  endif       
115
    endif
116
endfor
116
  endfor
117
order = length(c_ex)-1;
117
  order = length(c_ex)-1;
118
118
119
all_numeric = true;
119
  all_numeric = true;
120
for i=1:(order+1)
120
  for i=1:(order+1)
121
  if isempty(c_ex{i})
121
    if isempty(c_ex{i})
122
      c_ex(i) = vpa(0);
122
      c_ex(i) = vpa(0);
123
  endif
123
    endif
124
  cvar=findsymbols(c_ex{i});
124
    cvar=findsymbols(c_ex{i});
125
  ncvar = length(cvar);
125
    ncvar = length(cvar);
126
  if ncvar
126
    if ncvar
127
      all_numeric=false;
127
      all_numeric=false;
128
      for j=1:ncvar
128
      for j=1:ncvar
129
          if disp(cvar{j})==disp(x)
129
        if disp(cvar{j})==disp(x)
130
              printf("Possibly this error occurs because two symbols with the same name\n");
130
          printf("Possibly this error occurs because two symbols with the same name\n");
131
              printf("are different to GiNaC. Make sure the free variable exists as a\n");
131
          printf("are different to GiNaC. Make sure the free variable exists as a\n");
132
              printf("symbol in your workspace.\n");
132
          printf("symbol in your workspace.\n");
133
              error("The symbolic expression is not a polynomial.")
133
          error("The symbolic expression is not a polynomial.")
134
          endif
134
        endif
135
      endfor
135
      endfor
136
  endif
136
    endif
137
endfor
137
  endfor
138
138
139
c_ex = reverse(c_ex);
139
  c_ex = reverse(c_ex);
140
140
141
if all_numeric
141
  if all_numeric
142
  for i=1:(order+1)
142
    for i=1:(order+1)
143
      c(1,i)=to_double(c_ex{i});
143
      c(1,i)=to_double(c_ex{i});
144
  endfor
144
    endfor
145
else
145
  else
146
  c = c_ex;
146
    c = c_ex;
147
endif
147
  endif
148
149
endfunction
148
150
149
%!shared
151
%!shared
150
% symbols
152
% symbols
151
% x=sym("x"); y=sym("y");
153
% x=sym("x"); y=sym("y");
152
%!test
154
%!test