*****************************************************************************
Geordie La Forge @ http://MeAmI.org
*
* Zero Nine October Two-Thousand and Nine
* M. Michael Musatov
* Algorithm and Natural Pattern Engineering Section
* Symmetry Engineering Department
* College of Earth
* The State University
* Los Angeles, California
*
* M.S. thesis
* Advisor. Dr. Ivan Musatov
*
* Copywrite=(C)==©2009 M.Michael Musatov
* [http://meami.org]
* All Rights Reserved in Perpetuity
* These statement apply to All Intellectual Property=(IP)
*
* 'SUBROUTINE BUDAN' (TM)
* Subroutine uses the Fourier-Budan Theorem to determine
* the number of roots that the alpha polynomial has on the
* interval [u,v].
*
* PARAMETERS: iu = lower bound of alpha interval
* iv = uppper bound of alpha interval
* VARIABLES: coefficient = coefficient of alpha polynomial
* dcoeff = coefficient of polynomial
derivatives
* deriv = derivatives of alpha polynomial
* fvapor = the alpha polynomial
* ia,ib = # of sign changes for derivative
series
* ivapor = alpha = vapor fraction
* jsign,ksign = flags for derivative sign
change
* numroot = number of zeros on the
interval
*****************************************************************************
SUBROUTINE BUDAN(J,Ncomp,coefficient,nroot)
IMPLICIT REAL*8(a--ho--z)
INTEGER p
DIMENSION dcoeff(0:100,0:100), coefficient(0:100), deriv
(0:100)
PARAMETER(iu = 0, iv = 1)
C <i>DATA (coefficient(l), l = 0,Ncomp-1) /-1.,1.,-2.,3.,-4.,5./</
i>
OPEN(unit=2,file='test',status= 'unknown')
REWIND(Unit=2)
ia = 0
ib = 0
do 0500 ivapor = iu, iv, 1
*
* Evaluate the polynomial function at the endpoints iu and iv
*
fvapor = 0d0
do 0600 p = 1, Ncomp
fvapor = fvapor + coefficient(Ncomp-p)*ivapor**(Ncomp-p)
0600 continue
write(2,*) 'fvapor = ',fvapor
write(2,*) ' '
*
* Calculate coefficients of first derivative
*
do 1000 n = Ncomp-1, 0, -1
dcoeff(0,n) = coefficient(n)
write (2,*) 'dcoeff(0,',n,') = ',dcoeff(0,n)
1000 continue
write(2,*) ' '
* Calculate coefficients of 2nd- and higher-order derivatives
* as multiples of those of the first derivative
do 1500 m = 1, Ncomp-1
do 2000 n = Ncomp-m, 1, -1
dcoeff(mn-1) = n*dcoeff(m-l,n)
write (2,*) 'dcoeff(',m,',',n-1,') = ',
@dcoeff(m,n-1)
2000 continue
write(2,*) ' '
1500 continue
*
* Evaluate the derivative series at the endpoints iu and iv
*
do 3000 n = 1, Ncomp-1
deriv(m) = 0.d0
do 4000 n = Ncomp-m, 1, -1
term = dcoeff(m,n-1)*ivapor**(n-1)
if( (n-I) .EQ. 0 ) term = dcoeff(m,n-1)
deriv(m) = deriv(m) + term
write(2,*) 'inter deriv(',m,') = ',deriv(m)
4000 continue
write(2,*) 'total deriv(',m,') = ',deriv(m)
write(2,*) ' '
3000 continue
*
* Count the sign changes between the terms of the series
*
if(fvapor. LT. 0.) then
ksign = 0
else
ksign = 1
end if
write(2,*) 'ksign = ',ksign,' for fvapor'
do 5000 i = 1, Ncomp-1
if(deriv(i) .LT. 0.) then
jsign = 0
else
jsign = 1
end if
write(2,*) 'jsign = ',jsign,' for deriv(',i,')'
*
* Increment A or B, depending upon the endpoint under
evaluation
*
if(ivapor .EQ. iu) then
if(ksign .NE. jsign) then
ia = ia + I
write(2,*) 'ia = ',ia,' for deriv(',i,')'
end if
end if
if(ivapor .EQ. iv) then
if(ksign .NE. jsign) then
ib = ib + 1
write(2,*) 'ib = ',ib,' for deriv(',i,')'
end if
end if
ksign = jsign
write(2,*) 'ksign = ',ksign,' after deriv(',i,')'
write(2,*) ' '
5000 continue
0500 continue
*
* Pass a flag to calling program to indicate root conditions
*
write(2,*) 'ia =',ia,' and ib = ',ib
numroot = ia -ib
write(2,*) 'numroot = ',numroot
write(2,6000) Ncomp-1, numroot, iu, iv, J
6000 format("This polynomial of order ',i3,' has ',i3,' zeros on the
in
@terval [',i2,',',i2,'] for J = ',i3)
CLOSE(unit=2)
return
end
*****************************************************************************
Geordie La Forge @ http://MeAmI.org
SUBROUTINE BUDAN(J,Ncomp,coefficient,nroot)
IMPLICIT REAL*8(a--ho--z)
INTEGER p
fvapor = 0d0
1500 continue
4000 continue
3000 continue
if(fvapor. LT. 0.) then
ia = ia + 1