\section{Algoritmo Newton-Raphson per sistemi di equazioni nonlineari} \subsection{Iterazione base} Consideriamo un sistema lineare di $n$ equazioni \begin{equation} \label{eqNR1} \ve{R}\left(\ve{u}\right)=\ve{F}\left(\ve{u}\right) \end{equation} nelle $n$ componenti incognite del vettore $\ve{u}$, con \[ \ve{R}:\ve{u} \rightarrow \mathbb{R}^n \;,\; \ve{u} \in C \subseteq \mathbb{R}^n \] \[ \ve{F}:\ve{u} \rightarrow \mathbb{R}^n \;,\; \ve{u} \in C \subseteq \mathbb{R}^n \] funzioni vettoriali di variabile vettoriale. Nel caso specifico della soluzione di sistemi di equazioni derivate dagli equilibri nodali di strutture discretizzate con metodo FEM, ho \begin{description} \item[$\ve{u}$:] vettore contenente le componenti di spostamento/rotazione nodale dalla configurazione indeformata (incognite); \item[$\ve{F}\left(\ve{u}\right)$:] vettore contenente le componenti di forza/coppia nodale applicate dall'esterno sul sistema, supposte note per una data configurazione della struttura \footnote{la dipendenza delle forze esterne dalla configurazione è stata inserita in questa trattazione per amor di generalità, nonché per includere entro la stessa fenomeni nonlineari specifici quali l'instabilità dei rotori o la pressurizzazione di membrane}; \item[$\ve{R}\left(\ve{u}\right)$:] vettore contenente le componenti di azione nodale necessarie a mantenere la struttura in equilibrio nello stato deformativo associato al vettore spostamenti nodali $\ve{u}$, \emph{ovvero} vettore contenente le componenti di azione nodale associate (uguali e contrarie) alle reazioni elastiche della struttura costretta in stato deformato. Nel caso particolare di sistema elastico lineare $\ve{R}\left(\ve{u}\right) = \ma{K}\ve{u}$, con $K$ matrice di rigidezza. \end{description} Si nota che tale interpretazione dei termini dell'equazione \ref{eqNR1} è appropriata nel caso le condizioni al contorno siano alle sole forze. Nel caso in cui siano definiti vincoli di spostamento nodale imposto alcune coppie di termini coniugati $R_l\left(\ve{u}\right),F_l\left(\ve{u}\right)$ risulteranno modificate in quanto all'equazione di equilibrio $l$-esima si sostituisce l'identità cinematica tra spostamenti incognito ed imposto. Una scrittura alternativa prevede la definizione e l'annullamento di un termine di residuo \begin{equation} \label{eqNR1r} \ve{r}\left(\ve{u}\right)=\ve{0} \end{equation} con \begin{equation} \label{eqdefres} \ve{r}\left(\ve{u}\right) = \ve{R}\left(\ve{u}\right)-\ve{F}\left(\ve{u}\right) \end{equation} Tale scrittura permette di riassumere in un unico termine le variazioni in $\ve{u}$ di forze e reazioni elastiche, per cui risulta vantaggioso procedere con tale notazione. Il metodo N-R è costruito a partire dallo sviluppo in serie di Taylor al primo ordine dell'equazione \ref{eqNR1r} nell'intorno di un punto di iterato $i$-esimo $\ve{u}^i$, ossia \begin{equation} \label{taylorexpansion} \ve{r} \left( \ve{u}^\ast \right) = \ve{r} \left( \ve{u}^i \right) + \ma{J}_r \left( \ve{u^i} \right) \cdot \left( \ve{u}^\ast-\ve{u}^i \right) + o \left( \ve{u}^\ast-\ve{u}^i \right) =\ve{0}. \end{equation} In tale espressione\footnote{ qui utilizzo la notazione per cui $\left[ A \right]_{i,j}$ è l'elemento alla $i$-esima riga e $j$-esima colonna di $A$} \begin{equation} \left[ \ma{J}_r \left( \ve{u}^i \right) \right]_{l,m} = \left[ \ma{J}_r^i \right]_{l,m} = \left. \frac{\partial r_l}{\partial u_m} \right|_{\ve{u}=\ve{u}^i} , \quad l,m=1 \ldots n \end{equation} è lo Jacobiano della funzione residuo $r$ calcolato al punto $\ve{u}^i$, mentre $\ve{u}^\ast$ è la soluzione esatta incognita; data la natura della funzione residuo risulta inoltre \begin{equation} \ma{J}_r^i = \ma{J}_R^i - \ma{J}_F^i \end{equation} ove \begin{equation} \left[ \ma{J}_R^i \right]_{l,m} = \left. \frac{\partial R_l}{\partial u_m} \right|_{\ve{u}=\ve{u}^i} ,\quad \left[ \ma{J}_F^i \right]_{l,m} = \left. \frac{\partial F_l}{\partial u_m} \right|_{\ve{u}=\ve{u}^i} , \quad l,m=1 \ldots n \end{equation} Notiamo che nel caso $\ve{F}$ sia costante in $\ve{u}$ si ha \begin{equation} \ma{J}_F^i=\ma{0} \Rightarrow \ma{J}_r^i=\ma{J}_R^i \end{equation} Trascurando il termine di ordine superiore, l'identità \ref{taylorexpansion} non risulterà più strettamente verificata; sostituendo tuttavia alla soluzione esatta $\ve{u}^\ast$ un meno pretenzioso termine di iterato successivo $\ve{u}^{i+1}$ ottengo la forma \begin{equation} \ma{J}_r^i \left( \ve{u}^{i+1}-\ve{u}^i \right) = - \ve{r} \left( \ve{u}^i \right) \end{equation} da cui, supponendo $\ma{J}_r^i$ non singolare \footnote{la notazione $\ve{x} = \ma{A} \backslash \ve{b}$ indica la soluzione di un sistema lineare di matrice $\ma{A}$, vettore termini noti $\ve{b}$, e l'assegnazione del risultato al vettore delle incognite $\ve{x}$.} \begin{equation} \label{getnewiterval} \ve{u}^{i+1} = \ve{u}^i - \ma{J}_r^i \backslash \ve{r} \left( \ve{u}^i \right). \end{equation} In notazione equivalente, posto \begin{equation} \begin{array}{ccc} \ve{\delta u}^{i} \equiv \ve{u}^{i+1} -\ve{u}^{i}, & \ma{K}^i \equiv \ma{J}_r^i, & \ve{r}^i \equiv \ve{r} \left( \ve{u}^i \right) \end{array}, \end{equation} abbiamo l'espressione \begin{equation} \label{getnewiterval2} \ve{\delta u}^{i} = - \ma{K}^i \backslash \ve{r}^i \end{equation} nella quale è evidenziato il ruolo di matrice di rigidezza \emph{tangente} che $\ma{J}_r^i$ ricopre, in analogia al caso lineare. Fornito infine un vettore di primo tentativo $\ve{u}^0$, la \ref{getnewiterval} definisce una successione iterativa di termini che \emph{potenzialmente} converge alla soluzione esatta $\ve{u}^\ast$. Dovendo limitare ad un numero finito le iterazioni di \ref{getnewiterval} o \ref{getnewiterval2}, occorre definire dei criteri di convergenza (accettazione risultato, fine iterazione) nella forma di \begin{itemize} \item Convergenza ai residui: stop iterazione implicato da $ \left\Vert \ve{r} \left( \ve{u}^{i+1} \right) \right\Vert \le \epsilon_r$, con $\epsilon_r$ errore (assoluto) ammesso ai residui; \item Convergenza alle incognite (\emph{``agli spostamenti''} nello specifico FEM): stop iterazione implicato da $ \left\Vert \ve{u}^{i+1} - \ve{u}^{i} \right\Vert = \left\Vert \ve{\delta u}^{i} \right\Vert \le \epsilon_u$, con $\epsilon_u$ errore (assoluto) ammesso alle incognite. \end{itemize} L'implementazione di valori di tolleranza \emph{relativa} per incognite e residui è ovviamente possibile (e in generale auspicabile) una volta individuati opportuni termini di riferimento in valore assoluto per forze e spostamenti. La determinazione di tali valori di tolleranza ammissibile è un aspetto cruciale dell'analisi nonlineare in quanto la richiesta di una eccessiva precisione aumenta inutilmente i tempi di calcolo, mentre tolleranze molto lasche possono produrre soluzioni inaccurate. Una possibile implementazione dell'algoritmo Newton-Raphson, iterazione base è riportata come Algorithm \ref{NRbase}. \begin{algorithm} \caption{Iterato N-R base}\label{NRbase} \begin{algorithmic} \Require $\ve{u}^0$ \Require $\ve{r}\left(\ve{u^i}\right)$ calcolabile $\forall \ve{u}^i$ di iterato \Require $\ma{J_r}\left(\ve{u^i}\right)$ calcolabile $\forall \ve{u}^i$ di iterato, e non singolare \State $i \gets 0$ \While{$ \left( \left\Vert \ve{r} \left( \ve{u}^{i} \right) \right\Vert > \epsilon_F \right) \vee \left( i \geq 1 \wedge \left( \left\Vert \ve{u}^{i} - \ve{u}^{i-1} \right\Vert > \epsilon_u \right) \right) $} \State $\ma{K}^i \gets \ma{J_r}\left(\ve{u}^i \right)$ \State $\ve{r}^i \gets \ve{r} \left( \ve{u}^i \right)$ \State $\ve{\delta u}^{i} \gets - \ma{K}^i \setminus \ve{r}^i $ \State $\ve{u}^{i+1} \gets \ve{u}^{i} + \ve{\delta u}^{i} $ \State $i \gets i+1$ \EndWhile \State $\ve{u}^\ast \gets \ve{u}^{i}$ \end{algorithmic} \end{algorithm} \begin{comment} \subsection{Doppia iterazione, definizione della storia del carico} , ossia conservativo. Conservativi sono ad esempio i sistemi elastici sottoposti a grandi deformazioni, i sistemi in cui il carico è funzione degli spostamenti, i sistemi con contatti monolateri in assenza di attrito. Non conservativi sono invece i sistemi presentanti fenomeni dissipativi quali deformazioni plastiche, o interfacce di contatto con attrito. Sotto l'ipotesi di conservatività uno stato del sistema elastico è univocamente\footnote{In assenza di moti rigidi o \emph{assimilati}, ossia moti a variazione di energia potenziale elastica nulla; in generale la soluzione agli spostamenti è unica a meno di quote arbitrarie in tali moti.} definito una volta fornito un valore (vettore) per il termine noto. Nei sistemi a comportamento non lineare non è in generale sufficiente definire una condizione di sollecitazione (carichi e/o cedimenti imposti) per caratterizzare univocamente lo stato deformativo della struttura. Esistono infatti una molteplicità di fenomeni che possono rendere tale stato funzione in quanto la configurazione finale è funzione del percorso col quale, a partire da uno stato iniziale tipicamente scarico/indeformato, tale valore viene raggiunto. L'evoluzione del termine noto - e quindi l'evoluzione del sistema - può essere convenientemente descritta in forma parametrica, introducendo un parametro tempo $t$ monotonicamente crescente. Avremo quindi, su di un intervallo di tempo $t\in\left[t_a,t_b\right]$ \begin{equation} \label{eqNRt} \ve{R}\left(\ve{u}\left( t \right) \right)=\ve{F}\left( t \right) \end{equation} Non potendo gestire tale relazione in forma continua, si provvede a discretizzare l'intervallo temporale mediante un numero finito di campionamenti $t_j$, con $j=0 \ldots m$, $t_0 =t_a$ e $t_m =t_b$. Si definiscono di conseguenza $m+1$ valori di termine noto $\ve{F}_j = \ve{F}\left( t_j \right)$, ed una sequenza di iterati base Newton-Rhapson nella forma descritta (ad esempio) dallo pseudocodice Algorithm \ref{NRstruct}, che restituiscono $m$ soluzioni $\ve{u}_j^\ast$, $j=1 \ldots m$\footnote{L'apice $\square^\ast$ è utilizzato per indicare l'avvenuta -- o meglio, \emph{dichiarata} -- convergenza dell'iterato N-R base}. Si nota che è necessario fornire una soluzione $\ve{u}_0^\ast$ per $t_0=t_a$ in equilibrio con le forze iniziali $\ve{F}_0=\ve{F}\left(t_a\right)$, che in generale è presa a sistema scarico ed indeformato; la soluzione $\ve{u}_{j-1}^\ast$ giunta a convergenza per il valore di termine noto $\ve{F}_{j-1}$ è inoltre utilizzata come partenza per la successiva iterazione base N-R, a termine noto $\ve{F}_{j}$. \begin{algorithm} \caption{Algoritmo N-R per applicazioni strutturali (parziale..)}\label{NRstruct} \begin{algorithmic} \Require $\ve{u}^\ast_0$ in equilibrio con $\ve{F}_0$, potenzialmente ambedue $=\ve{0}$ \Require $\ve{R}\left(\ve{u}^i_j\right)$ calcolabile $\forall \ve{u}^i_j$ di iterato \Require $\ma{J_R}\left(\ve{u}^i_j\right)$ calcolabile $\forall \ve{u}^i_j$ di iterato, e non singolare \For{$j = 1 \to m$} \State $\ve{u}_j^0 \gets \ve{u}_{j-1}^\ast$ %\Comment{Parto dall'ultima configurazione in convergenza} \State $i \gets 0$ %\Comment{Indice di ciclo iterazione N-R} \State $conv \gets $ False \While{$\left(\left(i\leq imax\right)\; \wedge \; \left(!\;conv\right)\right)$ } \State $\ma{K}^i_j \gets \ma{J_R}\left(\ve{u}^i_j \right)$ \State $\ve{R}^i_j \gets \ve{R} \left( \ve{u}^i_j \right)$ \State $\ve{\Delta u}^{i}_j \gets \ma{K}^i_j \setminus \left( \ve{F}_j - \ve{R}^i_j \right)$ \State $\ve{u}_j^{i+1} \gets \ve{u}_j^{i} + \ve{\Delta u}^{i}_j $ \State $i \gets i+1$ \If{$ \left( \left\Vert \ve{F}_j - \ve{R} \left( \ve{u}_j^{i} \right) \right\Vert \leq \epsilon_F \right) \wedge \left( \left\Vert \ve{u}_j^{i} - \ve{u}_j^{i-1} \right\Vert \leq \epsilon_u \right) $} \State{$conv \gets$ True} %\Else % \State $i \gets i+1$ \EndIf \EndWhile \If{$conv=$True} \State $\ve{u}_j^\ast \gets \ve{u}_j^{i}$ \Else \State{Break.}\Comment{"Houston, we've had a problem."} \EndIf \EndFor \end{algorithmic} \end{algorithm} Poiché l'evoluzione del sistema entro l'intervallo temporale discreto \[t_{\left( j-1 \right) } \rightarrow t_{j}\] è demandata all'iterazione N-R base, ai cui step intermedi non è associato significato \emph{fisico} -- tant'è che non è rispettata la condizione di equilibrio nodale, e poichè tale evoluzione risulterebbe rilevante ai fini della soluzione nel caso siano presenti fenomeni non conservativi, conviene in presenza di questi ultimi utilizzare discretizzazioni temporali molto fini pena accumulo di errore nella soluzione. \end{comment} \subsection{Caso unidimensionale: algoritmo di Newton e soluzione grafica} blablabla, vedi Figura \ref{fig:nr1d}. \begin{figure} \centering \includegraphics[scale=1]{NR-1d.pdf} \caption{Costruzione grafica per l'iterato N-R, caso N=1} \label{fig:nr1d} \end{figure} \subsection{Caso bidimensionale: accenno di soluzione grafica} blablabla, vedi Figura \ref{fig:nr2d}a e Figura \ref{fig:nr2d}b. \begin{figure} \centering \includegraphics[scale=0.8]{NR-2d_v1.pdf} \caption{Costruzione grafica per l'iterato N-R, caso N=2} \label{fig:nr2d} \end{figure}