[Math] Calculating stability and order of implicit midpoint scheme

numerical methodsordinary differential equations

Consider solving $y'(t) = f(t,y(t))$ by the implicit midpoint method:
$$
y_{n+1} = y_n + h \cdot f \left(t_n + \frac{h}{2},\frac{y(t_n) + y(t_{n+1})}{2} \right).
$$
I want to determine the order and regions of stability for this method.

  1. Order. My first idea was to do the following: Substitute the exact solution.
    \begin{align*}
    &y(t_{n+1}) – \left[ y(t_n) + h \cdot f \left(t_n + \frac{h}{2},\frac{y(t_n) + y(t_{n+1})}{2} \right) \right] \\
    =& [ y(t_n) + hy '(t_n) + \frac{1}{2} h^2 y''(t_n) + O(h^3) ] \\
    – & [ y(t_n) + h \left[ y'(t_n) + \frac{1}{2} h y''(t_n) + O(h^2) \right]
    \end{align*}
    which is $O(h^3)$. But I don't think I can substitute in for the last part $f(t,y(t)) = y'(t)$ since $\frac{y(t_n) + y(t_{n+1})}{2} \neq y \left( t_n + \frac{h}{2} \right)$ . . . It's only an approximation.

  2. Stability. Apply to the test problem $y'(t) = f(t,y(t)) = \lambda y(t)$. This is done here (Determine a stability region?), but again $f(\cdot, y(\cdot))$ is not of the form $f(t, y(t))$, so I don't know why the solution is valid. In particular, why is the test problem $y' =f(t,y(t)) = \lambda y(t)$ applied to $f(t_{n+1/2} , (y_n + y_{n+1})/2)$ equal to $\lambda (y_n + y_{n+1})/2$?

The questions are similar and I probably have some misconception on numerics that is (hopefully) easy to clarify.

EDIT: Should I just think about $f \left(t_n + \frac{h}{2},\frac{y(t_n) + y(t_{n+1})}{2} \right)$ as $f(t_{n+1/2} , y_{n+1/2})$? If so I think my question is answered, but I would appreciate someone wiser in the field taking a look.

Best Answer

First a correction: The implicit midpoint method proceeds by solving in each step the implicit equation $$ y_{n+1} = y_n + h \cdot f \left(t_n + \frac{h}{2},\frac{y_n + y_{n+1}}{2} \right). $$ One gets the local truncation error by replacing $y_k$ with $y(t_k)$ for an exact solution of $y'(t)=f(t,y(t))$ that is close to the numerical solution, for instance satisfies $y(t_n)=y_n$ or $y(t_{n+1})=y_{n+1}$, and then computing the defect in the above equation.

To keep the symmetry of the situation, insert an exact solution with $y(t_{n+\frac12})=y_{n+\frac12}$. Then with Taylor expansions \begin{align} y(t_{n+1})-y(t_n)&-h·f\left(t_n + \frac{h}{2},\frac{y(t_n) + y(t_{n+1})}{2} \right) \\ &=y\left(t_{n+\frac12}+\frac{h}2\right)-y\left(t_{n+\frac12}-\frac{h}2\right) \\&\qquad-h·f\left(t_{n+\frac12},y(t_{n+\frac12})+\frac{y(t_n) - 2y(t_{n+\frac12})+ y(t_{n+1})}{2} \right) \\ &=h·y'(t_{n+\frac12})+\frac{h^3}{24}y'''(t_{n+\frac12})+O(h^5) \\&\qquad-h·f\left(t_{n+\frac12},y(t_{n+\frac12})\right)-\frac{h^3}{8}·\partial_yf\left(t_{n+\frac12},y(t_{n+\frac12})\right)y''(t_{n+\frac12})+O(h^5) \\ &=\frac{h^3}{24}\left(y'''(t_{n+\frac12})-3\partial_yf\left(t_{n+\frac12},y(t_{n+\frac12})\right)y''(t_{n+\frac12})\right)+O(h^5) \end{align} This confirms the local $O(h^3)$ truncation error.


For the test problem, just insert the function $f(t,y)=λy$ into the method step to get $$ y_{n+1} = y_n + h \cdot λ \left(\frac{y_n + y_{n+1}}{2} \right), \\~\\ \left(1-\frac{hλ}{2}\right)y_{n+1}=\left(1+\frac{hλ}{2}\right)y_n \\~\\ y_{n+1} = \left(\frac{1+\frac{hλ}{2}}{1-\frac{hλ}{2}}\right)y_n $$ and then check where $$\left|\frac{2+hλ}{2-hλ}\right|\le 1.$$

Related Question