# Unexpectedly Slow Convergence Implicit Euler

I'm solving the coupled ODE $$\left[\begin{array}{c}x^\prime(z)\\p_x^\prime(z)\end{array}\right] = C(z)\cdot\left[\begin{array}{c}x(z)\\p_x(z)\end{array}\right] = \left[\begin{array}{cc}0& A(z)\\B(z) & 0\end{array} \right]\cdot\left[\begin{array}{c}x(z)\\p_x(z)\end{array}\right]$$ using the implicit Euler method to get the transfer matrix $$M$$ of the system where $$\left[\begin{array}{c}x(z)\\p_x(z)\end{array}\right] = M(z)\cdot\left[\begin{array}{c}x(0)\\p_x(0)\end{array}\right]$$ The problem I'm running into is that some of the elements of $$M$$ converge at a surprisingly slow rate. I calculated the matrices for a variety of number of steps in the region $$n$$ and fit the region of small $$1/n$$ to a power law with a constant added to estimate relative error. Elements $$M_{11}$$ and $$M_{21}$$ have an exponent of ~1.5 while the others are more like 0.5. I would have expected that since the error in one step should behave like O($$1/n^2$$) and there are $$n$$ steps that the convergence should be at least O($$1/n$$) not O($$\sqrt{1/n}$$).

In my implementation, I have manually calculated the matrices $$(1-C(z_i))^{-1}$$ at every step and then repeatedly taken matrix products to find the final result.

Is there something that I'm missing that would give me this type of convergence, or is this just something expected from certain systems. For reference, my matrix of coefficients looks like this: There is a regular singular point at $$z=0$$ which gets skipped over due to implicit Euler only looking at the points $$z_{i+1}$$.

• Relative error with respect to what? Analytical solution? – Alone Programmer Feb 14 '20 at 0:15
• How do you compute $M(z)$? That involves a matrix exponential. – Wolfgang Bangerth Feb 14 '20 at 4:07
• The fact that either A or B is discontinuous seems like a likely culprit, but I don’t know for certain that it should give you a convergence rate of 1/2. – David Ketcheson Feb 14 '20 at 6:57
• I estimate relative error by fitting the function $\epsilon(h) = Ah^B + C$ to my data for small $h$. There isn't a good analytical solution for these cases. – ElectronsAndStuff Feb 14 '20 at 17:50
• $M(z)$ is computed from implicit Euler. For each time step you already have the matrices which transform $y(z)$ to $y(z + dz)$. You just need to multiply them all together. Alternatively, you can think of it as numerically solving the IVP for the initial conditions y = [1, 0] and y = [0, 1] separately and then because all solutions are a linear combo of that, you can pull out the matrix from the results. – ElectronsAndStuff Feb 14 '20 at 17:52

## 1 Answer

Due to your formulation, I call $$X(z) = \begin{bmatrix} x(z) \\ p_{x}(z) \end{bmatrix}$$ so your ODE is written in matrix form:

$$X^{'}(z) = C(z) X(z)$$

Where: $$C(z) = \begin{bmatrix} 0 & A(z) \\ B(z) & 0 \end{bmatrix}$$.

Your general formula by using backward Euler method is:

$$\frac{X(z+\Delta z) - X(z)}{\Delta z} = C(z+\Delta z) X(z+\Delta z)$$

So:

$$X(z+\Delta z) = (\mathbf{I}-\Delta z C(z+\Delta z))^{-1} X(z)$$

Where $$\mathbf{I} = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}$$. Also you have:

$$X(\Delta z) = (\mathbf{I}-\Delta z C(\Delta z))^{-1} X(0)$$ $$X(2\Delta z) = (\mathbf{I}-\Delta z C(2\Delta z))^{-1} X(\Delta z)$$

or:

$$X(n\Delta z) = (\mathbf{I}-\Delta z C(n\Delta z))^{-1} X((n-1)\Delta z)$$

or:

$$X(n\Delta z) = \prod_{i=1}^{n} (\mathbf{I} - \Delta z C(i\Delta z))^{-1} X(0)$$

So explicitly your $$M(n,\Delta z)$$ is:

$$M(n,\Delta z) = \prod_{i=1}^{n} (\mathbf{I} - \Delta z C(i\Delta z))^{-1}$$

As far as I understand, you are trying to plot convergence rate define as:

$$\varepsilon = \frac{||M(n,\Delta z) - M(n-1,\Delta z)||_{F}}{||M(n-1,\Delta z)||_{F}}$$

The operator $$||M(n,\Delta z)||_{F}$$ is just Frobenius norm.

We know that:

$$M(n,\Delta z) - M(n-1,\Delta z) = M(n-1,\Delta z) ((\mathbf{I} - \Delta z C(n\Delta z))^{-1}-\mathbf{I})$$

Also for Frobenius norm we have this inequality:

$$||AB||_{F} < ||A||_{F} ||B||_{F}$$

So:

$$||M(n,\Delta z) - M(n-1,\Delta z)||_{F} < ||M(n-1,\Delta z)||_{F} ||((\mathbf{I} - \Delta z C(n\Delta z))^{-1}-\mathbf{I})||_{F}$$

or:

$$\varepsilon < ||((\mathbf{I} - \Delta z C(n\Delta z))^{-1}-\mathbf{I})||_{F}$$

Now, you see that the behavior of this convergence rate strongly depends on behavior of $$A(z)$$ and $$B(z)$$, but as you know that Frobenius norm is just a Euclidean distance and is defined by square root and if $$(\mathbf{I} - \Delta z C(n\Delta z))^{-1}$$ is just is in the order of $$O(\frac{1}{n})$$, because of that square root your final convergence rate is in the order of $$O(\frac{1}{\sqrt{n}})$$. But for more accurate analysis you need to be more explicit about $$A(z)$$ and $$B(z)$$ I believe.

• If I understand you correctly, you're talking about the change in between the solution vectors $X(n\Delta z)$ and $X((n-1)\Delta z)$. The convergence rate I'm actually concerned with is of the family of solution vectors at a point $X(z) = M(z) \cdot X(0)$ as I change step size, not position. What I'm plotting is $(M_{ii}(h) - M_{ii}(0))/M_{ii}(0)$ for each element of the matrix and the step size $h$, but crucially for fixed final $z$. I'm approximating $M_{ii}(0)$ with a numerical solution with $h$ much smaller than what I study in my plots. Also, $n = z/h$ for the previous plot. – ElectronsAndStuff Feb 16 '20 at 21:19
• @ElectronsAndStuff So, basically you are plotting $M_{ij}(\Delta z)$ versus $\Delta z$? But, what's the unitless $n$ in your plot then? And I'm not sure what $M_{ij}(\Delta z)$ vs. $\Delta z$ suppose to tell you... – Alone Programmer Feb 16 '20 at 21:23
• Exactly, but I suppose in your notation I am also changing $n$ so that $n \Delta z = z$ remains constant. The setup is that I care about the solution to the IVP at some $z$. I'd like to know how bad the numerical error is as I change step step size so that I can choose a reasonable trade-off between accuracy and computational time. The term $n$ is the number of steps the solver takes and is just $n = z/\Delta z$. Edit: Sorry, I realized I have bad notation in my previous comment. By $M(h)$ in the relative error equation, I mean $M(z)$ as I change $h$, so it should be $M(z,h)$. – ElectronsAndStuff Feb 16 '20 at 21:31
• @ElectronsAndStuff OK, but you think $M_{ij}$s show you the numerical errors at that particular point? I don't think so. So, in summary: You don't have any analytical solution here cause obviously the solution of this system of ODEs can't be derived with hand, but you are interested in a solution of particular point to know how it behaves when you change $\Delta z$, am I right? If, I'm right, so why you are trying to monitor $M_{ij}$s? The more straightforward way is to just monitor $X$ at that particular point, which I call it $z^{*}$, and plot $X(z^{*})$ versus $\Delta z$, and you are done! – Alone Programmer Feb 16 '20 at 21:36
• The problem is that for my application, it's the $M_{ij}(z)$s that we care about, not $X(z)$ directly. Because of this, we are concerned with the numerical error of the matrix elements themselves, not just the solution vectors for a particular set of initial values. – ElectronsAndStuff Feb 16 '20 at 21:46