Condizionamento Degli Zeri Di Un Polinomio

Calcolo numerico

1.

Condizionamento degli zeri di un polinomio

Approccio al problema

Il condizionamento degli zeri (ξ) al variare dei coefficienti a  del polinomio è definito come: 

ε ry

ξ k − n −1 ⋅ a n = Cp = ,  ε rx p' (ξ ) Δ

occorre quindi in primo luogo determinare le radici del polinomio.  Radici. Per prima cosa il polinomio viene tabulato per poter isolare gli zeri per via  grafica. 

  Dai grafici di p ( x ) , a sinistra con l’asse delle ordinate in scala logaritmica del modulo di y, a destra in scala lineare si nota la presenza di tre radici intorno a -5, -1, e 60,  identificate rispettivamente con ξ1, ξ2 e ξ3.  Si osserva immediatamente come la rapidità della funzione intorno alla radice ξ3  renda  il  problema  particolarmente  malcondizionato:  ad  una  piccola  variazione  sulla  x   12 corrisponde  una  grande  variazione  della  y   (l’asse  y   è  visualizzato  con  scala  10 ).  Per  la  stessa radice, nel grafico logaritmico, si osserva una variazione di segno in cui il modulo  -12 della funzione è superiore ai 10  anche se il passo della y  è di sole 0,2 unità; Il picco verso il basso avrebbe dovuto tendere verso -∞.  Dall’ingrandimento  della  funzione  intorno  all’origine  si  osserva  la  presenza  di  flessi  in  prossimità di ξ1 e di ξ2, dei quali  occorre  tener  conto  per  la  determinazione della radici, in particolare  ξ1  potrebbe  coincidere  con un flesso. 

Roberto Patrizi, 09114657

pag 1 di 8 

Calcolo numerico

2.

Condizionamento degli zeri di un polinomio

Ricerca delle radici con Matlab

Sono stati implementati 2 algoritmi in linguaggio Matlab per la ricerca delle radici:  il  metodo  di  Newton  ed  il  metodo  di  bisezione,  con  i  rispettivi  script  riportati  in  appendice. Il calcolatore utilizzato è in grado di operare con numeri in virgola mobile con  una precisione massima di 17 cifre.  Radice ξ1. La  presenza  del  flesso  (p ’ ’ ( x ) =0)  in  prossimità  della  radice  rende  difficoltosa  la  ricerca  di  un  intervallo  [a,b],  all’interno  del  quale  deve  trovarsi  la  radice  ed  in  cui  la  derivata  seconda  non  è  mai  nulla.  Le  condizioni  sovra  esposte  garantiscono  la  convergenza  del  metodo  di  Newton a partire da un estremo di  Fourier (p ( x 0 ) ·p ’ ’ ( x 0 ) >0).  Non  è  neppure  possibile  garantire che  p ’ ( x )  sia  non nulla  in un intorno della radice, pertanto  non  sono  soddisfatte  le  condizioni per la convergenza del metodo di Newton. La radice potrebbe avere molteplicità 3,  in tal caso sarebbe possibile applicare il metodo di Newton modificato, per il quale: 

x n +1 = x n − 3

p( x n )   p' ( x n )

Con il metodo di bisezione invece è sufficiente scegliere un intervallo [a,b] per il  quale risulta p (a)·p (b)<0, ciò è verificato per a=-4.6, b=-4. L’output ottenuto è:  >>bisezpoly(A,a,b,epsilon); ------------------------------------------------------------------------------n a b m p(m) ------------------------------------------------------------------------------1 -4.600000000000000 -4.000000000000000 -4.300000000000000 5.133e+000 ... ... ... ... 46 -4.333329182771303 -4.333329182771286 -4.333329182771294 -3.492e-009 47 -4.333329182771294 -4.333329182771286 -4.333329182771291 -6.694e-010 48 -4.333329182771291 -4.333329182771286 -4.333329182771289 -3.900e-009 49 -4.333329182771289 -4.333329182771286 -4.333329182771287 2.037e-010 50 -4.333329182771289 -4.333329182771287 -4.333329182771288 0.000e+000 Raggiunta radice in 50 passi>> -16

in cui epsilon è stato impostato a 1·10 . Si assumerà come radice   ξ1=-4.333329182771288.

Roberto Patrizi, 09114657

pag 2 di 8 

Calcolo numerico

Condizionamento degli zeri di un polinomio

Radice ξ2. Dal grafico intorno alla seconda  radice,  si  osserva  la  presenza  di  una  distanza tangibile tra la radice ed il nullo della  derivata  seconda.  In  questo  caso  quindi  sono valide le ipotesi:  2

p (a) · p (b) < 0     con p (x ) X C [a, b]  p ’’(x ) K0, ∀x X [a, b] 

p ’(ξ) K0  in  cui  l’intervallo  [a,b]  è  ad  esempio  a=-0.7,  b=-0.5. Dunque è sufficiente determinare un  estremo di fourier dal quale iniziare ad iterare  con  il  metodo  di  newton  per  avere  la  convergenza.  Ad  esempio  posto  x 0=-0.7,  p (x 0)*p ’ ’ (x 0)>0, dunque  x 0 è un estremo di Fourier. L’esecuzione dello script newtonzp  da luogo al seguente output:  >> newtonzp(A,0,epsilon,delta) --------------------------------------n xi p(xi) --------------------------------------1 -0.735849056603774 3.955e+004 2 -0.601542910811362 1.599e+004 3 -0.590245920932620 1.148e+003 4 -0.590158659143148 8.732e+000 5 -0.590158653920796 5.225e-004 Raggiunto limite precisione macchina in 5 passi -16

10

in cui epsilon=1·10 , delta=1·10 . L’algoritmo si arresta alla 5 iterazione, quando  p (x i ) è ancora dell’ordine di 10-4. Con il metodo di bisezione invece:  >> bisezpoly(A,a,b,,epsilon) -----------------------------------------------------------------------------n a b m p(m) -----------------------------------------------------------------------------1 -.59999999999999998 .50000000000000000 0.55000000000000004 -3.907e+003 ... ... ... ... 47 -.59015865392079669 -.59015865392079525 -.59015865392079592 1.455e-011 48 -.59015865392079592 -.59015865392079525 -.59015865392079558 -1.455e-011 49 -.59015865392079592 -.59015865392079558 -.59015865392079569 0.000e+000 Raggiunta radice in 49 passi

La  precisione  della  radice  è  maggiore,  l’algoritmo  sembra  funzionare  meglio  del  metodo di Newton. Ciò probabilmente è dovuto ad errori di arrotondamento eseguiti nel  -4 calcolo del rapporto delle funzioni  p (x )/p ’(x ) (dell’ordine di 10  ). Infatti la differenza tra  x i -p (x i )/p ’(x i ) non genera problemi dovuti all’annullamento in quanto  x i  è svariati or5 dini di grandezza superiore del rapporto. Mentre p ’(x i ) si mantiene elevato ordine di 10 ,  p (x i ) decresce fino a rendere inattendibile il rapporto tra le funzioni.  ξ2=-0.59015865392079569

Roberto Patrizi, 09114657

pag 3 di 8 

Calcolo numerico

Condizionamento degli zeri di un polinomio

Radice ξ3. Per la terza radice è facile individuare le ipotesi per la convergenza del  metodo  di  Newton,  ed  è  anche  facile  individuare  l’estremo  di  Fourier  x 0,  per  il  quale  p (x 0)*p ’ ’ (x 0)>0. Tuttavia la ripidità della funzione nell’intorno della radice, da luogo ad  un  problema  malcondizionato:  non  appena  ci  discostiamo  dalla  radice  ξ3, otteniamo  un  valore  del  polinomio  molto  elevato  in  modulo.  Osserviamo  i  dati  ottenuti  con  -16 10 l’applicazione del metodo con x0=62, epsilon=1·10 , delta=1·10 :  newtonzp(A,x0,epsilon,delta) ---------------------------------------n xi p(xi) ---------------------------------------1 61.157300198834456 1.721e+012 2 61.085638246829994 1.254e+011 3 61.085147905914937 8.461e+008 4 61.085147883085845 3.939e+004 5 61.085147883085845 3.466e-005 Raggiunto limite precisione macchina in 5 passi

In  soli  5  passi  viene  raggiunto  il  limite  di  precisione  macchina  sulla  x,  quando  il  -5 polinomio è dell’ordine di 10 . Forzando l’algoritmo a continuare, si ottengono tutte righe analoghe alla quinta. Notevole il guadagno di 9 ordini di grandezza sulla y  in una sola iterazione, che implica che la funzione è ben linearizzabile intorno alla radice.  Occorre  osservare  che  né  con  il  metodo  di  bisezione,  né  con  altri  metodi  non  è  possibile ottenere una precisione maggiore: Infatti osserviamo i dati ottenuti con il metodo di bisezione:  >> bisezpoly(A,61,61.1,epsilon) ------------------------------------------------------------------------------n a b m p(m) ------------------------------------------------------------------------------1 61.000000000000000 61.100000000000001 61.049999999999997 -6.044e+010 43 61.085147883085831 61.085147883085853 61.085147883085838 -1.289e-002 44 61.085147883085838 61.085147883085853 61.085147883085845 3.466e-005 45 61.085147883085838 61.085147883085845 61.085147883085838 -1.289e-002 46 61.085147883085838 61.085147883085845 61.085147883085838 -1.289e-002

  Il metodo giunge alla stessa identica radice del metodo di newton in 44 iterazioni,  che tra l’altro è la migliore che si riesce ad ottenere, come visibile dai risultati.  Alla riga 45 si osserva che la media (a+b)/2, coincide con il valore a, risultato ovviamente privo di senso matematico. Ciò è dovuto all’errore di rappresentazione dei numeri del calcolatore utilizzato. Infatti eseguendo in matlab i comandi:  >> xi=61.085147883085838; %il valore di a alla riga 44 >> xi+eps(xi) ans = 61.085147883085845

otteniamo che all’iterazione 45 il numero successivo ad a che è possibile rappresentare è proprio b: non esistono numeri intermedi per il calcolatore. Ciò arresta di fatto  l’algoritmo, e impedisce di ottenere una soluzione più accurata, se non facendo uso di un  calcolatore in grado di eseguire calcoli in virgola mobile con un numero maggiore di cifre.  ξ3=61.085147883085845 

Roberto Patrizi, 09114657

pag 4 di 8 

Calcolo numerico

3.

Condizionamento degli zeri di un polinomio

Coefficienti di condizionamento

Perturbando il termine noto, il coefficiente di condizionamento che si ottiene ha  la forma:  Δ

Cp =

ε ry ε rx

=

an   p' (ξ )

che per le radici trovate assume il valore rispettivamente di   Cp1 = 5.393094715224958e+009  Cp2 = 0.39524133495742  Cp3 = 2.291931500642379e-008  Intuitivamente si osserva come al variare del termine noto, il polinomio venga traslato orizzontalmente, con il conseguente spostamento delle radici. Questo spostamento  è tanto maggiore, quanto più il polinomio interseca l’asse delle ascisse orizzontalmente.  Ciò si rispecchia nei coefficienti di condizionamento: per la radice 1 a tangenza orizzontale  il  Cp1  è  notevolmente  elevato,  al  contrario  per  la  radice  3  data  da  un’intersezione  quasi verticale del polinomio con l’asse delle ascisse, si ha un Cp3 particolarmente basso,  come era ragionevole aspettarsi.  Ciò significa che per la terza radice non è critico il termine noto, anche per variazioni piuttosto ampie dello stesso non si hanno consistenti scostamenti della radice. È invece critica la determinazione della stessa, per la quale, infatti, su ha un indice di condi18 zionamento della radice di circa 3·10 .  Per la prima radice si ha il discorso inverso, mentre la seconda è in un certo qual  modo più “equilibrata”, le difficoltà della determinazione della radice e la sua sensibilità a  variazioni del termine noto si equivalgono, e non comportano particolari problemi.  Perturbando  il  termine di quinto  grado,  per il  quale  k -n -1=4,  i  valori  dei coefficienti di condizionamento sono dati da: 

ξ 4 ⋅ a2 Cp =   p' (ξ ) Si ottengono quindi i valori  Cp1 = -9.629300430614109e+011  Cp2 = 0.02427769436537  Cp3 = -0.16158986377900  Anche in questo caso le radici più sensibili sono quelle in cui il polinomio interseca  l’asse delle ascisse con pendenza più bassa, vale a dire ξ1 in primo luogo (a tangenza orizzontale) e in misura minore ξ2. In generale le radici multiple, per le quali si annulla anche  la derivata prima sono malcondizionate rispetto a perturbazioni dei coefficienti. 

Roberto Patrizi, 09114657

pag 5 di 8 

Calcolo numerico

4.

Condizionamento degli zeri di un polinomio

Codice degli algoritmi

Poly.dat 27 -1296 -20025 -92663 -146320 -35555 -53742 -39546

function A1=deriva(A) % % A1=deriva(A) % %Se A è una matrice dei coefficienti di un polinomio, la funzione %restituisce la matrice A1 dei coefficienti della derivata del polinomio % n=length(A)-1; %è il grado del polinomio (il -1 serve per la presenza del termine noto) for i=1:n A1(i)=A(i)*(n+1-i); end;

  function Y=polycalc(A,X) % % Y=polycalc(A,X) % %Funzione per il calcolo del valore assunto da un polinomio in un punto %(o vettore di punti X). A è la matrice dei coefficienti del polinomio % n=length(A); Y=X*0; %Inizializzazione di Y for i=1:n Y=Y+A(i)*(X.^(n-i)); end;

    function m=bisezpoly(A,a,b,epsilon) % % m=bisezpoly(A,a,b,epsilon) % %Calcolo gli zeri di un polinomio tramite il metodo di Bisezione. A è la %matrice dei coefficienti del polinomio, a e b sono gli estremi %dell'intervallo al quale è applicato il metodo, epsilon è la precisione %che si vuole ottenere sulla determinazione della radice m % %Visualizzazione dei passaggi s='---------------------------------------------------------------' ; disp(s) fprintf(' n\t\t\t a\t\t\t\t\t b\t\t\t\t\t m\t\t\t\t p(m)\n') disp(s) n=1; %inizializzazione N=ceil((log(b-a)-log(epsilon))/log(2)); %numero di iterazioni richieste fa=polycalc(A,a); fb=polycalc(A,b);

Roberto Patrizi, 09114657

pag 6 di 8 

Calcolo numerico

Condizionamento degli zeri di un polinomio

if sign(fa)==sign(fb) %verifica intervallo fprintf('L''intervallo non è valido!' ) else while n
    function xi=newtonzp(A,x0,epsilon,delta) % % xi=newtonzp(A,x0,epsilon,delta) % %Calcolo gli zeri di un polinomio tramite il metodo di Newton. A è la %matrice dei coefficienti di un polinomio, x0 è il punto iniziale dal quale %far partire l'algoritmo, xi è la radice alla quale è giunto il metodo, %epsilon e delta sono le precisioni che si vuole raggiungere, %rispettivamente sulla x e sulla f(x) %Inizializzazione A1=deriva(A); xi=x0; x=xi; %Visualizzazione dei passaggi s='---------------------------------------'; disp(s) fprintf(' n\t\t\t xi\t\t\t\t\t p(xi)\n') disp(s)

Roberto Patrizi, 09114657

pag 7 di 8 

Calcolo numerico

Condizionamento degli zeri di un polinomio

%parte centrale dell'algoritmo: for n=1:100 %100=numero massimo di iterazioni p=polycalc(A,x); p1=polycalc(A1,x); r=p/p1; %controllo sul rapporto delle funzioni if abs(r) < eps fprintf('\nRaggiunto limite precisione macchina in %d passi',n) return end %calcolo radice successiva xi=x-(p/p1); fprintf('%2d\t %17.15f\t\t %8.3e\n' ,n,xi,p) %controllo precisione %if (abs(xi-x)<epsilon)|(abs(p)<delta) if (abs(p)<delta) fprintf('\nRaggiunta precisione richiestain %d passi',n) break ; end x=xi; end;

 

Roberto Patrizi, 09114657

pag 8 di 8