From octave-sources-request at bevo dot che dot wisc dot edu  Wed Jul 26 04:44:40 2000
Subject: invroots : inversion of polynomial 'roots'-function (incl demo)
From: Rolf Fabian <fabian at tu-cottbus dot de>
To: "'octave-sources UWISC'" <octave-sources at bevo dot che dot wisc dot edu>
Date: Tue, 25 Jul 2000 11:43:57 +0200

Hi,

the attached ,invroots' function  is something missing in octave's polynomial toolbox !

	P = invroots ( x0 )

it inverts the available ,roots' function in the following sense, that it
returns the normalized polynomial P of degree N
( i.e. highest exponents coefficient  set to 1), 
which	PASSES though ALL ROOTS
given by the N elements of input vector  x0


#------- snip -------- snip ------- snip --------#
#USAGE	P = invroots( x0 )
#
#	inversion of polynomial 'roots' function
#
#	returns coefficients of polynomial P ( at order N )
#		equivalent to the product
#	P(x) = ( x-x0(1) )*( x-x0(2) )* ... *( x-x0(N) )
#	which passes through all N roots given by vector x0
#NOTES	* roots(invroots(x0))	|- (permuted) x0
#	* invroots(roots(C*P))	|-	P
#	  where  C denotes an arbitrary scale factor
#	* always P(1)=1, P(2)=-sum(x0), P(N+1)=prod(-x0)
#ASSOC	roots
#EXA	roots(invroots(x0=(-5:5)+i))  |-  permuted x0
#AUTHOR	(C) Rolf Fabian 000725 <fabian at tu-cottbus dot de>

function P=invroots(x0)

	if	nargin!=1, error("argument number"); end

	#requested by 'roots.m', which defines roots(0)=[]
	if	isempty(x0), P=0; return; end

	[R,C]=size(x0);
	N=max([R,C]);
	if	all([R,C]!=1), error("input x0 must be a vector of roots"); end

	P=1;
	for	k=1:N
		P=sum([P,0;0,-x0(k)*P]);
	end

	if	R!=1
		P=reshape(P,R+1,C);
	end
#------- snip -------- snip ------- snip --------#

DEMO:

#------- snip -------- snip ------- snip --------#
#USAGE	invrootsX()
#	graphical 'invroots'-demo
#AUTHOR	(C) Rolf Fabian 000725 <fabian at tu-cottbus dot de>

P=invroots(x0=-2:0.1:2);

SAMPLES=5001;			#many samples necessary for best graph.
				#display of x-axis croosings at roots x0

t=linspace(-2.2,2.2,SAMPLES)';

disp('P=invroots(x0=-2:0.1:2)');
disp('result polynomial passing through all 41 equidistant roots x0:');
disp('P=');
disp(P.');

gset yrange [3e-6:1e10];
xlabel('x'); ylabel('y');
title('EXA  polynomial P(x)=invroots(x0=-2:0.1:2) evaluation (degree 41)');

semilogy(t, Z=polyval(P,t),';pos. branch;',t,-Z,';neg. branch;')

#------- snip -------- snip ------- snip --------#




-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------