/[escript]/branches/doubleplusgood/doc/user/nonlinearPDE.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 <F(u),v> :=\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 > = <F(u^{(\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{<F(u), ve_{i}>}{\|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}).

  ViewVC Help
Powered by ViewVC 1.1.26