We’re rewarding the question askers & reputations are being recalculated! Read more.
21

In general, all Krylov methods essentially seek a polynomial that is small when evaluated on the spectrum of the matrix. In particular, the $n$th residual of a Krylov method (with zero initial guess) can be written in the form $$ r_n = P_n (A) b $$ where $P_n$ is some monic polynomial of degree $n$ . If $A$ is diagonalizable, with $A=V\Lambda V^{-1}$, we ...


17

On norms As an addendum to Reid.Atcheson's answer, I would like to clarify some issues regarding norms. At the $n^{\mathrm{th}}$ iteration, GMRES finds the polynomial $P_n$ that minimizes the $2$-norm of the residual $$r_n = A x_n - b = \big(P_n(A) - 1 \big)b - b = P_n(A) b .$$ Suppose $A$ is SPD, so $A$ induces a norm and so does $A^{-1}$. Then $$\begin{...


15

The matrix you're referring to is positive definite. The eigenvalues of the matrix must be real, because symetric matrices are equal to their own conjugate transpose, and are thus Hermitian. All eigenvalues of Hermitian matrices are real. If all entries of the matrix are positive or zero and the matrix is weakly diagonally dominant, then all eigenvalues of ...


15

There are two main classes of solutions to be discussed in this regard. "Sufficiently" Smooth Solutions In Strang's classical paper it is shown that the Lax equivalence theorem (i.e., the idea that consistency plus stability implies convergence) extends to nonlinear PDE solutions if they have a certain number of continuous derivatives. Note that that ...


11

Iterative methods in a nutshell: Stationary methods are in essence fixed point iterations: To solve $Ax=b$, you pick an invertible matrix $C$ and find a fixed point of $$ x = x + Cb- CAx $$ This converges by Banach's fixed point theorem if $\|I-CA\|<1$. The various methods then correspond to a specific choice of $C$ (e.g., for Jacobi iteration, $C=D^{-1}...


11

In addition to Christian's answer, it's also worth noting that for linear convergence you have $e_{k+1} \le \lambda_1 e_k$ where you have $\lambda_1<1$ if the method converges. On the other hand, for quadratic convergence you have $e_{k+1} \le \lambda_2 e_k^2$ and the fact that a method converges does not necessarily imply that $\lambda_2$ must be smaller ...


10

It isn't so much that we want to compare the $p$-refinement and $h$-refinement errors directly, instead we want to compare the convergence properties (e.g. speed) of each refinement strategy. This requires more knowledge of the constant in the apriori error estimate. We'll illustrate by looking at the apriori error estimate of a discontinuous Galerkin ...


10

It is because we typically neglect higher order terms in error estimates. For example, we can show that $$ \|e\| \le C(u) h^2 + {\cal O}(h^3). $$ The point is that when $h$ is small, the cubic term is small and can be neglected. In fact, when $h$ is small, you can observe quadratic convergence. But whenever $h$ is not small (where "small" is relative to ...


10

Yes, this is the standard Aubin-Nitsche (or duality) trick. The idea is to use the fact that $L^2$ is its own dual space to write the $L^2$-norm as an operator norm $$\|u\|_{L^2} = \sup_{\phi\in L^2\setminus\{0\}} \frac{(u,\phi)}{\|\phi\|_{L^2}}.$$ We thus have to estimate $(u-u_h,\phi)$ for arbitrary $\phi\in L^2$. To do that, we "lift" $u-u_h$ to $H^1_0$ ...


10

Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.", Birkhäuser 2018, gives the following example due to Muller, specifically constructed to deliver incorrect results with floating-point evaluation: $$ {u_{0} = 2,\\ u_{1} = -4,\\ u_{n} = 111 - \frac{1130}{u_{n-1}} + \frac{3000}{u_{n-1}u_{n-2}},\>\>\>\>n \ge 2.} $$ ...


9

In practice, yes. While $e_k$ is still large, the rate coefficient $\lambda$ will dominate the error rather than the q-rate. (Note that these are asymptotic rates, so the statements you linked to only hold for the limit as $k\to\infty$.) For example, for first order methods in optimization you often observe an initially fast decrease in error, which then ...


9

For a single rational equation in the complex domain, the basin of attraction is fractal, the compelement of a so-called Julia set. http://en.wikipedia.org/wiki/Julia_set . For theory with some nice online figures, see, e.g., http://mathlab.mathlab.sunysb.edu/~scott/Papers/Newton/Published.pdf http://hera.ugr.es/doi/15019160.pdf Even the ''globalized'' ...


9

Let me first answer all the questions: What is the theoretical convergence rate for an FFT Poison solver? The theoretical convergence is exponential as long as the solution is sufficiently smooth. How fast should this energy converge? The Hartree energy $E_H$ should converge exponentially for a sufficiently smooth solution. If the solution is less ...


9

$\|x^{(k)}-x^*\|$ is the error in the $k$th term, call it $E_k$. For a "good" numerical method, we want the approximation to get closer and closer to the desired result so $E_k$ has to decrease to zero. If the error is guaranteed to reduce to at least a certain fraction $L$ of the previous step, you have linear convergence: $$E_{k+1} \le L E_k.$$ This ...


8

The conjugate gradient algorithm works for semidefinite problems and produces the minimal norm solution.


8

Due to the under-resolved boundary layer near the lid, the velocity in the adjacent cells is significantly lower than the lid. This section is showing you a trick to make the code run faster while still being stable. Increasing the Courant number would normally make the method unstable, but since the velocity in all interior cells is significantly less than ...


8

Let's examine the one-dimensional three-point stencil case in detail, because I think it's important to be clear just how this behaviour arises, and what it means to set a point to a certain value in a finite-difference grid when the underlying function is discontinuous. The equation will be $$ u''(x) = \rho(x). $$ Instead of using the interval $[-1,1]$ with ...


8

Yes: See Higham's book "Accuracy And Stability of Numerical Algorithms", second edition, chapter 25: Nonlinear Systems and Newton's Method. In particular, see the section on the "limiting residual" in terms of a condition number for the Jacobian. It may well be that since your system is ill-conditioned, that you quickly hit the limiting residual and your ...


7

The Hartree-Fock equations are the result of performing constrained Newton-Raphson minimization of the energy with respect to the parameter space of Slater determinants (I don't have my copy of Szabo-Ostlund at hand, but I believe this is pointed out in the derivation). Hence, HF-SCF will converge if your starting guess is in a convex region around a minimum....


7

One cannot conclude from a residual how accurate the solution is. Between the best and the worst case in norm, there is a factor of exactly the condition number. More precisely, if the residual norm is r and the error norm is e then $\|A\|^{-1}\le e/r \le \|A^{-1}\|$, and both bounds are attainable. Taking the quotient of the bounds proves the claim. The ...


7

Convergence is a statement about the asymptotic limit $h\to0$. Therefore, it's not possible to "prove" convergence by any finite number of computations. The best you can do with computations is motivate the belief that an algorithm does converge. To prove convergence, you need to...write a proof.


7

If analytic techniques are disallowed but the periodic structure is known, here is one approach. Let $$g(x) = \frac{\cos x}{2-\cos x}$$ be periodic with period $2 \pi$, so that $$g(x) = \sum_j w_j e^{ijx}$$ where $$w_j = \frac{1}{2\pi}\int_0^{2\pi} g(x) e^{-ijx} dx$$ Thus, $$\begin{aligned} f(x) &= \sum_{k \ge 1} \frac{g(kx)}{k^p} \\ &= \sum_{k ...


7

The Taylor-Hood approximation of the Stokes flow is a mixed finite element method, for which error estimates generally have the form $$ \|u-u_h\|_V + \|p-p_h\|_M \leq C (\inf_{w_h\in V_h}\|u-w_h\|_V + \inf_{q_h\in M_h}\|p-q_h\|_M), \tag{1} $$ where $(u,p)\in V\times M$ is the exact solution and $(u_h,p_h)\in V_h\times M_h$ is the approximation. For the ...


7

If your matrix is symmetric, positive definite, the CG method may converge slowly, but it converges for $n\to\infty$. The only reason it does not converge on a computer are round-off errors, in particular if the condition number of the matrix, the quotient of largest and smallest eigenvalue is large. Experience is, that with double precision arithmetic, ...


7

The two normalization formulas result in two different algorithms, that they both "normalize" a vector is not so relevant. As an example, consider the following transition matrix: $$M = \begin{pmatrix}0&1\\1&0\end{pmatrix}.$$ Starting with $r=(1,0)$, consecutive normalized vectors $r$ will be $(0,1)$ and $(1,0)$, so the algorithm clearly will not ...


7

The issues you're running into now are not a failing of Newton-Raphson, but a question of coupling. You're doing iterated sequential coupling -- solving each equation sequentially and then iterating until (hopeful) convergence. No solver choice in place of NR is going to fix this lack of convergence, as long as you are doing iterated sequential coupling. ...


7

Standard terminology in nonlinear (i.e., derivative-based) optimization is that "global convergence" means "convergence to a stationary point no matter where you start from" (where the limit may in fact depend on the starting point). This is in contrast to local convergence, which requires that you start sufficiently close to a stationary point; if your ...


7

Unfortunately, convergence of GMRES does not have a clear dependence on the distribution of eigenvalues. It was proved by Greenbaum, Ptak and Strakos in 1996 that you can construct examples with an arbitrary spectrum and an arbitrary convergence history: that is, give me any $n$ nonzero complex numbers, and any decreasing sequence $\|r_k\|$, and I can ...


6

The interpretation is qualitatively correct. Note that linear and quadratic convergence are with respect to the worst case, the situation in a particular algorithm can be better than what you get from the worst case analysis given by Wolfgang Bangerth, though the qualitative situation usually corresponds to this analysis. In concrete algorithms (e.g., in ...


6

For many partial differential equations arising in nature, particularly with strong nonlinearities or anisotropies, the choice of an appropriate preconditioner can have a large effect on whether the iterative method converges rapidly, slowly, or not at all. Examples of problems that are known to have fast and effective preconditioners include strongly ...


Only top voted, non community-wiki answers of a minimum length are eligible