Il 31/12/2012 15:37, Tommaso Russo, Trieste ha scritto:
> 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.
no, me la sono ammaccata, ma ho tamponato, prima di
ipotizzarla. Oddio la bisezione mi aveva causato degli stack
esauriti quando usavo la ricorsione, ma poi ho lavorato e ho
tramutato una semplicissima sintassi ricorsiva in una
complicata e macchinosa sintassi iterativa, che per� usava
solo poca memoria.
In particolare la ricorsiva consumava l'integrale delle
potenze di due, la iterativa, solo l'ultima potenza di due
Per dire, al 5 stadio di profondit�
Ricorsiva
2^0 + 2^1 + 2^2 + 2^3 + 2^4 + 2^5 = 2^6 - 1 = 63
Iterativa
2^5 = 32.
Praticamente consuma la met� di memoria (di solo
tracciamento, senza contare le variabili locali). Ed � anche
pi� veloce.
Ma cmq sia, ho scoperto una cosa sgradevole che ti racconto
dopo ... anzi subito.
Ebbene, ho cambiato completamente strategia e non ricerco
pi� nemmeno le radici del polinomio.
Perch� mi sono accorto che in un sistema non lineare, non
saprei cosa farmene delle soluzioni.
Mi spiego meglio, oltre alle VARIABILI, uno o pi�, al limite
tutti i termini COSTANTI di ciascun polinomio possono
essere, variazionalmente, ad ogni iterazione, modificati
dagli altri polinomi. Quindi le soluzioni singole non
rappresentano per forza soluzioni del sistema nel suo complesso.
Ed � vero che le soluzioni FINALI, a convergenza avvenuta,
di ogni polinomio, non rappresentano soluzioni delle
rispettive formulazioni iniziali (perch� tutte le costanti
sono state variate ad ogni iterazione).
Provo a rispiegare questo punto.
Avendo scoperto che per polinomi di grado superiore non
esistono procedure deterministiche di previsione delle
soluzioni (come quelle dell'algebra lineare stile Gauss
Jordan che mi aveva spiegato Elio un po' di tempo fa),
allora avevo implicitamente puntato su una procedura
evolutiva, variazionale.
E l'effetto collaterale di ci� � che ciascun polinomio,
anche se non varia struttura formale, pure varia una o pi�
d'una o tutte le costanti, che si coevolvono insieme.
Quindi la ricerca delle soluzioni esatte di ognuno
singolarmente l'ho abbandonata : ho fatto delle prove con
polinomi generati da radici note, e queste soluzioni non
soddisfavano il sistema completo, ma solo il padre della
soluzione.
A me invece interessava risolvere il sistema.
Ho optato per una soluzione ispirata dalla fisica
sottostante : dai concetti di compensazione degli scarti
dall'equilibrio, e la rispettiva proporzionalit� al
Logaritmo del rapporto tra il quozione attuale di reazione e
il valore della sua costante di equilibrio.
Questi valori di Ln(Q_r / K_eq) rappresentano differenze
di energia libera. Ebbene io ho ordinato la priorit� di ogni
equazione in maniera proporzionale al suo scartamento, di
modo che la sua richiesta di compensazione avesse priorit�
proporzionale. Ho anche rinormalizzato i pesi, in modo che
tutto ci� diventasse adimensionale, diciamo.
Bene, converge tutto, fa qualche bizza giusto all'inizio,
con un po' di oscillazioni, un paio, poi converge tutto.
E' algoritmicamente lento, nel senso che occorre qualche
migliaio di iterazioni (Che per� sono molto rapide, perch�
ogni computazione � il computo di una formula analitica e
non pi� una ricerca a sua volta, n� Newton Raphson n�
bisezione, niente : esegue il calcolo del polinomio come
espressione unica di costanti, il valore di partenza e la
variazione stimata).
Quindi normalmente conti sino a dieci, e visualizzazioen
grafica inclusa, completi diciamo 5000 iterazioni con gi�
notevole convergenza.
Cmq la nostra discussione mi � servita molto, anche per
superare un impasse psicologica.
Mi legger� il resto della risposta con calma, perch� al
momento ho imboccato una strada diversa di ricerca multipla
contemporanea.
Credo che l'intuizione di tipo "fisico" abbia consentito di
cercare le soluzioni nella loro direzione naturale, ossia
consente di sfruttare il fatto che ciascuna equazione in
realt� sa sempre dove vuole andare, e non deve cercare a
caso. Sa di essere troppa o troppo poca, e anche stimare
comparativamente di quanto.
A risentirci ...
A (s)proposito ... ti ho anche cercato su feisbuk, ma forse
sono stato invadente. Se � cos�, mi scuso, � stata una cosa
impulsiva, ma non mi pare che si possa ritirare una
richiesta. Ormai � nelle tue mani.
Ciao e grazie delle consulenze e della pazienza come al solito
SOVIET
>
> 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.
che figata, ottima considerazione.
>
> 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
>
>
--
1) Resistere, resistere, resistere.
2) Se tutti pagano le tasse, le tasse le pagano tutti
Soviet_Mario - (aka Gatto_Vizzato)
Received on Wed Jan 02 2013 - 15:43:26 CET
