Il 17/12/2012 18:30, Soviet_Mario ha scritto:
> Il 12/12/2012 01:39, Tommaso Russo, Trieste ha scritto:
>> Ma nessuno dei due!
>
> ...aggiungo qualche info sulla natura di questi polinomi, nel
> caso esista qualche suggerimento ancora più utile di quelli che già fai
> dopo.
> Si possono tutti descrivere come
>
> (a1-n1X)^n1*(a2-n2X)^n2*...*(az-nzX)^nz = K* (segue)...
> (b1-m1X)^m1*(b2-m2X)^m2*...*(bt-mtX)^mt
Insomma, se ho capito bene, la tua equazione ha la forma
P1(x) = K*P(2x)
dove P1 e P2 hanno alcune proprieta' interessanti: in particolare, zeri
multipli di molteplicita' n1...nz, m1..mt. Pero' quando porti entrambi i
membri a sinistra del segno =, ottieni un polinomio P
P(x) = P1(x) - K*P(2x) = 0
che queste proprieta' in generale non le conserva per nulla: in
particolare, non e' affatto detto che abbia anche lui zeri multipli (li
avrebbe se quelli di P1 e P2 coincidessero), anzi, non e' detto neanche
che abbia zeri... tanto vale quindi considerarlo un polinomio
"qualunque" a coefficienti reali, e cercarne gli zeri reali.
Per quanto detto sopra, quando hai paventato la "mera tangenza all'asse
X", non ti sarai per caso fasciato la testa prima di essertela rotta? La
tangenza ha di solito un ben preciso significato fisico.
Se il tuo polinomio P e' "abbastanza onesto", ed ha in generale solo
zeri reali singoli o complessi coniugati, allora il metodo che ti ho
suggerito funziona benissimo: te lo illustro meglio in un altro post. Se
invece P ha effettivamente problemi di tangenza, allora trovare i suoi
zeri multipli e' un problema niente affatto banale, che richiede
preliminarmente che tu definisca (in base a considerazioni fisiche) cosa
vuol dire, per te, "abbastanza piccolo da essere considerato nullo",
>> *Entrambi* i metodi garantiscono la convergenza solo se applicati a un
>> intervallo chiuso in cui e' certo che vi sia una e una sola radice di
>> f=0, e nel quale il prodotto f'*f" [non cambi mai di segno]
>
> urc ! ! Bella questa, le due derivate consecutive devono avere segno
> costante o che si inverte simultaneamente ???
No, scusa, ho scritto male: il prodotto deve proprio essere non nullo
ovunque (salvo che puo' annullarsi agli estremi).
Devono avere entrambe segno costante e non annullarsi: come fai a far
invertire di segno in uno stesso punto interno x_0 sia la f' che la f"?
Dovresti avere una f che da crescente diventa decrescente o viceversa, e
nello stesso punto da convessa diventa concava o viceversa...
> ... Ma che significato
> geometrico ha il segno del prodotto delle due derivate ?
Significa che nell'intervallo che contiene lo zero la f (supposta
liscia) dev'essere sempre crescente (o decrescente) e convessa (o concava).
In questo caso, e' sempre possibile determinare nell'intervallo uno
starting point da cui l'algoritmo convergera' allo zero. La ragione e'
intuitiva dalla figura:
http://commons.wikimedia.org/wiki/File:Newtonsches_N%C3%A4herungsverfahren.png?uselang=en
Se questa condizione non e' soddisfatta (e anche se e' soddisfatta, ma
scegli lo starting point dalla parte sbagliata) puo' succedere di tutto:
che l'algoritmo rimanga intrappolato nei dintorni di un minimo, come e'
successo a te:
http://commons.wikimedia.org/wiki/File:NewtonsMethodConvergenceFailure.svg?uselang=en
oppure che finisca *proprio* sul minimo, dove la f' e' 0 e l'algoritmo
termina con un errore "divide by zero";
oppure che trovi una derivata prima talmente piccola che vada a finire,
lontano, in prossimita' di *un altro* zero, o in un punto talmente
distante che solo per tornare indietro ci vogliono milioni di iterazioni:
http://en.wikipedia.org/wiki/File:Secant_method.svg
(questo e' il metodo delle secanti, ma con le tangenti succede lo stesso.)
> Mi è capitato invece pure che non convergessero a velocità ragionevoli,
> o che trovassero dei minimi locali e rimanessero intrappolate lì.
Infatti :-)
>> Immagino che n sia al piu' dell'ordine del centinaio,
>
> no no, direi della decina normalmente.
Beh, questo semplifica le cose.
>> ... sopratutto visto che delle derivate qui interessano solo gli zeri, e
>> quindi si puo' normalizzare a 1 il coefficente del termine di grado
>> massimo di ogni polinomio derivata.
>
> cioè, raccolgo il termine del massimo grado ? Non capisco in cosa mi
> semplifichi i calcoli.
Evita possibili overflow. Avevo ipotizzato che tu dovessi trattare
polinomi anche di 300-simo grado. Se derivi 300 volte un polinomio di
grado 300 senza precauzioni, finisci che hai moltiplicato il coefficente
del termine di massimo grado per 300!, che vale circa 3*10^614, al di
fuori della rappresentabilita' della doppia precisione (che va da
10^-308 a 10^308).
Se normalizzi a 1 (o anche, meglio, al valore iniziale) il termine di
grado massimo a ogni derivazione, i coefficenti rimangono vicini agli
ordini di grandezza iniziali.
Certo, se invece il tuo grado massimo e' 7 o anche 10, dato che 10! vale
circa 3 milioni, il rischio non lo corri.
>> puoi usare il teorema di Cauchy: tutti gli zeri reali di un
>> polinomio
>>
>> a_0 + a_1 x + a_2 x^2 ... + a_n x^n
>>
>> sono compresi nell'intervallo [-a, a] con
>>
>> a = 1 + max (| a_i/a_n |) , i=/=n.
>
> uhm ... qui mi sono perso completamente.
> che valore assume la i ?
> è un contatore iterativo (0<=i<n) con cui scandire tutti i rapporti, e
> di cui selezionare il massimo ?
Si', la seconda che hai detto. In altre parole, se normalizzi a 1 il
coefficiente del termine di grado massimo, a e' il massimo valore
assoluto degli altri coefficienti, piu' 1.
Insomma,
real function CauchyRadius(N,a)
C note: a_0 = a(0) is known term,
C a_i = a(i) is coefficient of x^i
dimension a(0:N)
CauchyRadius = 1.
do i=0,N-1
CauchyRadius = max (CauchyRadius, abs(a(i)/a(N))+1.)
end do
return
end
http://captainblack.wordpress.com/2009/03/08/cauchys-upper-bound-for-the-roots-of-a-polynomial/
>> Parti dalla derivata n-1esima, trova il suo zero x_0 (c'e' di sicuro,
>> altrimenti il polinomio non e' di grado n ma inferiore), e verifica se
>> cade all'inteno di [lower,upper]. Se cade al di fuori, mantieni un solo
>> intervallo
>
> cioè mantengo inalterato l'intervallo che avevo predeterminato ? Scusa
> le domande banali, ma non ho chiaramente intuito dove si debba andare a
> parare.
Dunque, cerco di rispiegare meglio....
[CUT vari tentativi]
niente da fare, un algoritmo e' piu facile scriverlo che DEscriverlo.
Ho avuto un po' di tempi morti durante le feste e ho buttato giu' la
routine.
Te la allego e descrivo in un altro post.
>> (suggerimento: per calcolare
>> i valori dei polinomi usare l'algoritmo di Horner, o almeno conservare
>> x^n per il calcolo di x^(n+1), il che rispetto ad Horner si limita a
>> raddoppiare il tempo necessario).
>
> Sin al fatto di generare iterativamente la potenza pian piano c'ero arrivato per intuizione, ma non ad inglobare la costante del grado inferiore ad ogni step ! Che figata ... grazie della segnalazione !
Attenzione pero', non innamorarti di un metodo elegante quando non ne
hai bisogno! L'algoritmo che salva ad ogni passo x^n-1 per usarlo nel
calcolo di x^n:
function poly(N,coeff,x)
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
rispetto all'algoritmo di Horner, fa il doppio di moltiplicazioni; ma la
sua correttezza (almeno nel limite della rappresentabilita' dei
risultati intermedi) risulta evidente a prima vista. In un programma
come quello che ti serve, il calcolo dei polinomi richiede si' e no 10
ms di cpu: se per ridurli a 5 ti incasini con la scrittura di un
algoritmo meno evidente, rischi di perdere molto piu' tempo, in debug,
di quello che guadagnerai in tutta la vita. Ricordati i suggerimenti di
"The elements of programming style" di Kernighan e Plauger
* "MAKE IT RIGHT BEFORE YOU MAKE IT FASTER"
* "KEEP IT RIGHT BEFORE YOU MAKE IT FASTER"
* "MAKE IT CLEAR BEFORE YOU MAKE IT FASTER"
* "DON'T SACRIFICE CLEARITY FOR SMALL GAINS IN 'EFFICENCY'"
L'algoritmo di Horner lascialo a chi deve calcolare polinomi per 48 ore:
a lui si' che vale la pena di lavorare con un po' di attenzione ridurre
il tempo di esecuzione a 24.
> P.S. non so se alla fine lo farò, perché in fondo è uno sfizio, e non so
> manco se e quanto ti possa interessare (presumo una epsilon piccola a
> piacere), ma quando/se farò funzionare al meglio questo programma
> "didattico" di risoluzione di multi equilibri, vorrei vedere se per caso
> su "Chimica nella Scuola" accettano software oltre ad articoli di solo
> testo, e se si, ovviamente ti citerei volentieri come coautore e
> consulente matematico.
Ma figurati, per un suggerimento sul metodo? In questi casi si usa al
piu' un ringraziamento...
> Ma non sono sicurissimo di perdere interesse (già scarso) a questo
> aspetto. 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.
Ma allora ti chiederei di fare qualche prova con la routine che ti mando
a parte, tanto per verificare di non esserti fasciato la testa per
niente... se poi invece effettivamente P(x) e' molto malcondizionato, ne
riparliamo.
Ciao
--
TRu-TS
Buon vento e cieli sereni
Received on Mon Dec 31 2012 - 15:37:34 CET
