[Generico] Consigli di sistemi non lineari da studiare e controllare
Salve ragazzi devo fare un progetto di dimanica e controllo non lineare. Devo controllare un sistema non lineare con vari tipologia di controllo e avevo intenzione di fare un modello con il DeepMPC. Avete da consigliarmi qualche modello non lineare da controllare? E avete da passarmi qualche materiale di DeepMPC? Preferibilemte qualche modello non lineare con al massimo 3 variabili di stato? Anche script di matlab toolbox da consigliare ecc...?
Non posso portare questi modelli perchè bannati:
• brusselator,
• chemostato,
• oscillatore di Van der Pol,
• modello di FitzHugh-Nagumo,
• modello preda-predatore (Lotka-Volterra),
• levitatore magnetico.
• pendolo planare
Non posso portare questi modelli perchè bannati:
• brusselator,
• chemostato,
• oscillatore di Van der Pol,
• modello di FitzHugh-Nagumo,
• modello preda-predatore (Lotka-Volterra),
• levitatore magnetico.
• pendolo planare
Risposte
Ci sono un paio di strade percorribili.
Se al variare dei parametri in un range assegnato il punto di equilibrio (o di funzionamento) non varia troppo, allora si possono linearizzare le equazioni attorno al punto di equilibrio nominale e quindi trovare come tale punto varia parametricamente.
Qualora ci sia una variabile principale che contraddistingue il punto di funzionamento, si può assegnare ai restanti parametri dei valori nominali, trovare la curva che descrive il punto di funzionamento al variare del parametro principale e quindi determinare le curve di correzione per tener conto degli altri parametri.
Ad esempio, se ammettiamo che il parametro principale sia T e che le altre grandezze, eccetto C_A, possano vedersi come variazioni rispetto a dei valori nominali, si può calcolare la curva C_A=C_A(T) per tali valori nominali e poi, fissandosi su un punto a T (e quindi C_A) assegnato, vedere i delta di C_A al variare degli altri parametri presi singolarmente. Si otterrano così n curve di correzione e se le variazioni sono piccole basterà sommare gli effetti di ciascuna curva.
Però hai proprio necessità di fare questa valutazione?
Se al variare dei parametri in un range assegnato il punto di equilibrio (o di funzionamento) non varia troppo, allora si possono linearizzare le equazioni attorno al punto di equilibrio nominale e quindi trovare come tale punto varia parametricamente.
Qualora ci sia una variabile principale che contraddistingue il punto di funzionamento, si può assegnare ai restanti parametri dei valori nominali, trovare la curva che descrive il punto di funzionamento al variare del parametro principale e quindi determinare le curve di correzione per tener conto degli altri parametri.
Ad esempio, se ammettiamo che il parametro principale sia T e che le altre grandezze, eccetto C_A, possano vedersi come variazioni rispetto a dei valori nominali, si può calcolare la curva C_A=C_A(T) per tali valori nominali e poi, fissandosi su un punto a T (e quindi C_A) assegnato, vedere i delta di C_A al variare degli altri parametri presi singolarmente. Si otterrano così n curve di correzione e se le variazioni sono piccole basterà sommare gli effetti di ciascuna curva.
Però hai proprio necessità di fare questa valutazione?
"ingres":
Ci sono un paio di strade percorribili.
Se al variare dei parametri in un range assegnato il punto di equilibrio (o di funzionamento) non varia troppo, allora si possono linearizzare le equazioni attorno al punto di equilibrio nominale e quindi trovare come tale punto varia parametricamente.
Qualora ci sia una variabile principale che contraddistingue il punto di funzionamento, si può assegnare ai restanti parametri dei valori nominali, trovare la curva che descrive il punto di funzionamento al variare del parametro principale e quindi determinare le curve di correzione per tener conto degli altri parametri.
Ad esempio, se ammettiamo che il parametro principale sia T e che le altre grandezze, eccetto C_A, possano vedersi come variazioni rispetto a dei valori nominali, si può calcolare la curva C_A=C_A(T) per tali valori nominali e poi, fissandosi su un punto a T (e quindi C_A) assegnato, vedere i delta di C_A al variare degli altri parametri presi singolarmente. Si otterrano così n curve di correzione e se le variazioni sono piccole basterà sommare gli effetti di ciascuna curva.
Però hai proprio necessità di fare questa valutazione?
Beh si le vedo anche in latri progetti per l'esame e quindi penso che sia necessario, mo ci provo e vediamo.
Comuqnue per tenerti informato ce l'ho fatta anche con lo sviluppo in serie di taylor facendo uno sviluppo in serie di taylor di entrambe le variabili intorno ai loro punti di equilibrio, vedi se ha senso logico o invece non ha senso logico. Perchè sto sviluppo che mi è uscito mi puzza un pò.
Allora ho preso r(t) e lo sviluppo per entrambe le variabili:
r(t)=k0e−ERT(t)CA(t)
ovvero riscritta per maggiore chiarezza in funzione di x1 e x2
r(t)=k0e−ERx2(t)x1(t)
CODICE:
Questo è lo sviluppo in serie di taylor credo sia sensato vedo termini esponenziali che rimangono se faccio la derivata parziale in funzione della variabile x1 e in funzione della variabile x2.
CODICE:
[/code]
Secondo te ha senso perchè questo è il risultato che esce soluzione parametrica che è un matrice di 2x5 soluzione rispetto ad x1 è una 1x5 e quella $x_2$ lo stesso. Però è strana perchè la soluzione paramentrica è in funzione di una z che non capisco cos'è. Ne prendo una a campione
ho messo i punti altrimenti era illegibile
Allora ho preso r(t) e lo sviluppo per entrambe le variabili:
r(t)=k0e−ERT(t)CA(t)
ovvero riscritta per maggiore chiarezza in funzione di x1 e x2
r(t)=k0e−ERx2(t)x1(t)
CODICE:
clear all close all clc syms x_1 x_2 a b c d k C_af T_f equ1= (a)*(C_af-x_1)-x_1*k*exp(-(b/x_2)) ==0; equ2= (a)*(T_f-x_2)-(c)*x_1*exp(-(b/x_2))-(d)*(x_2) ==0; equE= x_1*k*exp(-(b/x_2)); syms R E equE= subs(equE, b, E/R) load('soluzione_x1.mat') load('soluzione_x2.mat') T_equ= taylor(equE, [x_1 x_2], [sol_x1 sol_x2]) syms z equ1= (a)*(C_af-x_1)-z ==0; equ2= (a)*(T_f-x_2)-(c)*z-(d)*(x_2) ==0; set_equ=[equ1, equ2]; equZ= subs(set_equ, z, T_equ) [sol_x1 sol_x2]=solve(set_equ,[x_1,x_2])
Questo è lo sviluppo in serie di taylor credo sia sensato vedo termini esponenziali che rimangono se faccio la derivata parziale in funzione della variabile x1 e in funzione della variabile x2.
CODICE:
9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R) + k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053) - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082)^3*((0.0000000000007642938227379998291357512212758*E^2)/R^2 - (0.00000000035272177890299932004420613223581*E)/R + (0.000018780888661162956802355423359175*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*((0.000000081390591949835131801250299071242*E)/R - (0.00000000017636088945149966002210306611791*E^2)/R^2)*(x_2 - 230.75011756565803129780229007082)^2 - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082)^4*((0.0000000000015285876454759996582715024425516*E)/R - (0.0000000000000033122142289722836809600451822844*E^2)/R^2 + (0.000018780888661162956802355423359175*E*((0.00000000000025476460757933327637858374042527*E^2)/R^2 - (0.00000000017636088945149966002210306611791*E)/R + (0.000018780888661162956802355423359175*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R - (0.000000081390591949835131801250299071242*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082)^5*((0.000000000000000014354117189256991842248954755283*E^2)/R^2 - (0.0000000000000066244284579445673619200903645688*E)/R - (0.000000081390591949835131801250299071242*E*((0.00000000000025476460757933327637858374042527*E^2)/R^2 - (0.00000000017636088945149966002210306611791*E)/R + (0.000018780888661162956802355423359175*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R + (0.00000000035272177890299932004420613223581*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R + (0.000018780888661162956802355423359175*E*((0.0000000000007642938227379998291357512212758*E)/R - (0.0000000000000011040714096574278936533483940948*E^2)/R^2 + (0.000018780888661162956802355423359175*E*((0.000000000000063691151894833319094645935106317*E^2)/R^2 - (0.000000000058786963150499886674034355372636*E)/R + (0.000018780888661162956802355423359175*E*((0.0000000033912746645764638250520957946351*E)/R - (0.0000000000029393481575249943337017177686318*E^2)/R^2))/R))/R - (0.000000081390591949835131801250299071242*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R) + (0.00018776973796794374142047734691718*E*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082))/R - 1.0*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053)*(x_2 - 230.75011756565803129780229007082)^4*((0.0000000000015285876454759996582715024425516*E)/R - (0.0000000000000033122142289722836809600451822844*E^2)/R^2 + (0.000018780888661162956802355423359175*E*((0.00000000000025476460757933327637858374042527*E^2)/R^2 - (0.00000000017636088945149966002210306611791*E)/R + (0.000018780888661162956802355423359175*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R - (0.000000081390591949835131801250299071242*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 1.0*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053)*(x_2 - 230.75011756565803129780229007082)^3*((0.0000000000007642938227379998291357512212758*E^2)/R^2 - (0.00000000035272177890299932004420613223581*E)/R + (0.000018780888661162956802355423359175*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 1.0*k*exp(-(0.0043336922665508861443061157835321*E)/R)*((0.000000081390591949835131801250299071242*E)/R - (0.00000000017636088945149966002210306611791*E^2)/R^2)*(x_1 - 9.997915506321765158317363850[code]0053)*(x_2 - 230.75011756565803129780229007082)^2 + (0.000018780888661162956802355423359175*E*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053)*(x_2 - 230.75011756565803129780229007082))/R
[/code]
Secondo te ha senso perchè questo è il risultato che esce soluzione parametrica che è un matrice di 2x5 soluzione rispetto ad x1 è una 1x5 e quella $x_2$ lo stesso. Però è strana perchè la soluzione paramentrica è in funzione di una z che non capisco cos'è. Ne prendo una a campione
a*root(- 156250...000*R^5*T_f*a*k - 19899...3*E^5*a*c*k + 1562500.....0000000*R^5*a*k*z + 15625000...0000*R^5*d*k*z - 156250....00000*R^5*T_f*a^2*exp((216....603*E)/(500...0*R)) + 156...00*R^5*a^2*z*exp((216...03*E)/(500...00*R)) - 373....000*E^5*a*c*k*z^2 + 161....000*E^5*a*c*k*z^3 -35094....00*E^5*a*c*k*z^4
ho messo i punti altrimenti era illegibile
Devo guardare a fondo quello che hai postato, ma al momento sono del tuo stesso avviso, ovvero che non sembra granchè convincente.
Dovrebbe trattarsi alla fine dello sviluppo di come variano x1=C_A e x2=T attorno al loro punto di equilibrio in funzione dei parametri, a, b, c, d, k per fissati C_af e Tf , e quindi mi aspetterei, al primo ordine, una formula del tipo:
$Delta x_1= alpha_1*Delta a + beta_1*Delta b + gamma_1*Delta c + delta_1*Delta d+ epsilon1*Delta k$
$Delta x_2= alpha_2*Delta a + beta_2*Delta b + gamma_2*Delta c + delta_2*Delta d+ epsilon2*Delta k$
e non certo una dipendenza polinomiale da z che vedo peraltro definita ad un certo punto come come
(a)*(C_af-x_1)-z ==0
ovvero
z=a*(C_af-x_1)

Dovrebbe trattarsi alla fine dello sviluppo di come variano x1=C_A e x2=T attorno al loro punto di equilibrio in funzione dei parametri, a, b, c, d, k per fissati C_af e Tf , e quindi mi aspetterei, al primo ordine, una formula del tipo:
$Delta x_1= alpha_1*Delta a + beta_1*Delta b + gamma_1*Delta c + delta_1*Delta d+ epsilon1*Delta k$
$Delta x_2= alpha_2*Delta a + beta_2*Delta b + gamma_2*Delta c + delta_2*Delta d+ epsilon2*Delta k$
e non certo una dipendenza polinomiale da z che vedo peraltro definita ad un certo punto come come
(a)*(C_af-x_1)-z ==0
ovvero
z=a*(C_af-x_1)

"ingres":
Devo guardare a fondo quello che hai postato, ma al momento sono del tuo stesso avviso, ovvero che non sembra granchè convincente.![]()
Dovrebbe trattarsi alla fine dello sviluppo di come variano x1=C_A e x2=T attorno al loro punto di equilibrio in funzione dei parametri, a, b, c, d, k per fissati C_af e Tf , e quindi mi aspetterei, al primo ordine, una formula del tipo:
$Delta x_1= alpha_1*Delta a + beta_1*Delta b + gamma_1*Delta c + delta_1*Delta d+ epsilon1*Delta k$
$Delta x_2= alpha_2*Delta a + beta_2*Delta b + gamma_2*Delta c + delta_2*Delta d+ epsilon2*Delta k$
e non certo una dipendenza polinomiale da z che vedo peraltro definita ad un certo punto come come
(a)*(C_af-x_1)-z ==0
ovvero
z=a*(C_af-x_1)
no l'ho modificato e l'ho addirittura tolto e mi esce comunque non capisco perchè però non ho fissato Tf e Caf ora ci provo ma non credo cambi molto. Comunque ho fissato tutti i parametri tranne l'esponente b del termine esponenziale, ovvero:
$k_0e^(-E/(Rx_2(t)))x_1(t)=k_0e^(-b/(x_2(t)))x_1(t) $
clear all close all clc % Parametri Matlab F = 1; % Tasso di flusso volumetrico m^3/h V = 1; % Volume reattore m^3 dH = 5960; % Entalpia kcal/kmol k= 34930800; % Costante non termica 1/h rho_cp = 500; % Densita moltiplicata per la capacità termica UA = 150; % Scambio termico alla superficie T_f = 300; % Temperatura del flusso di alimentazione in ingresso K C_af=10; % Concentrazione del reagente A nel flusso di alimentazione in ingresso, misurata in kmol/m^3 T_c = 292; a=F/V; c=dH*k/(rho_cp); d=UA/(V*rho_cp); syms x_1 x_2 b equ1= (a)*(C_af-x_1)-x_1*k*exp(-(b/x_2)) ==0; equ2= (a)*(T_f-x_2)-(c)*x_1*exp(-(b/x_2))-(d)*(x_2) ==0; equE= x_1*k*exp(-(b/x_2)); load('soluzione_x1.mat') load('soluzione_x2.mat') T_equ= taylor(equE, [x_1 x_2], [sol_x1 sol_x2]) equ1a= (a)*(C_af-x_1)-T_equ ==0; equ2a= (a)*(T_f-x_2)-(c)*T_equ-(d)*(x_2) ==0; set_equ=[equ1a, equ2a]; [sol_x1 sol_x2]=solve(set_equ,[x_1,x_2])
Quello che mi esce è questo
(7.68e-413*root(1...0*b^9*z^9 - 8...0*b^8*z^9 + 2...0*b^7*z^9 - 3...0*b^6*z^9 + 2...0*b^5*z^9 - 9...0*b^4*z^9 + 1...0*b^3*z^9 - 1....0*b^2*z^9 - 2...0*b^9*z^8 + 1...0*b^8*z^8 - 5...0*b^7*z^8 + 7...0*b^6*z^8 - 5...0*b^5*z^8 + 2...0*b^4*z^8 - 4...0*b^3*z^8 + 2...0*b^2*z^8 + 2....0*b^9*z^7 - 1...00*b^8*z^7 + 4...0*b^7*z^7 - 7...0*b^6*z^7 + 5...0*b^5*z^7 - 2...0*b^4*z^7 + 4...0*b^3*z^7 - 3...0*b^2*z^7 - 1....0*b^9*z^6 + 9...0*b^8*z^6 - 2...0*b^7*z^6 + 4....0*b^6*z^6 - 3...0*b^5*z^6 + 1....0*b^4*z^6 - 3....0*b^3*z^6 + 2....0*b^5*z^5*exp((2...3*b)/5...0) + 1...0*b^4*z^5*exp((2...3*b)/5....0) - 2....0*b^3*z^5*exp((2...3*b)/5...0) + 1....0*b^2*z^5*exp((2...3*b)/5....0) - 1....0*b*z^5*exp((2...3*b)/5...0) + 3...0*b^8*z^5 + 1....0*b^7*z^5 - 1....0*b^6*z^5 + 1....0*b^5*z^5 - 6...0*b^4*z^5 + 1...0*b^3*z^5 - 1...0*b^2*z^5 + 7...0*b*z^5 + 4....0*b^5*z^4*exp((2...3*b)/5...0) - 2...0*b^4*z^4*exp((2...0) + 3...0*b^3*z^4*exp((2...0) - 1...0*b^2*z^4*exp((2...3*b)/5...0) + 2...0*b*z^4*exp((2...3*b)/5...0) - 9...0*b^9*z^4 + 7...0*b^8*z^4 - 2...0*b^7*z^4 + 4....0*b^6*z^4 - 3...0*b^5*z^4 + 1....0*b^4*z^4 - 3....0*b^3*z^4 + 3...0*b^2*z^4 -... Output truncated. Text exceeds maximum line length for Command Window display.
E penso che continui all'infinito non c'è un modo per renderla più compatta e leggibile? Tipo una riduzione ai minimi termini mettendo in evidenza tutto ciò che è possibile?
Devo vedere con Matlab per bene lo script che mi hai allegato, anche se una soluzione che parte con 7.68e-413 non mi sembra che possa funzionare.
Comunque, per altra via, ho calcolato ad es. i coefficienti del primo ordine di influenza sulla soluzione di equilibrio per a=F/V che dovrebbero essere:
$alpha_1 = 0.33$
$alpha_2 = 3.42$
Guarda se riesci a riottenere questo risultato.
Comunque, per altra via, ho calcolato ad es. i coefficienti del primo ordine di influenza sulla soluzione di equilibrio per a=F/V che dovrebbero essere:
$alpha_1 = 0.33$
$alpha_2 = 3.42$
Guarda se riesci a riottenere questo risultato.
"ingres":
Devo vedere con Matlab per bene lo script che mi hai allegato, anche se una soluzione che parte con 7.68e-413 non mi sembra che possa funzionare.
Comunque, per altra via, ho calcolato ad es. i coefficienti del primo ordine di influenza sulla soluzione di equilibrio per a=F/V che dovrebbero essere:
$alpha_1 = 0.33$
$alpha_2 = 3.42$
Guarda se riesci a riottenere questo risultato.
Aspetta non ho capito come li hai ottenuti
Semplicissimo. Ho variato F/V di una quantità $Delta$ e ho trovato i nuovi valori di equilibrio. Quindi ho fatto il rapporto tra la variazione del punto di equilibrio e il $Delta$.
"ingres":
Semplicissimo. Ho variato F/V di una quantità $Delta$ e ho trovato i nuovi valori di equilibrio. Quindi ho fatto il rapporto tra la variazione del punto di equilibrio e il $Delta$.
Ah ok alla fine stai facendo un sorta algoritmo evoluzionistico patter serch. Effittivamente sarebbe interessante creare uno codice patter search per creare un dominio di variazioni parametriche che non mi fanno scomparire il punto di equilibrio. Creando una sorta di dominio sicuro una sorta di LAS euristica (Local asymptotic stability). Sempre se non ho detto fesserie e non ha alcun senso quello che ho pensato. Però io devo cercare la formula parametrica "bruta" della serie $x=a/b$.
"MrChopin":
stai facendo un sorta algoritmo evoluzionistico patter serch
Beh veramente il calcolo fatto era molto più "casareccio"

"MrChopin":
devo cercare la formula parametrica "bruta" della serie x=a/b.
Cosa intendi per x=a/b?
"ingres":
[quote="MrChopin"]stai facendo un sorta algoritmo evoluzionistico patter serch
Beh veramente il calcolo fatto era molto più "casareccio"

"MrChopin":
devo cercare la formula parametrica "bruta" della serie x=a/b.
Cosa intendi per x=a/b?[/quote]
Hai ragione non si capisce nulla e ho scritto delle cose senza senso. Allora ti faccio un esempio pratico di un progetto di un mio amico questo è il suo sistema di equazione che rappresenta un sistema economico:
$ { ( dot(x_1)=(1-x_1/(x_(1max)))(ax_2-bx_1+c) ),( dot(x_2)=(1-x_2/(x_(2max)))(dx_1-ex_2+f)):} $
utilizza la seguente istruzione matlab:
a=2,b=3,e=2,d=4, X1M =6, X2M =4 syms x_1 x_2 [x_1 x_2 ]= solve ([(b/X1M)*( x_1 ^2) -b*x_1 -(a/X1M)*x_1*x_2+a*x_2 ==0 ,(e/X2M )*( x_2 ^2) -e*x_2 -(d/X2M )*x_1 *x_2 +d*x_1 ==0 , x_1 >=0 , x_2 >=0 , x_1 <= X1M ,x_2 <= X2M ],[x_1 ,x_2 ])
per ottenere i seguenti punti di equilibri
$ x_1=( (0), (6), (8/3), (6) ) $
$ x_2=( (0), (4), (4), (12) ) $
Quindi si ottengono i seguenti punti di equilibrio:
$ x_(eq1)=(0,0) $ Sella
$ x_(eq2)=(8/3,4) $ Attrattore
$ x_(eq3)=(6,4) $ Sella
$ x_(eq4)=(6,12) $
esclude il punto $ x_(eq4)=(6,12) $ perchè non rispetta i vincoli che $x_2$ non può essere fisicamente superiore a 6 e viene escluso.
ora si passa all'analisi di stabilità strutturale al variare del parametro $ a $ usa questo script:
syms x_1 x_2 a [x_1 x_2 ]= solve ([(b/X1M)*( x_1 ^2) -b*x_1 -(a/X1M)*x_1*x_2+a*x_2 ==0 , (e/X2M)*( x_2 ^2) -e*x_2 -d/X2M)*x_1*x_2+d*x_1 ==0 , x_1 >=0 , x_2 >=0 , x_1 <= X1M ,x_2 <= X2M ],[x_1 ,x_2 ])
e ottiene i seguenti punti di equlibrio:
$ x_1=( (0), ((4a)/3), (6) ) $
$ x_2=( (0), (4), (4) ) $
ottenendo in funzione della variazione di $ a $ :
$ x_(eq1)=(0,0) $ Attrattore per $09/2$
$ x_(eq2)=(8/3,4) $ Sella per $09/2$
$ x_(eq3)=(6,4) $ Sella per $09/2$
E così via per gli altri parametri. A me invece uscita sta cosa senza senso di prima
(7.68e-413*root(1...0*b^9*z^9 - 8...0*b^8*z^9 + 2...0*b^7*z^9 - 3...0*b^6*z^9 + 2...0*b^5*z^9 - 9...0*b^4*z^9 + 1...0*b^3*z^9 - 1....0*b^2*z^9 - 2...0*b^9*z^8 + 1...0*b^8*z^8 - 5...0*b^7*z^8 + 7...0*b^6*z^8 - 5...0*b^5*z^8 + 2...0*b^4*z^8 - 4...0*b^3*z^8 + 2...0*b^2*z^8 + 2....0*b^9*z^7 - 1...00*b^8*z^7 + 4...0*b^7*z^7 - 7...0*b^6*z^7 + 5...0*b^5*z^7 - 2...0*b^4*z^7 + 4...0*b^3*z^7 - 3...0*b^2*z^7 - 1....0*b^9*z^6 + 9...0*b^8*z^6 - 2...0*b^7*z^6 + 4....0*b^6*z^6 - 3...0*b^5*z^6 + 1....0*b^4*z^6 - 3....0*b^3*z^6 + 2....0*b^5*z^5*exp((2...3*b)/5...0) + 1...0*b^4*z^5*exp((2...3*b)/5....0) - 2....0*b^3*z^5*exp((2...3*b)/5...0) + 1....0*b^2*z^5*exp((2...3*b)/5....0) - 1....0*b*z^5*exp((2...3*b)/5...0) + 3...0*b^8*z^5 + 1....0*b^7*z^5 - 1....0*b^6*z^5 + 1....0*b^5*z^5 - 6...0*b^4*z^5 + 1...0*b^3*z^5 - 1...0*b^2*z^5 + 7...0*b*z^5 + 4....0*b^5*z^4*exp((2...3*b)/5...0) - 2...0*b^4*z^4*exp((2...0) + 3...0*b^3*z^4*exp((2...0) - 1...0*b^2*z^4*exp((2...3*b)/5...0) + 2...0*b*z^4*exp((2...3*b)/5...0) - 9...0*b^9*z^4 + 7...0*b^8*z^4 - 2...0*b^7*z^4 + 4....0*b^6*z^4 - 3...0*b^5*z^4 + 1....0*b^4*z^4 - 3....0*b^3*z^4 + 3...0*b^2*z^4 -... Output truncated. Text exceeds maximum line length for Command Window display.