I am trying to solve a linear system of ODEs of the form: $$ \frac{du}{dt} = A u, \quad u(0)=k$$ where $A$ is a 2x2 matrix and $u(t)$ is a 2x1 column vector. I want to solve this numerically, using Euler forward scheme, which means that I will use finite differences to approximate the above equation this way: $$ \frac{u[n+1]-u[n]}{dt} = Au[n],$$ using $u[0] = k$ as the initial condition. This means that the iteration matrix of the Euler method is given by: $$ u[n+1] = Bu[n], \quad B = I + A*dt, $$ where $I$ is the 2x2 identity matrix. I am having trouble studying the 'stability' of Euler forward in this context. For me, 'stability' means that the difference between the numerical solution and the exact analitical solution doesn't increase over time. I know that Lax equivalence theorem states that, sometimes, stability+consistency = convergence. So, I was wondering wether the following condition: $$ \text{Iteration method is convergent i.e. } \rho(B) < 1,$$ is necessary and sufficient to guarantee that the Euler forward scheme is 'stable', in the sense described above. In other words, does 'convergence' imply 'stability'? Is the 'stability' in Lax theorem the same 'stability' I am refering to? For me 'convergence' of an iterative mehod is quivalent to $\rho <1$, but is this 'convergence' the same as the one in Lax theorem? Moreover, does Lax theorem apply in this context, or is there some reason why the theorem cannot be applied here?


1 Answer 1


The solution of $\frac{du}{dt} = Au$ is $u(t) = \exp(tA)u(0)$, and explicit Euler approximates $\exp(tA)$ using $\lim_{n\to\infty} \left(I+\frac{t}{n}A\right)^n$. Of course in practice you cannot compute this for $n\to \infty$ so you choose a finite $n$. Then the explicit Euler step size is really $\tau = t/n$. But if $\tau$ is too large then $\|I+\tau A\| > 1$ and any initial error will explode. If you consider the $2$-norm and $A$ is symmetric negative semi-definite, then you need to choose $\tau \in \left(0, \frac{2}{\rho(A)}\right)$ where $\rho(A)$ is the spectral radius of $A$ (maximum absolute eigenvalue).

Of course, for a $2\times 2$ matrix you can compute the solution directly using the eigendecomposition, then you don't need any time stepping scheme.

Edit: A clarification for the stability criterion. Let $A$ be a symmetric negative semi-definite matrix (in the complex case it is sufficient that it is normal negative semi-definite afaik). Since it's symmetric and real it always has an eigendecomposition $A=QDQ^T$. You want for any error $e$ to not increase with the iterations, i.e. $\|e_{i+1}\|<\|e_i\|$, but $$\|e_{i+1}\| = \|(I+\tau A) e_i\| \leq \|I+\tau A\| \|e_i\|.$$ Then $\forall e_i, \,\|e_{i+1}\|\leq \|e_i\| \iff \|I+\tau A\|\leq 1$.

If you pick the $2$-norm the criterion $\|I+\tau A\|_2\leq 1$ is equivalent to $\|Q(I+\tau D)Q^T\|_2 = \|I+\tau D\|_2\leq 1$. Thus you get the condition $|1+\tau \lambda_i|\leq 1$ for every eigenvalue $\lambda_i$. For $1+\tau \lambda_i\geq 0$ you get $1+\tau\lambda_i \leq 1$ which becomes $\tau \geq 0$ (taking into account that $\lambda_i\leq 0$). For $1+\tau\lambda_i <0$ you get $-1\leq 1+\tau \lambda_i$ which becomes $\tau \leq -\frac{2}{\lambda_i}$. Note that if you have an eigenvalue $\lambda_k$ which is positive then you would get a condition $\tau\leq 0$ but from the negative ones you have $\tau \geq 0$ so this process is stable (i.e. does not magnify any components of the error) only if $A$ is not indefinite. In practice this is usually the case, e.g. $A$ approximating the Laplacian $A\approx \Delta$. In either case, if for all eigenvalues you have $\lambda_i\leq 0$ then you get the constraints $\tau \in \left[0, \frac{2}{|\lambda_i|}\right]$. To satisfy all of those you need to take the intersection of these constraints $\tau \in \left[0, \frac{2}{\max_i |\lambda_{i}|}\right]$. Finally you have to omit $\tau=0$ as you cannot have $\tau = t/n=0$ unless $t=0$ (and a step with $\tau=0$ would change nothing anyways). You should also omit $\tau = \frac{2}{\rho(A)}$ as it does not guarantee an error decrease and you can get an oscillating solution (i.e. typically you want a strict inequality $\|I-\tau A\|<1$ so that $\|e_{i+1}\| < \|e_i\|$).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.