docs:math:newton_krylov_hookstep
Differences
This shows you the differences between two versions of the page.
| Both sides previous revision Previous revision Next revision | Previous revision | ||
|
docs:math:newton_krylov_hookstep [2014/12/01 09:07] gibson [The Newton algorithm] |
docs:math:newton_krylov_hookstep [2014/12/01 09:16] (current) gibson [Problem definition] |
||
|---|---|---|---|
| Line 21: | Line 21: | ||
| Let | Let | ||
| - | + | \begin{eqnarray*} | |
| - | <latex> $ \begin{align*} | + | f^t(u) = \text{the time-t map of Navier-Stokes computed by DNS}.\\ |
| - | f^t(u) &= \text{the time-t map of Navier-Stokes computed by DNS}.\\ | + | \sigma = \text{a symmetry of the flow} |
| - | \sigma &= \text{a symmetry of the flow} | + | \end{eqnarray*} |
| - | \end{align*} $ </latex> | + | |
| We seek solutions of the equation | We seek solutions of the equation | ||
| - | <latex> | + | \begin{eqnarray*} |
| G(u,\sigma,t) = \sigma f^t(u) - u = 0 | G(u,\sigma,t) = \sigma f^t(u) - u = 0 | ||
| - | </latex> | + | \end{eqnarray*} |
| //t// can be a fixed or a free variable, as can σ. For plane | //t// can be a fixed or a free variable, as can σ. For plane | ||
| Couette the possible continuously varying variables in σ are | Couette the possible continuously varying variables in σ are | ||
| //a_x,a_z// in the phase shifts //x → x + a_x L_x// and //z → z + a_z L_z//. | //a_x,a_z// in the phase shifts //x → x + a_x L_x// and //z → z + a_z L_z//. | ||
| + | |||
| ===== Discretization ===== | ===== Discretization ===== | ||
| Line 51: | Line 51: | ||
| between continuous fields u and the discrete representation {u_m}. | between continuous fields u and the discrete representation {u_m}. | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| \{u_m\} &= P(u) \\ | \{u_m\} &= P(u) \\ | ||
| u_m &= P_m(u), \;\;\; 1 \leq m \leq M | u_m &= P_m(u), \;\;\; 1 \leq m \leq M | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| Then //x = (u, σ, t)// is an //M+3// dimensional vector | Then //x = (u, σ, t)// is an //M+3// dimensional vector | ||
| - | <latex> | + | \begin{eqnarray*} |
| x_m = (u_1, u_2, ..., u_M, ax, az, t) | x_m = (u_1, u_2, ..., u_M, ax, az, t) | ||
| - | </latex> | + | \end{eqnarray*} |
| and the discrete equation to be solved is | and the discrete equation to be solved is | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| G(x) &= 0 \\ | G(x) &= 0 \\ | ||
| G_m(x) &= P_m(\sigma f^t(u) - u), \;\;\; 1 \leq m \leq M | G_m(x) &= P_m(\sigma f^t(u) - u), \;\;\; 1 \leq m \leq M | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| This is //M// equations in //M+3// unknowns. The indeterminacy follows from the | This is //M// equations in //M+3// unknowns. The indeterminacy follows from the | ||
| Line 81: | Line 80: | ||
| Suppose we start with an initial guess //x// near //x*//. Let | Suppose we start with an initial guess //x// near //x*//. Let | ||
| - | \begin{eqnarray} \begin{align*} | + | \begin{eqnarray*} |
| x^* &= x + dx_N \\ | x^* &= x + dx_N \\ | ||
| u^* &= u + du_N \\ | u^* &= u + du_N \\ | ||
| \sigma^* &= \sigma + d\sigma_N \\ | \sigma^* &= \sigma + d\sigma_N \\ | ||
| t^* &= t + dt_N | t^* &= t + dt_N | ||
| - | \end{align*} \end{eqnarray} | + | \end{eqnarray*} |
| In these equations the //N// subscript signifies the Newton step. Then | In these equations the //N// subscript signifies the Newton step. Then | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| G(x^*) &= G(x + dx_N) \\ | G(x^*) &= G(x + dx_N) \\ | ||
| 0 &= G(x) + DG(x) dx_N + O(|dx_N|^2) | 0 &= G(x) + DG(x) dx_N + O(|dx_N|^2) | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| Drop higher order terms and solve the Newton equation | Drop higher order terms and solve the Newton equation | ||
| - | <latex> | + | \begin{eqnarray*} |
| DG(x) \, dx_N = -G(x) | DG(x) \, dx_N = -G(x) | ||
| - | </latex> | + | \end{eqnarray*} |
| Because //DG// is a //M x (M+3)// matrix, we have to supplement this | Because //DG// is a //M x (M+3)// matrix, we have to supplement this | ||
| Line 105: | Line 104: | ||
| //dx_N//. | //dx_N//. | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| (du_N, du/dt) &= 0 \\ | (du_N, du/dt) &= 0 \\ | ||
| (du_N, du/dx) &= 0 \\ | (du_N, du/dx) &= 0 \\ | ||
| (du_N, du/dz) &= 0 | (du_N, du/dz) &= 0 | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| //( , )// here signifies an inner product, so these are orthogonality | //( , )// here signifies an inner product, so these are orthogonality | ||
| Line 116: | Line 115: | ||
| dimensional system | dimensional system | ||
| - | <latex> | + | \begin{eqnarray*} |
| A \, dx_N = b | A \, dx_N = b | ||
| - | </latex> | + | \end{eqnarray*} |
| ===== Krylov subspace solution of the Newton equation ===== | ===== Krylov subspace solution of the Newton equation ===== | ||
| Line 128: | Line 127: | ||
| with finite differencing and time integration: | with finite differencing and time integration: | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| DG(x) dx &= 1/e \; (G(x + \epsilon dx) - G(x)) \text { where } \epsilon |dx| << 1. \\ | DG(x) dx &= 1/e \; (G(x + \epsilon dx) - G(x)) \text { where } \epsilon |dx| << 1. \\ | ||
| &= 1/e \; P[(\sigma+d\sigma) f^{t+dt}(u+du) - (u+du) - (\sigma f^{t}(u) - u)] \\ | &= 1/e \; P[(\sigma+d\sigma) f^{t+dt}(u+du) - (u+du) - (\sigma f^{t}(u) - u)] \\ | ||
| &= 1/e \; P[(\sigma+d\sigma) f^{t+dt}(u+du) - \sigma f^{t}(u) - du] | &= 1/e \; P[(\sigma+d\sigma) f^{t+dt}(u+du) - \sigma f^{t}(u) - du] | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| The last three entries of //A dx// are simple evaluations of the inner products //(du, du/dt), | The last three entries of //A dx// are simple evaluations of the inner products //(du, du/dt), | ||
| Line 142: | Line 141: | ||
| algorithm to a constrained minimization algorithm called the | algorithm to a constrained minimization algorithm called the | ||
| "hookstep", specially adapted to work with GMRES. | "hookstep", specially adapted to work with GMRES. | ||
| + | |||
| ===== The hookstep algorithm ===== | ===== The hookstep algorithm ===== | ||
| Line 155: | Line 155: | ||
| matrices //Q_n, Q_{n-1},// and //H_n// such that | matrices //Q_n, Q_{n-1},// and //H_n// such that | ||
| - | <latex> | + | \begin{eqnarray*} |
| A Q_n = Q_{n-1} H_n | A Q_n = Q_{n-1} H_n | ||
| - | </latex> | + | \end{eqnarray*} |
| where //Q_n// is //M x n//, //Q_{n-1}// is //M x (n+1)//, and //H_n// is //(n+1) x n//. | where //Q_n// is //M x n//, //Q_{n-1}// is //M x (n+1)//, and //H_n// is //(n+1) x n//. | ||
| Line 165: | Line 165: | ||
| constrained optimization. Namely, we minimize | constrained optimization. Namely, we minimize | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| || Q_{n-1}^T (A \, dx_H - b) || &= || Q_{n-1}^T A Q_n s - Q_{n-1}^T b || \\ | || Q_{n-1}^T (A \, dx_H - b) || &= || Q_{n-1}^T A Q_n s - Q_{n-1}^T b || \\ | ||
| &= || Q_{n-1}^T Q_{n-1} H_n s - Q_{n-1}^T b || \\ | &= || Q_{n-1}^T Q_{n-1} H_n s - Q_{n-1}^T b || \\ | ||
| &= || H_n s - Q_{n-1}^T b || | &= || H_n s - Q_{n-1}^T b || | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| subject to //||dx_H||^2 = || s ||^2 ≤ δ^2//. Note that the quantity in | subject to //||dx_H||^2 = || s ||^2 ≤ δ^2//. Note that the quantity in | ||
| Line 178: | Line 178: | ||
| Then the RHS of the previous equation becomes | Then the RHS of the previous equation becomes | ||
| - | <latex> $ \begin{align*} | + | \begin{eqnarray*} |
| || H_n s - Q_{n-1}^T b || &= || U D V^T s - Q_{n-1}^T b || \\ | || H_n s - Q_{n-1}^T b || &= || U D V^T s - Q_{n-1}^T b || \\ | ||
| &= || D V^T s - U^T Q_{n-1}^T b || \\ | &= || D V^T s - U^T Q_{n-1}^T b || \\ | ||
| &= || D \hat{s} - \hat{b} || | &= || D \hat{s} - \hat{b} || | ||
| - | \end{align*} $ </latex> | + | \end{eqnarray*} |
| where // \hat{s} = V^T s// and //\hat{b} = U^T Q_{n-1}^T b//. Now we need to | where // \hat{s} = V^T s// and //\hat{b} = U^T Q_{n-1}^T b//. Now we need to | ||
| Line 192: | Line 192: | ||
| The solution is | The solution is | ||
| - | <latex> | + | \begin{eqnarray*} |
| \hat{s}_i = (\hat{b}_i d_i)/(d^2_i + \mu) | \hat{s}_i = (\hat{b}_i d_i)/(d^2_i + \mu) | ||
| - | </latex> | + | \end{eqnarray*} |
| (no Einstein summation) for the //μ// such that //||\hat{s}||^2 = δ^2//. | (no Einstein summation) for the //μ// such that //||\hat{s}||^2 = δ^2//. | ||
| This value of //μ// can be found with a 1d Newton search for the zero of | This value of //μ// can be found with a 1d Newton search for the zero of | ||
| - | <latex> | + | \begin{eqnarray*} |
| \Phi(\mu) = ||\hat{s}(\mu)||^2 - \delta^2 | \Phi(\mu) = ||\hat{s}(\mu)||^2 - \delta^2 | ||
| - | </latex> | + | \end{eqnarray*} |
| A straight Newton search a la //μ → μ - Φ(μ)/Φ'(μ)// is suboptimal | A straight Newton search a la //μ → μ - Φ(μ)/Φ'(μ)// is suboptimal | ||
| Line 209: | Line 209: | ||
| solution //\hat{s}(μ)//, we compute the hookstep solution //s// from | solution //\hat{s}(μ)//, we compute the hookstep solution //s// from | ||
| - | <latex> | + | \begin{eqnarray*} |
| dx_H = Q_n s = Q_n V \hat{s} | dx_H = Q_n s = Q_n V \hat{s} | ||
| - | </latex> | + | \end{eqnarray*} |
| ===== Determining the trust-region radius δ ===== | ===== Determining the trust-region radius δ ===== | ||
docs/math/newton_krylov_hookstep.1417453657.txt.gz · Last modified: 2014/12/01 09:07 by gibson
Page Tools
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 3.0 Unported
