Mercurial > hg > octave-max
changeset 1025:f558749713f1
[project @ 1995-01-11 20:52:10 by jwe]
| author | jwe |
|---|---|
| date | Wed, 11 Jan 1995 20:52:10 +0000 |
| parents | 56520a75b5b3 |
| children | 9fc405c8c06c |
| files | scripts/polynomial/compan.m scripts/polynomial/conv.m scripts/polynomial/deconv.m scripts/polynomial/poly.m scripts/polynomial/polyderiv.m scripts/polynomial/polyinteg.m scripts/polynomial/polyreduce.m scripts/polynomial/polyval.m scripts/polynomial/polyvalm.m scripts/polynomial/residue.m scripts/polynomial/roots.m |
| diffstat | 11 files changed, 416 insertions(+), 202 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/polynomial/compan.m +++ b/scripts/polynomial/compan.m @@ -1,6 +1,25 @@ -function A = compan(c) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. -# compan (c) +function A = compan (c) + +# usage: compan (c) +# # Compute the companion matrix corresponding to polynomial vector c. # # In octave a polynomial is represented by it's coefficients (arranged @@ -24,26 +43,24 @@ # # SEE ALSO: poly, roots, residue, conv, deconv, polyval, polyderiv, polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. - if(nargin != 1) - usage ("compan(vector)"); + if (nargin != 1) + usage ("compan (vector)"); endif - if(!is_vector(c)) + if(is_matrix (c)) error("compan: expecting a vector argument."); endif - # Ensure that c is a row vector. +# Ensure that c is a row vector. + if(rows(c) > 1) c = c.'; endif - n = length(c); - A = diag(ones(n-2,1),-1); - A(1,:) = -c(2:n)/c(1); + n = length (c); + A = diag (ones (n-2, 1), -1); + A (1, :) = -c (2:n) /c (1); endfunction
--- a/scripts/polynomial/conv.m +++ b/scripts/polynomial/conv.m @@ -1,23 +1,42 @@ +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function y = conv (a, b) +# usage: conv (a, b) +# # Convolve two vectors. +# # y = conv (a, b) returns a vector of length equal to length (a) + # length (b) -1. +# # If a and b are polynomial coefficient vectors, conv returns the # coefficients of the product polynomial. # # SEE ALSO: deconv, poly, roots, residue, polyval, polyderiv, polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. if (nargin != 2) - usage (" conv(a,b)"); + usage ("conv(a, b)"); endif - if (is_matrix(a) || is_matrix(b)) + if (is_matrix (a) || is_matrix (b)) error("conv: both arguments must be vectors"); endif @@ -37,14 +56,14 @@ # Use the shortest vector as the coefficent vector to filter. if (la < lb) if (ly > lb) - x = [b zeros (1, ly - lb)]; + x = [b, zeros (1, ly - lb)]; else x = b; endif y = filter (a, 1, x); else if(ly > la) - x = [a zeros (1, ly - la)]; + x = [a, zeros (1, ly - la)]; else x = a; endif
--- a/scripts/polynomial/deconv.m +++ b/scripts/polynomial/deconv.m @@ -1,5 +1,25 @@ +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function [b, r] = deconv (y, a) +# usage: deconv (y, a) +# # Deconvolve two vectors. # # [b, r] = deconv (y, a) solves for b and r such that @@ -12,13 +32,10 @@ # SEE ALSO: conv, poly, roots, residue, polyval, polyderiv, # polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. if (nargin != 2) - usage (" deconv (y,a)"); + usage ("deconv (y, a)"); endif if (is_matrix (y) || is_matrix (a))
--- a/scripts/polynomial/poly.m +++ b/scripts/polynomial/poly.m @@ -1,13 +1,37 @@ +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function y = poly (x) +# usage: poly (x) +# # If A is a square n-by-n matrix, poly (A) is the row vector of # the coefficients of det (z * eye(n) - A), the characteristic # polynomial of A. +# # If x is a vector, poly (x) is a vector of coefficients of the # polynomial whose roots are the elements of x. -# Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993 -# Copyright Dept of Probability Theory and Statistics TU Wien +# Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) Dec 24, 1993. + + if (nargin != 1) + usage ("poly (x)"); + endif m = min (size (x)); n = max (size (x)); @@ -18,7 +42,7 @@ elseif (m == n) v = eig (x); else - usage ("poly(x), where x is a vector or a square matrix"); + usage ("poly (x), where x is a vector or a square matrix"); endif y = [1, zeros (1, n)];
--- a/scripts/polynomial/polyderiv.m +++ b/scripts/polynomial/polyderiv.m @@ -1,27 +1,43 @@ -function p = polyderiv(p) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. -# polyderiv(c) +function p = polyderiv (p) + +# usage: polyderiv (p) +# # Returns the coefficients of the derivative of the polynomial whose -# coefficients are given by vector c. +# coefficients are given by vector p. # # SEE ALSO: poly, polyinteg, polyreduce, roots, conv, deconv, residue, # filter, polyval, polyvalm -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. - if(nargin != 1) - usage ("polyderiv(vector)"); + if (nargin != 1) + usage ("polyderiv (vector)"); endif - if(is_matrix(p)) - error("argument must be a vector"); + if (is_matrix (p)) + error ("argument must be a vector"); endif - lp = length(p); - if(lp == 1) + lp = length (p); + if (lp == 1) p = 0; return; elseif (lp == 0) @@ -29,6 +45,6 @@ return; end - p = p(1:(lp-1)) .* [(lp-1):-1:1]; + p = p (1:(lp-1)) .* [(lp-1):-1:1]; endfunction
--- a/scripts/polynomial/polyinteg.m +++ b/scripts/polynomial/polyinteg.m @@ -1,39 +1,55 @@ -function p = polyinteg(p) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. -# polyinteg(c) +function p = polyinteg (p) + +# usage: polyinteg (p) +# # Returns the coefficients of the integral the polynomial whose coefficients -# are represented by the vector c. +# are represented by the vector p. # # The constant of integration is zero. # # SEE ALSO: poly, polyderiv, polyreduce, roots, conv, deconv, residue, # filter, polyval, polyvalm -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. if(nargin != 1) - usage ("polyinteg(vector)"); + usage ("polyinteg (vector)"); endif - if(is_matrix(p)) - error("argument must be a vector"); + if (is_matrix (p)) + error ("argument must be a vector"); endif - lp = length(p); + lp = length (p); - if(lp == 0) + if (lp == 0) p = []; return; end - if(rows(p) > 1) - # Convert to column vector + if (rows (p) > 1) +# Convert to column vector p = p.'; endif - p = [ p 0 ] ./ [lp:-1:1 1]; + p = [ p, 0 ] ./ [ lp:-1:1, 1 ]; endfunction
--- a/scripts/polynomial/polyreduce.m +++ b/scripts/polynomial/polyreduce.m @@ -1,30 +1,47 @@ -function p = polyreduce(p) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. -# polyreduce(c) +function p = polyreduce (p) + +# usage: polyreduce(c) +# # Reduces a polynomial coefficient vector to a minimum number of terms, # i.e. it strips off any leading zeros. # # SEE ALSO: poly, roots, conv, deconv, residue, filter, polyval, polyvalm, # polyderiv, polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. + + index = find (p == 0); - index = find(p==0); + index = find (index == 1:length(index)); - index = find(index == 1:length(index)); - - if (length(index) == 0) + if (length (index) == 0) return; endif if(length(p)>1) - p = p(index(length(index))+1:length(p)); + p = p (index (length (index))+1:length(p)); endif - if(length(p)==0) + if (length (p) == 0) p = 0; endif + endfunction
--- a/scripts/polynomial/polyval.m +++ b/scripts/polynomial/polyval.m @@ -1,5 +1,25 @@ -function y = polyval(c,x) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. +function y = polyval (c, x) + +# usage: polyval (c, x) +# # Evaluate a polynomial. # # In octave, a polynomial is represented by it's coefficients (arranged @@ -16,27 +36,25 @@ # SEE ALSO: polyvalm, poly, roots, conv, deconv, residue, filter, # polyderiv, polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. - if(nargin != 2) - usage ("polyval(c,x)"); + if (nargin != 2) + usage ("polyval (c, x)"); endif - if(is_matrix(c)) - error("poly: first argument must be a vector."); + if(is_matrix (c)) + error ("poly: first argument must be a vector."); endif - if(length(c) == 0) + if (length (c) == 0) y = c; return; endif - n = length(c); - y = c(1)*ones(rows(x),columns(x)); + n = length (c); + y = c (1) * ones (rows (x), columns (x)); for index = 2:n - y = c(index) + x .* y; + y = c (index) + x .* y; endfor + endfunction
--- a/scripts/polynomial/polyvalm.m +++ b/scripts/polynomial/polyvalm.m @@ -1,5 +1,25 @@ -function y = polyvalm(c,x) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. +function y = polyvalm (c, x) + +# usage: polyvalm (c, x) +# # Evaluate a polynomial in the matrix sense. # # In octave, a polynomial is represented by it's coefficients (arranged @@ -17,25 +37,22 @@ # SEE ALSO: polyval, poly, roots, conv, deconv, residue, filter, # polyderiv, polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. if(nargin != 2) - usage ("polyvalm(c,x)"); + usage ("polyvalm (c, x)"); endif - if(is_matrix(c)) + if (is_matrix (c)) error("poly: first argument must be a vector."); endif - if(!is_square(x)) + if(! is_square (x)) error("poly: second argument must be a square matrix."); endif [v, d] = eig(x); - y = v * diag(polyval(c,diag(d))) * v'; + y = v * diag (polyval (c, diag (d))) * v'; endfunction
--- a/scripts/polynomial/residue.m +++ b/scripts/polynomial/residue.m @@ -1,6 +1,25 @@ -function [r, p, k, e] = residue(b,a,toler) +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. -# [r p k e] = residue(b,a) +function [r, p, k, e] = residue (b, a, toler) + +# usage: [r, p, k, e] = residue (b, a) +# # If b and a are vectors of polynomial coefficients, then residue # calculates the partial fraction expansion corresponding to the # ratio of the two polynomials. The vector r contains the residue @@ -17,23 +36,24 @@ # the r, p, and e vectors) and N is the length of the k vector. # # [r p k e] = residue(b,a,tol) +# # This form of the function call may be used to set a tolerance value. # The default value is 0.001. The tolerance value is used to determine # whether poles with small imaginary components are declared real. It is # also used to determine if two poles are distinct. If the ratio of the -# imaginary part of a pole to the real part is less than tol, the imaginary -# part is discarded. If two poles are farther apart than tol they are -# distinct. +# imaginary part of a pole to the real part is less than tol, the +# imaginary part is discarded. If two poles are farther apart than tol +# they are distinct. # # Example: -# b = [1 1 1]; -# a = [1 -5 8 -4]; +# b = [1, 1, 1]; +# a = [1, -5, 8, -4]; # -# [r p k e] = residue(b,a) +# [r, p, k, e] = residue (b, a) # # returns # -# r = [-2 7 3]; p = [2 2 1]; k = []; e = [1 2 1]; +# r = [-2, 7, 3]; p = [2, 2, 1]; k = []; e = [1, 2, 1]; # # which implies the following partial fraction expansion # @@ -43,10 +63,7 @@ # # SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg -# Author: -# Tony Richardson -# amr@mpl.ucsd.edu -# June 1994 +# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994. # Here's the method used to find the residues. # The partial fraction expansion can be written as: @@ -78,13 +95,13 @@ # # s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2) # -# The coefficients of the polynomials on the right are stored -# in a row vector called rhs, while the coefficients of the -# polynomial on the left is stored in a row vector called lhs. -# If the multiplicity of any root is greater than one we'll -# also need derivatives of this equation of order up to the -# maximum multiplicity minus one. The derivative coefficients -# are stored in successive rows of lhs and rhs. +# The coefficients of the polynomials on the right are stored in a row +# vector called rhs, while the coefficients of the polynomial on the +# left is stored in a row vector called lhs. If the multiplicity of +# any root is greater than one we'll also need derivatives of this +# equation of order up to the maximum multiplicity minus one. The +# derivative coefficients are stored in successive rows of lhs and +# rhs. # # For our example lhs and rhs would be: # @@ -109,153 +126,169 @@ # # We then solve for the residues using matrix division. - if(nargin < 2 || nargin > 3) - usage ("residue(b,a[,toler])"); + if (nargin < 2 || nargin > 3) + usage ("residue (b, a [, toler])"); endif if (nargin == 2) - # Set the default tolerance level toler = .001; endif - # Make sure both polynomials are in reduced form. - a = polyreduce(a); - b = polyreduce(b); +# Make sure both polynomials are in reduced form. + + a = polyreduce (a); + b = polyreduce (b); - b = b/a(1); - a = a/a(1); + b = b / a(1); + a = a / a(1); - la = length(a); - lb = length(b); + la = length (a); + lb = length (b); - # Handle special cases here. - if(la == 0 || lb == 0) +# Handle special cases here. + + if (la == 0 || lb == 0) k = r = p = e = []; return; elseif (la == 1) - k = b/a; + k = b / a; r = p = e = []; return; endif - # Find the poles. - p = roots(a); - lp = length(p); +# Find the poles. + + p = roots (a); + lp = length (p); - # Determine if the poles are (effectively) real. - index = find(abs(imag(p) ./ real(p)) < toler); - if (length(index) != 0) - p(index) = real(p(index)); +# Determine if the poles are (effectively) real. + + index = find (abs (imag (p) ./ real (p)) < toler); + if (length (index) != 0) + p (index) = real (p (index)); endif - # Find the direct term if there is one. - if(lb>=la) - # Also returns the reduced numerator. - [k, b] = deconv(b,a); - lb = length(b); +# Find the direct term if there is one. + + if (lb >= la) +# Also returns the reduced numerator. + [k, b] = deconv (b, a); + lb = length (b); else k = []; endif - if(lp == 1) - r = polyval(b,p); + if (lp == 1) + r = polyval (b, p); e = 1; return; endif - # We need to determine the number and multiplicity of the roots. - # D is the number of distinct roots. - # M is a vector of length D containing the multiplicity of each root. - # pr is a vector of length D containing only the distinct roots. - # e is a vector of length lp which indicates the power in the partial - # fraction expansion of each term in p. +# We need to determine the number and multiplicity of the roots. +# +# D is the number of distinct roots. +# M is a vector of length D containing the multiplicity of each root. +# pr is a vector of length D containing only the distinct roots. +# e is a vector of length lp which indicates the power in the partial +# fraction expansion of each term in p. - # Set initial values. We'll remove elements from pr as we find - # multiplicities. We'll shorten M afterwards. - e = ones(lp,1); - M = zeros(lp,1); +# Set initial values. We'll remove elements from pr as we find +# multiplicities. We'll shorten M afterwards. + + e = ones (lp, 1); + M = zeros (lp, 1); pr = p; - D = 1; M(1) = 1; + D = 1; + M(1) = 1; - old_p_index = 1; new_p_index = 2; - M_index = 1; pr_index = 2; - while(new_p_index<=lp) - if(abs(p(new_p_index) - p(old_p_index)) < toler) - # We've found a multiple pole. - M(M_index) = M(M_index) + 1; - e(new_p_index) = e(new_p_index-1) + 1; - # Remove the pole from pr. - pr(pr_index) = []; + old_p_index = 1; + new_p_index = 2; + M_index = 1; + pr_index = 2; + + while (new_p_index <= lp) + if (abs (p (new_p_index) - p (old_p_index)) < toler) +# We've found a multiple pole. + M (M_index) = M (M_index) + 1; + e (new_p_index) = e (new_p_index-1) + 1; +# Remove the pole from pr. + pr (pr_index) = []; else - # It's a different pole. - D++; M_index++; M(M_index) = 1; - old_p_index = new_p_index; pr_index++; +# It's a different pole. + D++; + M_index++; + M (M_index) = 1; + old_p_index = new_p_index; + pr_index++; endif new_p_index++; endwhile - # Shorten M to it's proper length - M = M(1:D); +# Shorten M to it's proper length - # Now set up the polynomial matrices. + M = M (1:D); + +# Now set up the polynomial matrices. MM = max(M); - # Left hand side polynomial - lhs = zeros(MM,lb); - rhs = zeros(MM,lp*lp); - lhs(1,:) = b; + +# Left hand side polynomial + + lhs = zeros (MM, lb); + rhs = zeros (MM, lp*lp); + lhs (1, :) = b; rhi = 1; dpi = 1; mpi = 1; - while(dpi<=D) + while (dpi <= D) for ind = 1:M(dpi) - if(mpi>1 && (mpi+ind)<=lp) + if (mpi > 1 && (mpi+ind) <= lp) cp = [p(1:mpi-1); p(mpi+ind:lp)]; - elseif (mpi==1) - cp = p(mpi+ind:lp); + elseif (mpi == 1) + cp = p (mpi+ind:lp); else - cp = p(1:mpi-1); + cp = p (1:mpi-1); endif - rhs(1,rhi:rhi+lp-1) = prepad(poly(cp),lp); + rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp); rhi = rhi + lp; endfor - mpi = mpi + M(dpi); + mpi = mpi + M (dpi); dpi++; endwhile - if(MM > 1) + if (MM > 1) for index = 2:MM - lhs(index,:) = prepad(polyderiv(lhs(index-1,:)),lb); + lhs (index, :) = prepad (polyderiv (lhs (index-1, :)), lb); ind = 1; for rhi = 1:lp - cp = rhs(index-1,ind:ind+lp-1); - rhs(index,ind:ind+lp-1) = prepad(polyderiv(cp),lp); + cp = rhs (index-1, ind:ind+lp-1); + rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp); ind = ind + lp; endfor endfor endif - # Now lhs contains the numerator polynomial and as many derivatives as are - # required. rhs is a matrix of polynomials, the first row contains the - # corresponding polynomial for each residue and successive rows are - # derivatives. +# Now lhs contains the numerator polynomial and as many derivatives as +# are required. rhs is a matrix of polynomials, the first row +# contains the corresponding polynomial for each residue and +# successive rows are derivatives. - # Now we need to evaluate the first row of lhs and rhs at each distinct - # pole value. If there are multiple poles we will also need to evaluate - # the derivatives at the pole value also. +# Now we need to evaluate the first row of lhs and rhs at each +# distinct pole value. If there are multiple poles we will also need +# to evaluate the derivatives at the pole value also. - B = zeros(lp,1); - A = zeros(lp,lp); + B = zeros (lp, 1); + A = zeros (lp, lp); dpi = 1; row = 1; - while(dpi<=D) + while (dpi <= D) for mi = 1:M(dpi) - B(row) = polyval(lhs(mi,:),pr(dpi)); + B (row) = polyval (lhs (mi, :), pr (dpi)); ci = 1; for col = 1:lp - cp = rhs(mi,ci:ci+lp-1); - A(row,col) = polyval(cp,pr(dpi)); + cp = rhs (mi, ci:ci+lp-1); + A (row, col) = polyval (cp, pr(dpi)); ci = ci + lp; endfor row++; @@ -263,7 +296,8 @@ dpi++; endwhile - # Solve for the residues. - r = A\B; +# Solve for the residues. + + r = A \ B; endfunction
--- a/scripts/polynomial/roots.m +++ b/scripts/polynomial/roots.m @@ -1,12 +1,31 @@ +# Copyright (C) 1995 John W. Eaton +# +# 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function r = roots (v) +# usage: roots (v) +# # For a vector v with n components, return the roots of the # polynomial v(1) * z^(n-1) + ... + v(n-1) * z + v(n). -# Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993 -# Copyright Dept of Probability Theory and Statistics TU Wien - - [nr, nc] = size(v); +# Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993. + + [nr, nc] = size (v); if (nr <= 1 && nc <= 1) r = []; return; @@ -23,14 +42,14 @@ f = find (v); m = max (size (f)); if (m > 0) - v = v(f(1):f(m)); + v = v (f(1):f(m)); l = max (size (v)); if (l > 1) A = diag (ones (1, l-2), -1); - A(1,:) = -v(2:l) ./ v(1); + A (1, :) = -v (2:l) ./ v (1); r = eig (A); if (f(m) < n) - r = [r; zeros (n - f(m), 1)]; + r = [r; zeros (n - f (m), 1)]; endif else r = zeros (n - f(m), 1);