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:10] gibson [The Newton algorithm] |
docs:math:newton_krylov_hookstep [2014/12/01 09:16] 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 119: | Line 118: | ||
A \, dx_N = b | A \, dx_N = b | ||
\end{eqnarray*} | \end{eqnarray*} | ||
+ | |||
===== Krylov subspace solution of the Newton equation ===== | ===== Krylov subspace solution of the Newton equation ===== | ||
Line 127: | 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 141: | 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 154: | 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 164: | 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 177: | 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 191: | 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 208: | 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 δ ===== |