Il 17/12/2012 18:30, Soviet_Mario ha scritto:
> Tutto sommato vorrei solo uno strumento di risoluzione
> simultanea, per valutare rapidamente in vari contesti una miriade di
> formule approssimate che usiamo per trattare gli equilibri.
Dunque, in calce incollo la routine (fortran 95, testata con gcc f95 su
Linux), completa di routine di servizio e programma di test, di cui ti
ho parlato nel post precedente (non so quale dei due arrivera' prima :-)
Dato che spesso copiaincollare codice risulta in programmi
incompilabili, ho messo il file anche qui:
http://trusso.freeshell.org/SovietMario/DoubleZeriRealiPolinomiReali.f
Il modo migliore di capire l'algoritmo (che avevo cercato di descriverti
senza troppa chiarezza) e' di leggerla.
Come dicevo, la routine funziona egregiamente per trovare zeri singoli e
distinti di polinomi "onesti" :-), anche tirando un po' il loro range di
variazione. Per esempio, calcolando il polinomio
(x-.01)*(x-1)*(x-7.3)*(x-18)*(x-100)*(x-10000)*(x-568777)*(x^4+5)
che ha ovviamente 7 radici reali e 4 complesse, e si espande come (non
lo faccio a mano, uso maXima)
1000*x^11-578903310*x^10+5760878110833*x^9-720035853390518*x^8+15866453941479473*x^7-
90040929443886828*x^6+75664925530631965*x^5-4347552244952590*x^4+79332269707372365*x^3-
450204632746851390*x^2+378180605700389000*x-3736864890000000
dando in pasto al programma il grado, 11, e poi i coefficienti (in
ordine canonico):
1000 -578903310 5760878110833 -720035853390518 15866453941479473
-90040929443886828 +75664925530631965 -4347552244952590
79332269707372365 -450204632746851390 378180605700389000 -3736864890000000
vengono trovate correttamente le 7 radici reali.
Nel caso di zeri multipli, anche questo algoritmo (come tutti) da' dei
problemi: la ragione sta principalmente nel fatto che gli zeri multipli
cadono ad un estremo dell'intervallo, che e' stato determinato nel passo
precedente con precisione macchina (16 cifre decimali, nel nostro caso),
dove la derivata era molto prossima allo zero ma non necessariamente
proprio zero. Non sapendo a priori cosa, *nello specifico problema*,
possa essere considerato "sufficientemente vicino" allo zero, la routine
di bisezione agli estremi testa se il polinomio e' esattamente eguale a
zero; se non lo e', non lo considera tale. E, a seconda del segno che il
polinomio assume in quel punto e all'altro estremo dell'intervallo, puo'
trovare per bisezione un altro zero molto vicino, o non trovarne alcuno.
Un altro motivo che puo' far mancare uno zero di tangenza e'
l'arrotondamento effettuato nel calcolo del polinomio: a seconda dei
casi, uno zero multiplo puo' non essere trovato per niente, o essere
determinato come 2 o piu' zeri distinti, anche se molto vicini. Inoltre,
nella costruzione stessa del polinomio a partire dall'espressione
(x-x1)^n1 * (x-x2)^n2... ...(x-xm)^nm
si verificano arrotondamenti, per cui i coefficienti del polinomio che
poi si da' in pasto alla routine non sono quelli voluti, ma coefficienti
di *un altro* polinomio, che magari gli zeri cercati proprio non li ha.
Nei casi piu' semplici, in cui il calcolo del polinomio riesce a dare un
risultato esatto, gli zeri vengono trovati, e con la corretta
molteplicita': per esempio, partendo da
(x-1.5)^10
cioe'
1024*x^10-15360*x^9+103680*x^8-414720*x^7+1088640*x^6-1959552*x^5+2449440*x^4-2099520*x^3+
1180980*x^2-393660*x+59049
dando in input i coefficienti
1024 -15360 +103680 -414720 +1088640 -1959552 +2449440 -2099520 +1180980
-393660 +59049
si trovano proprio 15 zeri pari a 1.5000000;
ma in altri casi, dove gli arrotondamenti acquistano importanza, gli
zeri vengono trovati con molteplicita' difettosa o non trovati del
tutto. Es:
(x-1)^4*(x-0.07)^4*(x^2+3)
cioe'
100000000*x^10-428000000*x^9+1014940000*x^8-1863897200*x^7+2375011201*x^6-1780284404*x^5
+694076809*x^4-121925216*x^3+10512019*x^2-440412*x+7203
grado 10 e coefficienti
100000000 -428000000 +1014940000 -1863897200 +2375011201 -1780284404
+694076809 -121925216 +10512019 -440412 +7203
gli zeri 1 e 0.07 vengono determinati, ma con molteplicita' 2 anziche' 4;
sostituendo 0.07 con 0.007,
(x-1)^4*(x-0.007)^4*(x^2+3)
ciooe'
1000000000000*x^10-4028000000000*x^9+9112294000000*x^8-16253177372000*x^7+
19450651490401*x^6-12536716357604*x^5+3341607973609*x^4-87554106416*x^3+898509619*x^2-4144812*x
+7203
grado 10 e coefficienti
1000000000000 -4028000000000 9112294000000 -16253177372000
19450651490401 -12536716357604 3341607973609 -87554106416 898509619
-4144812 7203
gli zeri a 1 non vengono piu' trovati.
Come dicevo nell'altro post, non credo che i casi che tratti tu siano
cosi' patologici, per cui questa routine dovrebbe esserti utile e
sufficiente: in ogni caso, per vedere cosa dovresti aspettarti, ti
conviene tirar fuori per prima cosa il polinomio e plottarlo.
Se invece effettivamente lo sono, dovresti fornirmi qualche caso reale
su cui ragionare.
Ciao e buon anno
! (fixed width characters are better)
!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
RECURSIVE Subroutine findrealzeroes(N,coeff, Nzeroes,zeroes)
C note: a_0 = coeff(0) is known term,
C a_i = coeff(i) is coefficient of x^i
implicit double precision (a-h,o-z)
implicit integer (i-n)
C ( yes, I'm a VERY old Fortran programmer :-)
dimension coeff(0:N), zeroes(1:N)
double precision,dimension(:),allocatable :: endpoints, dcoeff
logical zerofound
if (coeff(N).eq.0.) STOP "*** Polynomial is NOT of degree N ***"
print *, "**** degree: ",N," ****"
Cauchy = CauchyRadius(N,coeff)
print*, "Cauchy bounds: ",-Cauchy,Cauchy
if (N.eq.1) then
C straight line with non zero a_1
Nzeroes = 1
zeroes(1) = - coeff(0)/coeff(1)
else
C allocate space for derivate polynomial and intervals
allocate ( endpoints(0:N), dcoeff(0:N-1) )
C derive polynomial and normalize derivate polynomial
C so that max degree coeffincent remains the same
do i=0,N-1
dcoeff(i) = ((i+1) * coeff(i+1)) / N
enddo
print *, "Normalized derivate:"
print *, (dcoeff(i), i=N-1,0,-1)
C initialize interval endpoints to -Cauchy,Cauchy,...,Cauchy.
endpoints(0) = -Cauchy
do i=1,N
endpoints(i)=Cauchy
end do
C RECURSIVELY find zeroes of derivate polynomial
call findrealzeroes(N-1,dcoeff, NDrealzeroes,endpoints(1) )
Nintervals = NDrealzeroes + 1
C here endpoints (0:NDrealzeroes+1) contains
C -Cauchy, 1st zero, 2nd zero ... last zero, Cauchy. (ordered)
print *, Nintervals," intervals endpoints, and values there:"
print *, (endpoints(i), i=0,Nintervals)
print *, (poly(N,coeff,endpoints(i)), i=0,Nintervals)
C find *at most one* zero for each interval:
Nzeroes = 0
do i = 1, Nintervals
call bisect(N,coeff,endpoints(i-1),endpoints(i),
$ zerofound,foundzero)
if (zerofound) then
Nzeroes = Nzeroes + 1
zeroes(Nzeroes) = foundzero
endif
end do
C cleanup:
deallocate ( endpoints, dcoeff )
endif
print *, "Degree ",N," real zeroes found, and relative values:"
print *, (zeroes(i), i=1,Nzeroes)
print *, (poly(N,coeff,zeroes(i)), i=1,Nzeroes)
return
end
double precision function CauchyRadius(N,a)
C note: a_0 = a(0) is known term,
C a_i = a(i) is coefficient of x^i
implicit double precision (a-h,o-z)
implicit integer (i-n)
dimension a(0:N)
CauchyRadius = 1.
do i=0,N-1
CauchyRadius = max (CauchyRadius, abs(a(i)/a(N))+1.)
end do
return
end
subroutine bisect(N,coeff,xinit,xend, zerofound,foundzero)
implicit double precision (a-h,o-z)
implicit integer (i-n)
dimension coeff(0:N)
logical zerofound, nearlysame
C simple bisection method: other methods
C would be faster, but in this case who cares...
zerofound = .false.
x1 = xinit
x2 = xend
poly1 = poly(N,coeff,x1)
poly2 = poly(N,coeff,x2)
C print *, "entered bisect: N = ",N,
C $ " x1 = ",x1,poly1," x2 = ",x2,poly2
if (poly1*poly2.gt.0.) return
if(poly1.eq.0. .and. poly2.eq.0.)then
zerofound = .true.
foundzero = x1
return
endif
if(poly1.eq.0.)then
zerofound = .true.
foundzero = x1
return
endif
if(poly2.eq.0.)then
zerofound = .true.
foundzero = x2
return
endif
do while (.not. nearlysame(x1,x2))
xm = (x1+x2)/2.
polym = poly(N,coeff,xm)
if (polym.eq.0.) then
zerofound = .true.
foundzero = xm
return
else if (poly1*polym.gt.0.) then
x1 = xm
poly1 = polym
else
x2 = xm
poly2 = polym
endif
end do
zerofound = .true.
foundzero = xm
return
end
double precision function poly(N,coeff,x)
implicit double precision (a-h,o-z)
implicit integer (i-n)
dimension coeff(0:N)
poly = coeff(0)
xpower = 1.
do i = 1,N
xpower = xpower * x
poly = poly + coeff(i)*xpower
end do
return
end
logical function nearlysame(x1,x2)
implicit double precision (a-h,o-z)
implicit integer (i-n)
xm = x1 + (x2-x1)/2.
C this is needed to avoid a too-optimizing compiler computing
C x1 - (x1 + (x2-x1)/2.) as
C x1 - x1 + (x2-x1)/2., i.e. (x2-x1)/2. ;-)
nearlysame = ((x1-xm).eq.0. .or. (x2-xm).eq.0.)
return
end
program testfindrealzeroes
implicit double precision (a-h,o-z)
implicit integer (i-n)
parameter (Nmax=100)
dimension coeff(0:Nmax), zeroes(1:Nmax)
logical EOF /.false./
do while (.not.EOF)
print *, "--------------------------------------------------"
print *, "enter degree "
READ (*, *, IOSTAT=IST) N
EOF = IST.LT.0
if(EOF) STOP "bye"
if (N.gt.Nmax) then
print *, "max degree is ",Nmax,
$ ". Enlarge ad recompile if needed"
STOP
endif
if (N.le.0) then
STOP "degree should be positive"
endif
print *, "enter ", N+1, " polynomial coefficents a_n...a0:"
READ (*, *, IOSTAT=IST) (coeff(i), i=N,0,-1)
EOF = IST.LT.0
IF (.NOT.EOF) THEN
print *, "polynomial coefficients:"
print *, (coeff(i), i=N,0,-1)
call findrealzeroes(N,coeff, Nzeroes,zeroes)
print *, "found zeroes, and value of polynomial there:"
print *, (zeroes(i), i=1,Nzeroes )
print *, (poly(N,coeff,zeroes(i)), i=1,Nzeroes)
ENDIF
end do
STOP "bye"
end
!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
--
TRu-TS
Buon vento e cieli sereni
Received on Mon Dec 31 2012 - 18:27:03 CET