# Contents of /branches/doubleplusgood/doc/user/nonlinearPDE.tex

Revision 4345 - (show annotations)
Fri Mar 29 07:09:41 2013 UTC (5 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 14933 byte(s)
Spelling fixes

 1 \chapter{Non-Linear Partial Differential Equations} 2 \label{APP: NEWTON} 3 4 The solution $u_i$ is given as a solution of 5 the nonlinear equation 6 \begin{equation} \label{APP NEWTON EQU 40} 7 \int_{\Omega} v_{i,j} \cdot X_{ij} + v_{i} \cdot Y_{i} \; dx 8 + \int_{\partial \Omega} v_{i} \cdot y_{i} \; ds = 0 9 \end{equation} 10 for all smooth $v_i$ with $v_i=0$ where $q_i>0$ and 11 \begin{equation} \label{APP NEWTON EQU 40b} 12 u_i=r_i \mbox{ where } q_i>0 13 \end{equation} 14 where $X_{ij}$ and $Y_i$ are non-linear functions of the solution $u_k$ and its gradient $u_{k,l}$ 15 and $y_i$ is a function of solution $u_k$. For further convenience we will use the 16 notation 17 \begin{equation} \label{APP NEWTON EQU 40c} 18 :=\int_{\Omega} v_{i,j} \cdot X_{ij} + v_{i} \cdot Y_{i} \; dx 19 + \int_{\partial \Omega} v_{i} \cdot y_{i} \; ds 20 \end{equation} 21 for all smooth $v$ on $\Omega$. If one interprets $F(u)$ as defined above as a functional 22 over the set of admissible functions $v$ 23 equation~(\ref{APP NEWTON EQU 40}) can be written in compact formulation 24 \begin{equation} \label{APP NEWTON EQU 40d} 25 F(u)= 0 26 \end{equation} 27 28 \section{Newton-Raphson Scheme} 29 This equation is iteratively solved by the Newton-Raphson method\index{Newton-Raphson method}, see \cite{Kelley2004a}. 30 Starting with the initial guess 31 $u^{(0)}$ the sequence 32 \begin{equation} \label{APP NEWTON EQU 43} 33 u^{(\nu)}= u^{(\nu-1)} - \delta^{(\nu-1)} 34 \end{equation} 35 for $\nu=1,2,\ldots \;$ generates the (general) Newton-Raphson iteration for the 36 solution $u$. The correction $\delta^{(\nu-1)}$ is the solution of the linear problem 37 \begin{equation} \label{APP NEWTON EQU 43b0} 38 < \fracp{F}{u^{(\nu-1)}} \delta^{(\nu-1)} ,v > = 39 \end{equation} 40 for all smooth $v$ on $\Omega$ with $v_i=0$ where $q_i>0$. 41 where 42 \begin{equation} \label{APP NEWTON EQU 43b} 43 < \fracp{F}{u} \cdot \delta ,v > = 44 \int_{\Omega} \left( \fracp{X_{ij}}{u_{k,l}} v_{i,j}\delta_{k,l} + 45 \fracp{X_{ij}}{u_{k}} v_{i,j}\delta_{k} + \fracp{Y_{i}}{u_{k,l}} v_{i}\delta_{k,l} + 46 \fracp{Y_{i}}{u_{k}} v_{i}\delta_{k} \right) \; dx 47 + \int_{\partial \Omega} 48 \fracp{y_{i}}{u_{k}} v_{i}\delta_{k} \; ds 49 \end{equation} 50 It is assumed that the initial guess $u^{(0)}$ fulfills the constraint~(\ref{APP NEWTON EQU 40b}). 51 The $\delta^{(\nu-1)}$ has to fulfill the homogeneous constraint. 52 Notice that the calculation of $\delta^{(\nu-1)}$ requires the solution of a linear PDE 53 as presented in section~\ref{SEC LinearPDE}. 54 55 The Newton iteration should be stopped in the $k$-th step if for 56 all components of the solution the error of the Newton 57 approximation is lower than the given relative tolerance {rtol}: 58 \begin{equation}\label{APP NEWTON EQU 61} 59 \| u_{i} - u_{i}^{(\nu)} \|_{\infty} \le \mbox{rtol} \cdot \|u_{i} \|_{\infty} \; , 60 \end{equation} 61 for all components $i$ 62 where $\|. \|_{\infty}$ denoted the $L^{sup}$ norm. To measure the quality of the solution approximation 63 on the level of the equation we introduce the weak norm 64 \begin{equation}\label{APP NEWTON EQU 62} 65 \| F(u) \|_{i} := \sup_{v , v=0 \mbox{ where } q_{i}>0 } \frac{}{\|v\|_1} 66 \end{equation} 67 where $(e_{i})_{j}=\delta_{ij}$ and $\|v\|_1=\int_{\Omega} |v| \,dx$ is the $L^1$ norm of $v$ 68 \footnote{In practice a discretization method is applied to solve the update $\delta^{(\nu-1)}$. 69 In this case also an approximation of $\| F(u) \|_{i}$ is calculated taking the maximum over all 70 base function used to represent the solution $u$.}. 71 72 The stopping criterion (\ref{APP NEWTON EQU 61}) is changed to the level of 73 equation. We use the reasonable heuristic but mathematically 74 incorrect argument that the change on the level of the solution and 75 the change on the level of the equation are proportional: 76 \begin{equation}\label{APP NEWTON EQU 64} 77 \frac{ \| u_i - u^{(\nu)}_i \|_{\infty} }{ \| 0 - F(u^{(\nu)}) \|_{i} } = 78 \frac{ \| \delta^{(\nu)}_i \|_{\infty} }{ \| F(u^{(\nu)}) - F(u^{(\nu-1)}) \|_{i} } \; . 79 \end{equation} 80 where we assume that that component $\nu$ $F(u)$ is mainly controlled by component $\nu$ of the solution. 81 82 83 We assume that the term $F(u^{(\nu)})$ can be neglected versus 84 $F(u^{(\nu-1)})$ since $u^{(\nu)}$ is a better approximation, 85 and use the stopping criterion in the formulation: 86 \begin{equation} \label{APP NEWTON EQU 65} 87 \| F(u^{(\nu)}) \|_i \le 88 \frac{ \| F(u^{(\nu-1)})\|_{i} \cdot \|u_{i} \|_{\infty} } { \| \delta^{(\nu)}_i \| _{\infty} } 89 \,\mbox{\it rtol}\, =:\, \mbox{\it qtol}_i \; , 90 \end{equation} 91 which has to hold for all components $\nu$. 92 Now {\it qtol} defines a tolerance for the level of equation. This stopping criterion is not free of problems, because a 93 decrease of the defect $F(u^{(\nu)})$ coupled with a constant 94 correction $\delta^{(\nu)}$ suggests a good approximation. 95 But the quality of the approximation $u^{(\nu)}$ is ensured if the 96 Newton iteration converges quadratically. This convergence behavior 97 is given by the error estimation: 98 \begin{equation} 99 \| u - u^{(\nu)} \|_{\infty} \le C \; 100 \| u - u^{(\nu-1)} \|_{\infty}^2 101 \end{equation} 102 with a positive value $C$, see \cite{Kelley2004a}. Therefore a quadratic 103 convergence of the Newton iteration can be assumed if the corrections 104 of the current and the last step fulfill the following condition: 105 \begin{equation} \label{APP NEWTON EQU 66} 106 \max_{i} 107 \frac{\| \delta^{(\nu)}_i \|_{\infty} }{\| \delta^{(\nu-1)}_i \|_{\infty} } 108 < \frac{1}{5} \;, 109 \end{equation} 110 where the limit $\frac{1}{5}$ was found by a large number of experiments. 111 The approximation $u^{(\nu)}$ is accepted if the conditions 112 (\ref{APP NEWTON EQU 65}) and (\ref{APP NEWTON EQU 66}) hold. Consequently a 113 safe approximation requires at least two Newton steps. 114 115 To stop a divergent iteration, which occurs for a bad initial solution, 116 the norms of the defects for the $(k-1)$-th and $k$-th 117 Newton step are compared. Here we use the estimation 118 \begin{equation} \label{APP NEWTON EQU 67} 119 \| F(u^{(\nu)}) \|_i \le 120 \gamma \| F(u^{(\nu-1)}) \|_i 121 \end{equation} 122 for the defects. The value $0<\gamma<1$ depends on $F$ 123 and the distance of the initial guess $u^{(0)}$ to the true solution $u$. 124 Since the constant $\gamma$ is unknown and can be close to one, 125 convergence is assumed if the following (weaker) condition holds 126 for all components $i$: 127 \begin{equation} \label{APP NEWTON EQU 68} 128 \| F(u^{(\nu)}) \|_i 129 < \| F(u^{(\nu-1)}) \|_i \; . 130 \end{equation} 131 If condition (\ref{APP NEWTON EQU 68}) fails, divergence may begin. Therefore 132 under-relaxation is started. Beginning with $\omega=1$ 133 the Newton-Raphson iteration is computed by 134 \begin{equation} \label{APP NEWTON EQU 69} 135 u^{(\nu)} = u^{(\nu-1)} - \omega \, \delta^{(\nu)} 136 \end{equation} 137 instead of (\ref{APP NEWTON EQU 43}). If this new iteration fulfills the 138 condition (\ref{APP NEWTON EQU 68}) of decreasing defect, we accept it. Otherwise 139 we put $\omega \rightarrow \frac{\omega}{2}$, recompute $u^{(\nu)}$ from equation 140 (\ref{APP NEWTON EQU 69}) and retry condition (\ref{APP NEWTON EQU 68}) for the 141 new $u^{(\nu)}$, and so on until either condition (\ref{APP NEWTON EQU 68}) 142 holds or $\omega$ becomes to small ($\omega < \omega_{lim}=0.01$). In the latter case the 143 iteration gives up. The under-relaxation 144 converges only linearly for $\omega<1$. it is a rather robust 145 procedure. 146 which switches back to $\omega=1$ as soon as possible. 147 The price for the robustness is the additional computation of the 148 defects. 149 150 Due to the quadratic convergence near the solution the error decreases 151 rapidly. Then the solution will not change much and it will not be 152 necessary to mount a new coefficient matrix in each iteration step. 153 This algorithm is called the simplified Newton method. It converges 154 linearly by 155 \begin{equation} 156 \| u_i - u^{(\nu)}_i \|_{\infty} \le \gamma^{\nu} 157 \| u_i - u ^{(0)}_i \|_{\infty} \; , 158 \end{equation} 159 where $\gamma$ equals the $\gamma$ in the estimation (\ref{APP NEWTON EQU 67}). 160 If the general iteration converges quadratically 161 (condition (\ref{APP NEWTON EQU 66}) holds) and we have 162 \begin{equation} 163 \| F(u^{(\nu)}) \|_i < 0.1 164 \| F(u^{(\nu-1)}) \|_i 165 \end{equation} 166 for all components $\nu$, 167 we can expect $\gamma \le 0.1$. Then the simplified iteration produces 168 one digit in every step and so we change from the general to the 169 simplified method. The slow' convergence requires more iteration steps, 170 but the expensive mounting of the coefficient matrix is saved. 171 172 The lienar PDE~(\ref{APP NEWTON EQU 43b}) is solved with with a certain tolerance, namely when 173 the defect of the current 174 approximation of $\delta^{(\nu)}$ relative to the defect of 175 the current Newton iteration is lower than {LINTOL}. To avoid 176 wasting CPU time the FEMLIN iteration must be controlled by an 177 efficient stopping criterion. We set 178 \begin{equation} 179 \mbox{LINTOL} = 0.1 \cdot \max ( \left( \frac{\|\delta^{(\nu)}_i\|_{\infty}}{\|u_i^{(\nu)}\|_{\infty}} \right)^2 180 ,\min_{i} \frac{\mbox{qtol}_i}{\|F(u^{(\nu)})\|_{i}} ) 181 \end{equation} 182 but restrict {LINTOL} by 183 \begin{equation} 184 10^{-4} \le \mbox{LINTOL} \le 0.1 \quad \mbox{ and } \quad 185 \mbox{LINTOL}=0.1 \; \mbox{ for }\nu=0 \;. 186 \end{equation} 187 The first term means that it would be useless to compute digits by the 188 linear solver, which are overwritten by the next Newton step. In the 189 region of quadratic convergence the number of significant digits is 190 doubled in each Newton step, i.e.~later digits are overwritten by the 191 following Newton-Raphson correction, see \cite{Schoenauer1981a}. The second 192 mean that no digits should be computed which are in significance 193 below the prescribed tolerance {\it rtol}. The number $0.1$ is a safety factor' to take care 194 of the coarse norm estimations. Figure~\ref{APP NEWTON PIC 61} shows the workflow of 195 the Newton-Raphson update algorithm. 196 197 \begin{figure} 198 \begin{center} 199 { 200 \unitlength0.92mm 201 \begin{picture}(200,235) \thicklines 202 203 \put(-10,-065){\begin{picture}(170,280) \thicklines 204 205 \newsavebox{\BZar} 206 \savebox{\BZar}(0,0) { 207 \thicklines 208 \put(0,10){\line(-3,-1){30}} 209 \put(0,10){\line(3,-1){30}} 210 \put(0,-10){\line(-3,1){30}} 211 \put(0,-10){\line(3,1){30}} 212 \put(0,20){\vector(0,-1){10}} 213 \put(0,-10){\vector(0,-1){10}} 214 \put(30,0){\vector(1,0){25}} } 215 \newsavebox{\BZal} 216 \savebox{\BZal}(0,0) { 217 \thicklines 218 \put(0,10){\line(-3,-1){30}} 219 \put(0,10){\line(3,-1){30}} 220 \put(0,-10){\line(-3,1){30}} 221 \put(0,-10){\line(3,1){30}} 222 \put(0,20){\vector(0,-1){10}} 223 \put(0,-10){\vector(0,-1){10}} 224 \put(-30,0){\vector(-1,0){10}} } 225 226 \put(20,285){\framebox(60,20){\parbox{48mm} 227 {Start: \\ $\nu=0$ , $\omega=1$ \\ 228 calculate $F(u^{(0)})$ }} } 229 \put(50,285){\vector(0,-1){10}} 230 \put(20,265){\framebox(60,10){\parbox{48mm} 231 {next iteration: $\nu \leftarrow \nu+1$}} } 232 \put(50,265){\vector(0,-1){10}} 233 \put(20,240){\framebox(60,15){\parbox{54mm} 234 {\hspace*{2.00mm} Solve \\ 235 $\fracp{F}{ u^{(\nu-1)}} \delta^{(\nu)} = 236 F(u^{(\nu-1)})$}} } 237 \put(50,240){\vector(0,-1){10}} 238 \put(20,220){\framebox(60,10){\parbox{48mm} 239 {$\omega = \min(2\,\omega,1)$}} } 240 \put(50,220){\vector(0,-1){10}} 241 \put(120,225){\line(-3,-1){30}} 242 \put(120,225){\line(3,-1){30}} 243 \put(120,205){\line(-3,1){30}} 244 \put(120,205){\line(3,1){30}} 245 \put(110,214){$\omega \le \omega_{lim}$ ?} 246 \put(83,219){no} 247 \put(153,219){yes} 248 \put(090,215){\vector(-1,0){40}} 249 \put(150,215){\line(1,0){05}} 250 \put(20,200){\framebox(60,10){\parbox{48mm} 251 {$u^{(\nu)} = u^{(k-1)} - \omega\delta^{(\nu)}$}} } 252 \put(50,180){\usebox{\BZar}} 253 \put(30,179){$F(u^{(\nu)}) < F(u^{(\nu-1)})$ ?} 254 \put(83,184){no} 255 \put(55,165){yes} 256 \put(105,175){\framebox(30,10){\parbox{23mm} 257 {$\omega=\omega/2$}} } 258 \put(120,185){\vector(0,1){20}} 259 260 \put(20,145){\framebox(60,15){\parbox{48mm} 261 {evaluate stopping criteria~(\ref{APP NEWTON EQU 65}) \\ 262 and if not simplified mode~(\ref{APP NEWTON EQU 66}) } } } 263 \put(50,145){\vector(0,-1){10}} 264 265 \put(50,125){\usebox{\BZar}} 266 \put(34,120){\shortstack{stopping criterion \\ satisfied ?}} 267 \put(83,129){yes} 268 \put(55,110){no} 269 \put(105,117.5){\framebox(30,15){\parbox{23mm} 270 {Newton ends \\ successfully!}} } 271 \put(135,125){\vector(1,0){20}} 272 \put(50,095){\usebox{\BZal}} 273 \put(26,093.5){$F(u^{(\nu)}) < 0.1 F(u^{(\nu-1)})$ ?} 274 \put(12,099){no} 275 \put(55,080){yes} 276 \put(20,065){\framebox(60,10){\parbox{48mm} 277 {switch to simplified Newton!}} } 278 \put(20,070){\vector(-1,0){10}} 279 280 \put(155,215){\vector(0,-1){60}} 281 \put(140,145){\framebox(30,10){\parbox{23mm} 282 {Newton fails!}} } 283 \put(155,145){\vector(0,-1){070}} 284 \put(155,070){\oval(40,10)\makebox(0,0){END}} 285 286 \put(010,280){\line(0,-1){210}} 287 \put(010,280){\vector(1,0){40}} 288 289 290 \end{picture} } 291 \end{picture} } 292 \end{center} 293 294 \caption{\label{APP NEWTON PIC 61}Flow diagram of the Newton-Raphson algorithm} 295 \end{figure} 296 297 \section{Local Sensitivity Analysis} 298 If the coefficients 299 $X_{ij}$, $Y_i$ and $y_{i}$ in equation~(\ref{APP NEWTON EQU 40}) 300 depend on a vector of input factors $f_i$ index{input factor} and its gradient one is interested in how the solution $u_i$ 301 is changing if the input factors are changed. This problem is called a local sensitivity 302 analysis \index{local sensitivity analysis}. If $u(f)$ denotes the 303 solution of equation~(\ref{APP NEWTON EQU 40}) for the input factor $p$ 304 and $u(f+\alpha \cdot g)$ denotes the solution 305 for a perturbed value $f+\alpha \cdot g$ for input factor $f$ where 306 $q$ denotes the direction of perturbation and $\alpha$ is the small scaling factor. 307 The derivative of the solution in the direction $g$ is defined as 308 \begin{equation} 309 \fracp{u}{g} : = \lim_{\alpha \rightarrow 0} \frac{ u(f+\alpha \cdot g) - u(f)}{\alpha} 310 \end{equation} 311 In practice one needs to distinguish between the cases of a spatially constant 312 spatially variable. In the first case $g$ is set to a unit vector while 313 in a second case an appropriate function needs to be given for $g$. 314 315 The function $\fracp{u}{g}$ is calculated from solving the equation 316 \begin{align} \label{APP NEWTON EQU 100} 317 \int_{\Omega} \left( \fracp{X_{ij}}{u_{k,l}} v_{i,j} \left(\fracp{u_k}{g}\right)_{,l} + 318 \fracp{X_{ij}}{u_{k}} v_{i,j}\fracp{u_k}{g} + \fracp{Y_{i}}{u_{k,l}} v_{i}\left(\fracp{u_k}{g}\right)_{,l} + 319 \fracp{Y_{i}}{u_{k}} v_{i}\fracp{u_k}{g}\right) \; dx 320 + \int_{\partial \Omega} 321 \fracp{y_{i}}{u_{k}} v_{i}\fracp{u_k}{g}\; ds \\ 322 + \int_{\Omega} v_{i,j} \left( \fracp{X_{ij}}{f_{k,l}} g_{k,l} + \fracp{X_{ij}}{f_{k}} g_k \right) 323 + v_{i} \left( \fracp{Y_{i}}{f_{k,l}} g_{k,l} + \fracp{Y_{i}}{f_{k}} g_k \right) \; dx 324 + \int_{\partial \Omega} v_{i} 325 \fracp{y_{i}}{f_{k}} g_k\; ds 326 \end{align} 327 for all smooth $v$ on $\Omega$ with $v_i=0$ where $q_i>0$ 328 for the unknown sensitivity $\fracp{u}{g}$. 329 Notice that this equation is similar to the equation which needs to be solved for 330 the Newton-Raphson correction $\delta_k$, see equation~(\ref{APP NEWTON EQU 43b}).