Is implicit Euler forward integration necessarily more accurate than the explicit method? Currently, I'm trying to use the Van Der Pol osciallator to compare the integration accuracy between the implicit Euler forward and explicit Euler Forward. 
The explicit Euler method is: $x_{k+1} = x_{k} + hf(x_{k})$;
The implicit Euler method is: $x_{k+1} = x_{k} + hf(x_{k+1})$; 
where $f(x)$ is the ODE of the Van Der Pol oscillator and $h$ is the integration step. I changed the integration step and compared the integration results with the ode15s from scipy. The formula to calculate the error for the implicit and explicit method is the same. For example, the implicit method error could calculated by the following, 
np.linalg.norm(ieu_sol - ref_sol)/np.linalg.norm(ref_sol)

It looks like the implicit method is always more accurate than the explicit method (ignore the sign). My questions are the following: 
1) is it always true that the implicit method is more accurate than the explicit method?
2) How would the ODE change the accuracy? (in this example I tried the Van Der Pol oscillator but what about other dynamic systems?)
3) Are there any references you may suggest reading?
Thank you.
 A: *

*No. When there is sufficient differentiability available the global error satisfies an asymptotic error expansion of the form $$x(t) - y_h(t) = \alpha(t) h^p + \beta(t) h^q + O(h^r), \quad h \rightarrow 0, \quad h > 0,$$ where $x$ is the exact solution, $y_h$ is the computed approximation and the functions $\alpha$ and $\beta$ can be computed by solving differential equations which depend on exact solution and the scheme deployed to compute $y_h$. It is possible to engineer situations where either the explicit or the implicit Euler is superior to the other.

*The values of the coefficients of the asymptotic error expansion given above depends on the solution of the differential equation.

*The existence of the asymptotic error expansion given above is discussed in Hairer et. al "Solving Ordinary Differential Equations I: Nonstiff Problems", section II.8.

A: In general no, as the implicit method, if considered backward in time, is the explicit method.
However, the implicit method has a larger region of stability. Which means that when the exact solution is moving towards some attractor, then the implicit method is more faithful to that behavior than the explicit method. This effect will be more pronounced for larger step sizes, for very small step sizes the difference in error should reduce. 
