[Proposto] Risolvere numericamente un sistema non lineare

Messaggioda feddy » 04/07/2018, 17:46

Ho trovato questo esercizietto carino per chi si approccia alla risoluzione di sistemi non lineari con metodi numerici:

Ex:
Si risolva1 il seguente sistema non lineare:
\( \begin{cases} \int_{x}^{y} e^{-t^2}dt =1 \\ x^2 + y^2 =1 \end{cases} \)

usando un metodo di tipo Newton


Hint "computazionale":
Testo nascosto, fai click qui per vederlo
utilizzare un'approssimazione dell'integrale attraverso i comandi disponibili, i.e. trapz o quad in Octave, in Python non saprei

Note

  1. implementandolo in MatLab/Octave/Python/quello che volete oppure anche in modo "formale"
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2096 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: [Proposto] Risolvere numericamente un sistema non lineare

Messaggioda feddy » 06/07/2018, 15:15

Se mai dovesse servire a qualcuno, posto la soluzione in spoiler :)

Testo nascosto, fai click qui per vederlo
L'iterazione funzionale che si richiede di risolvere è la solita (la scrivo con i vettori solo la prima volta, poi è sottinteso): $F(\vec{x_k}) \vec{x_{k+1}}=-JF(\vec{x_{k}})$, dove ${x_{k}}_{k \in NN}$ è la successione a valori vettoriali delle "soluzioni di tentativo" che ha come primo termine, al solito, il guess iniziale $x_0$. Qui $JF_{i,j}(x)= \frac{d F_{i}(x)}{d x_{j}}$ indica la matrice jacobiana di $F$. Ovviamente $F: RR^2 \times RR^2 \rarr RR^2$.

La funzione $F=F(x,y)$, esplicitamente, è: $F(x,y)=[\int_{x}^{y} e^{-t^2}dt -1 , \quad x^2+y^2 -1]$ e va azzerata.
Calcoliamone dunque lo Jacobiano: ovviamente l'ultima riga è data dai termini $2x$ e $2y$. Per la prima riga invece va ricordato il teorema fondamentale del calcolo integrale e la nota proprietà $\int_{a}^{b} f= - \int_{b}^{a} f$. Dunque lo jacobiano è
\( JF(x,y)= \begin{bmatrix} -e^{-x^2} & e^{-y^2} \\ 2x & 2y \end{bmatrix} \)
. Una rapida occhiata mostra subito che la matrice è singolare per la coppia $(x,y)= (0,0)$, che comunque non è soluzione del sistema, visto che non soddisfa nemmeno una delle due equazioni. Pertanto un buon guess iniziale non sarà certamente $x_0=[0,0]$.

A livello computazionale, per definire $F$ serve utilizzare, come da hint, un'approssimazione dell'integrale definito (per esempio, mediante trapezi o altre formule di quadratura).

A livello pratico questo si può fare col comando
Codice:
F=@(x) F=@(x) [-1 + quad(@(t) exp(-t^2),x(1),x(2)) ; x(1)^2+x(2)^2-1];

Mentre lo Jacobiano si scrive
Codice:
JF=@(x) [-exp(-x(1)^2),exp(-x(2)^2);2*x(1),2*x(2)];


A questo punto è sufficiente usare un metodo di tipo Newton (ma anche newton inesatto va bene), e prendiamo come guess iniziale $x_0=(1,1)$.

Codice:
x0=[1;1];
res=-JF(x0)\F(x0);
tol=1e-16; %tolleranza
maxiter=100;
iter=0;
while(norm(res,inf)>tol && iter<maxiter)
   x0+=res;
   res=-JF(x0)\F(x0);
   iter++;
   err(iter)=norm(res,inf); %tengo traccia dell'errore
end


trovando la soluzione (mi fermo alla 5 cifra decimale)
$x=( -0.27462,0.96155)$
che si verifica essere tale.

Un plot in scala semilogaritmica (su y, visti gli ordini di grandezza) mostra come descresce l'errore
Immagine
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2099 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA


Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite