Generazione di una mesh strutturata su un cerchio

Si vuole rappresentare graficamente una funzione $z=f(x,y),$ sul cerchio di $\mathbb{R}^2$ di centro l'origine e raggio 1.
A tale proposito si vuole generare una mesh strutturata (ovvero una opportuna trasformazione di una griglia cartesiana), in modo da poter utilizzare le function di MATLAB mesh, surf, meshc,....
In alto la mesh cartesiana sul quadrato $(0,1)^2$ (20 x 20 punti), in basso la trasformazione della mesh cartesiana sul cerchio (sempre 20 x 20 punti).
Ogni punto $P_{ij}$ (con $i=1,...,nx$ e $j=1,...,ny$ ) del quadrato viene trasformato nel corrispondente $Q_{ij}$ del cerchio (sempre con $i=1,...,nx$ e $j=1,...,ny$ ).

Per generare tale griglia si devono risolvere due sistemi lineari $A{\bf x}={\bf b1},\ A{\bf y}={\bf b2}$ , in cui i vettori incogniti ${\bf x}$ e ${\bf y}$ contengono le coordinate (rispettivamente ascisse e ordinate) dei punti della griglia.
Indicando con $nx$ e $ny$ il numero di punti lungo la direzione $x$ e $y$ , rispettivamente, i sistemi lineari da risolvere hanno dimensione $n=nx*ny$ . Il vettore x, soluzione di A x=b1, contiene le ascisse dei punti del cerchio, ordinati (rispetto ai punti corrispondenti del quadrato) da sinistra a destra e dal basso verso l'alto, ovvero x(1) contiene l'ascissa del punto $Q_{1,1}$ , x(nx) contiene l'ascissa del punto $Q_{nx,1}$ , x(nx+1) contiene l'ascissa del punto $Q_{2,1}$ , ecc., x(n) contiene l'ascissa del punto $Q_{nx,ny}$ . Analogo discorso vale per il vettore y contenente le ordinate dei punti $Q_{ij}$ . Fissati $nx$ e $ny$ , la matrice A ed i termini noti b1 e b2 sono generati mediante la chiamata alla function griglia.m:
[A,b1,b2]=griglia(nx,ny).

Si chiede di:
  1. Analizzare le proprietà della matrice A (dominanza diagonale stretta, simmetria, definita positività) e dire a priori quali metodi iterativi possono essere utilizzati per la risoluzione dei sistemi A x=b1 e A y=b2.
  2. Risolvere i sistemi lineari A x=b1 e A y=b2 con il metodo di Jacobi ( help jacobi), prendendo toll=1.e-6, kmax = 500, x0=rand(n,1).
    Rappresentare su un grafico la storia di convergenza del metodo, ovvero il vettore err, contenente le quantità $\|x_{k+1}-x_k\|_2$ , per $k=0,1,...,$ fino a convergenza o fino a kmax.
    Nel caso in cui in kmax iterazioni non si abbia convergenza, si può portare a convergenza il metodo?
  3. Risolvere i sistemi lineari A x=b1 e A y=b2 con il metodo di Gauss-Seidel ( help gaussseidel), prendendo toll=1.e-6, kmax = 500, x0=rand(n,1).
    Rappresentare su un grafico la storia di convergenza del metodo, ovvero il vettore err, contenente le quantità $\|x_{k+1}-x_k\|_2$ , per $k=0,1,...,$ fino a convergenza o fino a kmax.
    Nel caso in cui in kmax iterazioni non si abbia convergenza, si può portare a convergenza il metodo?
  4. Stimare il raggio spettrale delle matrici di iterazione di Jacobi e Gauss-Seidel.
    I valori ottenuti confermano quanto trovato ai punti percedenti?
  5. Risolvere i sistemi lineari A x=b1 e A y=b2 con il metodo bicgstab ( help bicgstab) con toll=1.e-6, kmax = 500, x0=rand(n,1). Richiamare la function con 5 parametri di output e 7 parametri di input, ponendo il quinto ed il sesto uguali a [].
    Rappresentare su un grafico la storia di convergenza del metodo, ovvero il vettore resvec, contenente le quantità $\|r_k\|_2$ , per $k=0,1/2,1,3/2,...,$ fino a convergenza. Plottare le componenti di resvec con $k$ intero.
  6. Confrontare le storie di convergenza dei tre metodi e commentare i risultati ottenuti.
  7. Rappresentare graficamente la mesh e la funzione $z=f(x,y)=sin(\pi x)sin(\pi y)$ .
    A tale scopo bisogna riordinare i vettori x e y in forma di matrici di ny righe e nx colonne, seguendo l'ordinamento indicato sopra. Scrivere una function (converti.m) che, dato in input un vettore x di dimension n=nx*ny, dia la matrice xh di ny righe e nx colonne.
    Utilizzare il comando mesh per visualizzare la mesh, ed utilizzare il comando meshc per visualizzare la funzione $f(x,y)$ .