Example 11.2.3.

This is now illustrated by the following numerical example. The Olmstead model [343] represents the flow of a layer of viscoelastic fluid heated from below. The equations are

$$\begin{eqnarray*} \displaystyle \frac{\partial u}{\partial t} & = & (1-C)\frac{\... ...artial X^2}+Ru-u^3,\\ B\frac{\partial v}{\partial t} & = & u-v, \end{eqnarray*}$$

where $u$ represents the speed of the fluid and $v$ is related to viscoelastic forces. The boundary conditions are $u(0)=u(1)=0$ and $v(0)=v(1)=0$. After discretization with central differences with grid size $h=1/(n/2)$, the equations may be written as $\dot{x}=f(x)$ with $x=[u_1,v_1,u_2,v_2,\ldots,u_{n/2},v_{n/2}]^*$. We consider the Jacobian matrix $A=\partial f/\partial x$ for $n=1000$ with the parameter values $B=2$, $C=0.1$, and $R=4.7$, evaluated in the trivial steady state solution.

Table 11.2: Residual norms of the inexact rational Krylov method for the Olmstead problem

	Case 1		Case 2
$j$	$\Vert r_j\Vert$	$\Vert f_j\Vert$	$\Vert r_j\Vert$	$\Vert f_j\Vert$
1	$1.3\xx 10^0$	$1.3\xx 10^0$	$1.3\xx 10^0$	$1.3\xx,10^0$
2	$1.8\xx 10^{-1}$	$1.8\xx 10^{-1}$	$1.8\xx 10^{-1}$	$1.8\xx 10^{-1}$
3	$8.6\xx 10^{-3}$	$8.5\xx 10^{-3}$	$8.6\xx 10^{-3}$	$8.6\xx 10^{-3}$
4	$9.8\xx 10^{-4}$	$9.9\xx 10^{-4}$	$5.0\xx 10^{-6}$	$1.7\xx 10^{-6}$
5	$5.0\xx 10^{-5}$	$5.0\xx 10^{-5}$	$3.7\xx 10^{-9}$	$1.0\xx 10^{-13}$
6	$1.1\xx 10^{-5}$	$1.1\xx 10^{-5}$	$1.8\xx 10^{-11}$	$2.1\xx 10^{-20}$
7	$9.3\xx 10^{-7}$	$9.4\xx 10^{-7}$
8	$6.2\xx 10^{-8}$	$6.3\xx 10^{-8}$
9	$2.0\xx 10^{-9}$	$1.9\xx 10^{-9}$
10	$1.2\xx 10^{-10}$	$1.3\xx 10^{-10}$

Figure: Logarithmic plots of residual norms of the inexact rational Krylov method for the Olmstead problem. The circles denote $\Vert f_j\Vert$ and the bullets $\Vert r_j\Vert$.

We solved linear systems by GMRES with 30 iteration vectors. The method was restarted by Morgan's implicitly restarted GMRES [334] keeping the 15 smallest Ritz pairs in the basis until the residual norm satisfied $\Vert s_j\Vert \leq \tau \Vert r_{j-1}\Vert$ with $\tau=1.\,10^{-3}$. We used Algorithm 11.4 for computing the eigenvalue nearest $5$. We consider two cases. For Case 1, we used a fixed $\mu=5$. For Case 2, we used two steps with $\mu=5$ and the remaining steps solved the correction equation (11.5). From the results in Table 11.2 and Figure 11.1, we can see linear convergence for Case 1. Since the linear systems are solved more accurately ($\tau=10^{-3}$) than the speed of the eigenvalue solver, $f_j$ and $r_j$ converge to zero with the same speed. For Case 2, the convergence is linear for Steps 1 and 2, since $\mu_j$ is constant. From the third iteration on, $f_j$ converges quadratically to zero. The $\Vert r_j\Vert$ converge linearly in steps 1 and 2, then we have quadratic convergence on steps 3 and 4, and then linear convergence again with convergence ratio $\tau=10^{-3}$.

The residual norm $r_j$ has two terms. The first term decreases very rapidly when $\mu$ is close to an eigenvalue. The decrease of the second term (the residual of the linear system solver) is often much more difficult when the pole lies close to the spectrum. In practice, we may need to balance the number of outer iterations (the eigenvalue solver) and the number of inner iterations (the linear system solver) by selecting the most optimal $\mu$. This comment is also related to the end of §11.2.3.