Test di Hurwitz in Maxima

[DarwinNE]

Sono sempre rimasto abbastanza intrigato dai problemi di teoria delle reti elettriche. Spesso sono intimamente legati a problemi non banali di natura matematica, ed è una matematica applicata molto interessante.

Nei circuiti a parametri concentrati con un numero finito di elementi, le funzioni di rete sono del tipo:

dove e sono dei polinomi nella variabile complessa . Una cosa importante da sapere è se gli zeri dei due polinomi hanno parte reale negativa oppure no, dato che queste informazioni sono di solito legate alla stabilità del sistema, oppure al fatto che la funzione è a fase minima oppure no. In campo reale le cose sono semplici da verificare utilizzando la regola di Cartesio, ma in campo complesso la situazione è molto più spinosa.

Diventa quindi importante verificare questa proprietà per un generico polinomio, se possibile senza calcolare esplicitamente gli zeri del polinomio, dato che questo potrebbe essere un problema non banale per un grado elevato del polinomio da studiare. Esiste un test algebrico chiamato prova di Hurwitz che permette di determinare se questa proprietà è verificata oppure no, senza calcolare esplicitamente gli zeri. Un polinomio è hurwitziano se tutti i suoi zeri hanno parte reale strettamente negativa, se ce n'é qualcuno immaginario è hurwitizano in senso lato.

Qualche anno fa, avevo messo a punto un piccolo script Maxima che permette di verificare se un polinomio è hurwitziano oppure no.
Ne ho già parlato sul mio sito:
http://davbucci.chez-alice.fr/index.php ... aliano#fun
ma mi sono detto che sarebbe stato carino proporlo nuovamente qui, per stimolare una discussione e perché magari potrebbe essere utile a qualcuno.

Codice: Seleziona tutto

/* Perform the Hurwitz test on the M(p) polynomials using Maxima
   Version 1.0
   
   
   Written by Davide Bucci, March-April 2007
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

/* Qualche esempio di polinomio hurwitziano e non:
   -----------------------------------------------
   56*p^5+7*p^4+161*p^3+18*p^2+98*p+8,  Hurwitziano in senso lato
   p^3+2*p^2+4*p+3,                     Strettamente hurwitziano   
   p^5+2*p^4+2*p^3+4*p^2+11*p+10,       Non hurwitziano
   p^3+2*p^2+4*p+9,                     Non hurwitziano 
*/

hurwitz(poly,var):=block ([M,p],

  /* Definiamo una funzione che ci permettera' di lavorare sul nostro
     polinomio in maniera generale, utilizzando la variabile var */
     
  define(M(var),poly),
 
 
  if(debugmode) then print(M(p)),
  print("Prova di Hurwitz per un polinomio"),
  /* Otteniamo parte pari e dispari del polinomio */
  Pa(p):=expand(1/2*(M(p)+M(-p))),
  Di(p):=expand(1/2*(M(p)-M(-p))),

  /* Abbiamo bisogno di ottenere il grado piu' elevato al numeratore */
  pg:hipow(Pa(p), p),
  dg:hipow(Di(p), p),
  num:Di(p),
  den:Pa(p),
  nm:dg,
   
  if (pg>dg) then
     (num:Pa(p), den:Di(p), nm:pg),
     
 
  res:divide(num,den,p),
 
  /* Adesso possiamo effettuare le divisioni successive */
 
  hasdivided: (second(res)#0), 
  print(hasdivided),
 
  /* This list will contain the results of the successive divisions */
  reslist:[],
  ishurwitz:true,
 
  print("Calcolo delle divisioni successive"),
  for i:0 thru nm while (hasdivided and ishurwitz) do (
    if(debugmode) then (
      print("nuova iterazione"),
       print(i),
       print("divide"),
       print(num),
       print("by"),
       print(den)
     ),
     res:divide(num,den,p),
     if(debugmode) then (
       print("resulting"),
       print(res)
     ),
     
    hasdivided: (second(res)#0),
 
     reslist:append(reslist,[first(res)]),
 
    num:expand(den),
     den:second(res),
     
     /* Si controlla se il risultato della divisione e' un monomio di grado 1
       con coefficiente positivo */
     
     if (hipow(first(res),p)>1) then (
        ishurwitz:false
     ),
     if (coeff(first(res),p,1)<0) then (
        ishurwitz:false
     )
   
  ),
 
  if (length(reslist)=nm and ishurwitz) then
     print("Polinomio strettamente hurwitziano")
  elseif (ishurwitz) then
    print("Polinomio hurwitziano in senso lato")
  else
    print("Polinomio non hurwiztiano"),
 
  print(reslist),
 
  reslist
)$



Una volta salvato lo script in un file, per esempio hurwitz.max, lo si può utilizzare all'interno di Maxima:

Codice: Seleziona tutto

(%i5) load("hurwitz.max");
(%o5)                             Hurwitz.max
(%i6) hurwitz(56*p^5+7*p^4+161*p^3+18*p^2+98*p+8,p);
Prova di Hurwitz per un polinomio
    3
17 p  + 34 p # 0
Calcolo delle divisioni successive
Polinomio hurwitziano in senso lato
      7 p  17 p
[8 p, ---, ----]
      17    4
                                     7 p  17 p
(%o6)                          [8 p, ---, ----]
                                     17    4
(%i7)


La funzione hurwitz(N,p) richiede il polinomio N nella variabile p ed effettua il test, che è basato su una serie di divisioni successive fra la parte pari e la parte dispari del poliniomio di partenza. Il programma restituisce una lista contenente i risultati delle divisioni.
Una cosa che mi piacerebbe fare una volta o l'altra è preparare altri script di questo tipo che applichino alcune procedure di sintesi empirica per le reti lineari passive. Che ne pensate?

[RenzoDF]

Si, un interessante script DarwinNE, ma direi che una volta disponibile Maxima, forse con una fattorizzazione seguita dalla

------->>>> regola di Cartesio

si fa prima! ... normalmente poi con Maxima si riescono anche a trovare direttamente le radici

Se facciamo per esempio riferimento al polinomio usato in questo documento

RouthHurwitzAnalysis.pdf

(21.96 KiB) Scaricato 158 volte

, dove è riportato il procedimento manuale,
potremo risolvere con uno scanmap('factor, f(s)); e poi con regola di Cartesio ... e concludere dicendo che: ci sarà una radice doppia nell'origine, 2 coppie di radici complesse coniugate a parte reale nulla e, vista la doppia permanenza dei segni dell'ultimo fattore, ancora una coppia di radici coniugate con parte reale negativa

RH.gif (19.05 KiB) Visto 2057 volte

ed inoltre, volendo, in questo caso particolare, determinare le soluzioni con un semplice solve.

Sarebbe invece interessante, caro Davide, la tua idea di sintesi automatica con Cauer e Foster a partire dall' ammettenza /impedenza complessa, magari con un bell'articolo

[DarwinNE]

Si, un interessante script DarwinNE, ma direi che una volta disponibile Maxima, forse con una fattorizzazione seguita dalla

------->>>> regola di Cartesio

si fa prima! ... normalmente poi con Maxima si riescono anche a trovare direttamente le radici

E' vero! Oggigiorno c'è da dire che ci sono strumenti talmente potenti che trovare delle radici di un polinomio di grado non banale è relativamente semplice. Credo che nei casi che si presentano nella pratica sia sempre possibile attaccarli con algoritmi numerici, ma in qualche situazione c'è da stare attenti. Prova a fare qualche prova con il polinomio seguente, che pure è "solo" di grado 8:




L'unico modo di affrontare il problema è utilizzare la funzione "allroots" che calcola le radici numericamente ed è una cosa molto più pesante da un punto di vista computazionale del semplice test di Hurwitz.
Nell'esempio qui sotto, invece, factor non funziona e allroots mi pesca addirittura delle radici reali.

Codice: Seleziona tutto

a:expand((x+1+1e-3*%i)*(x+1-1e-3*%i));
b:expand((x+2+1e-3*%i)*(x+2-1e-3*%i));
c:expand((x+3+1e-3*%i)*(x+3-1e-3*%i));
d:expand((x+4+1e-3*%i)*(x+4-1e-3*%i));
e:expand((x+5+1e-3*%i)*(x+5-1e-3*%i));
f:expand((x+6+1e-3*%i)*(x+6-1e-3*%i));
g:expand((x+7+1e-3*%i)*(x+7-1e-3*%i));
h:expand((x+8+1e-3*%i)*(x+8-1e-3*%i));
p:expand(a*b*c*d*e*f*g*h);
allroots(p);


C'è da dire che gli esempi di cui sopra (il primo è expand(a*b*c*d)) sono delle varianti dei polinomi di Wilkinson, che sono ben conosciuti per essere molto ostici da trattare numericamente.
Per i problemi di sintesi empirica, provo a rifletterci.

[RenzoDF]

Si d'accordo ma dal punto di vista ingegnetistico "usare" coefficienti a 15 cifre significative è pura fantascienza !

[DarwinNE]

Si d'accordo ma dal punto di vista ingegnetistico "usare" coefficienti a 15 cifre significative è pura fantascienza !

Sono totalmente d'accordo, però se il problema è mal condizionato, invece delle 15 cifre significative dei coefficienti del polinomio, le soluzioni ne hanno al massimo una o due (vedi difatti cosa fa allroots! ).