[Proposto] Risolvere numericamente un sistema non lineare

[feddy]

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

Ex:

Si risolva[nota]implementandolo in MatLab/Octave/Python/quello che volete oppure anche in modo "formale"[/nota] 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":

utilizzare un'approssimazione dell'integrale attraverso i comandi disponibili, i.e. trapz o quad in Octave, in Python non saprei

[feddy]

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

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

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

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)$.

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