[Matlab] Aitken

[Xemitron]

Salve, avrei implementato il metodo di accelerazione di aitken con il seguente codice:

function[x,i] = aitken(x0,f,df,tol,max)
    
    f0=feval(f,x0);
    d=feval(df,x0);
    x1=x0-(f0/d);
    f1=feval(f,x1);   
    d=feval(df,x1);
    x2=x1-(f1/d);
    
    for i=1:max
        x0=(x1*x1-x0*x2)/(2*x1-x2-x0);
        if abs(x0-x2)<=tol*(1+abs(x2))
            break
        end
        f0=feval(f,x0);
        d=feval(df,x0);
        x1=x0-(f0/d);
        f1=feval(f,x1);   
        d=feval(df,x1);
        x2=x1-(f1/d);
    end
    x=x2;

Nell'eseguirla sulla funzione \((x-1)^2 e^{x-1}\) con criterio d'arresto \( |x_{i+1} − x_i | ≤ tol · (1 + |x_i|)\) e con tolleranze 1e-3, 1e-6, 1e-9 ottengo rispettivamente come numero di iterazioni richieste 4,5,5. Tuttavia con una tolleranza 1e-12 non ottengo nessun risultato, fermandosi al limite massimo di iterazioni. Qualcuno potrebbe aiutarmi a capire il perché? Grazie delle eventuali risposte.

[Quinzio]

Oltre a fare il copia-incolla del codice Matlab, vuoi essere cosi' gentile da spiegarci qual e' la sequenza che vorresti trattare, cosa cosa sono "f" e "df" che si vedono nel codice, fare un esempio numerico in modo che rispondere sia un po' piu' facile ? Grazie

[feddy]

@Quinzio

f

è la funzione di cui vuoi trovare lo zero, mentre

df

la sua derivata.

@Xemitron Non ho tempo per fare debugging ora, ma sono certo che potrai trovare in rete dei file MatLab che implementano Aitken

[Xemitron]

"Quinzio":Oltre a fare il copia-incolla del codice Matlab, vuoi essere cosi' gentile da spiegarci qual e' la sequenza che vorresti trattare, cosa cosa sono "f" e "df" che si vedono nel codice, fare un esempio numerico in modo che rispondere sia un po' piu' facile ? Grazie

Allora si, come detto già da feddy f e df sono rispettivamente funzione e derivata della funzione, x0 il punto iniziale, tol è la tolleranza da utilizzare e max il numero massimo di iterazioni, e la funzione implementa il metodo di accelerazione di aitken. L'ho testato sulla funzione che ho scritto nel primo post in questo modo:

f =

     Inline function:
     f(x) = ((x-1)^3)*exp(x-1)

>> df=fcnchk('exp(x-1)*(3*((x-1)^2)+(x-1)^3)')

df =

     Inline function:
     df(x) = exp(x-1)*(3*((x-1)^2)+(x-1)^3)

>> [x,i]=aitken(0,f,df,1e-3,100)

x =

     1.0000


i =

     4

>> [x,i]=aitken(0,f,df,1e-6,100)

x =

     1


i =

     5

>> [x,i]=aitken(0,f,df,1e-9,100)

x =

    1


i =

     5

>> [x,i]=aitken(0,f,df,1e-12,100)

x =

   NaN


i =

   100

Ecco, con i primi tre valori di tolleranza mi torna ma con l'ultimo no.

"feddy":@Quinzio

f

è la funzione di cui vuoi trovare lo zero, mentre

df

la sua derivata.

@Xemitron Non ho tempo per fare debugging ora, ma sono certo che potrai trovare in rete dei file MatLab che implementano Aitken

Ho provato ad usare versione trovate in rete ma il risultato è lo stesso. C'è un errore nel criterio d'arresto, dal momento che è l'unica cosa che cambia?

[feddy]

Quel NaN è il risultato di una division per $0$, che proviene dalla valutazione dello jacobiano in un punto in cui la funzione è piatta. Se guardi il grafico della funzione, ti rendi conto che vicino a $x=1$, cioè il tuo zero, la funzione è piatta. In pratica, lì si annulla anche lo jacobiano. Alla 5 iterazione, hai già che $x_0=1$ (esattamente), mentre

d

è $0$. Il tuo zero l'hai trovato, eccome, ma quando poi il programma arriva alla riga successiva alla valutazione dello jacobiano, cioè alla riga

x1=x0-(f0/d);

, **divide per 0**, da cui il primo NaN. Il resto ovviamente è da buttare.

Quindi, in sostanza, non hai prestato attenzione all'eventualità che avvenga una divisione per $0$.

[feddy]

Ho notato che ho scritto jacobiano, avevo in mente il caso di un sistema di equazioni. In questo caso intendo la classica derivata ovviamente.