implicit midpoint method for newton's laws

I've been banging my head against the way trying to understand the Wikipedia article for the implicit midpoint method. I find the notation confusing which isn't helping. My goal is to use the implicit midpoint method to time advance a physics simulation based on Newton's Laws.

Given $$y(t) = (x(t),v(t)) = y_{n}$$ and $$\dot{y}(t) = (\dot{x}(t),\dot{v}(t)) = (v(t),a(t)) = (v(t),\frac{F(x(t))}{m}) = \dot{y}_{n}$$

Assuming I know $x(0),v(0),F(0)$ I would like to use the implicit midpoint method to find $y(t+h)=y_{n+1}$

As far as I can understand the article the update rule would be

$$y_{n+1} = y_{n} + \frac{h}{2}(\dot{y}_{n} + \dot{y}_{n+1})$$

but I'm not even sure that's correct because the way they write it in the article makes the fact that both sides of the equation depend on $y_{n+1}$ very explicit. Whereas the way I've written it it's not clear that is the case.

I was thinking I'd use Newton's method to solve this implicit system of equations which I think means that I need to do something like write

$$G(y_{n+1}) = y_{n+1} - y_{n} - \frac{h}{2}(\dot{y}_{n} + \dot{y}_{n+1}) = 0$$

but I really don't think that's correct because for Newton's method to work I need to be able to guess values for $y_{n+1}$ and evaluate $G$ at those values. However I need a value for $\dot{y}_{n+1}$ to evaluate $G$.

EDIT:

It may be useful to show how I ended up with this form for the midpoint method. My understanding of the notation $f(t,y) = \dot{y}(t)$ is that the $t$ parameter of $f$ only matters if there is an explicit dependence on time. For example if $\dot{y}(t) = t + y(t)$ then

$$f (t + \frac{h}{2},\frac{1}{2}(y(t) + y(t+h)) = t + \frac{h}{2} + \frac{1}{2}(y(t) + y(t+h))$$

Since there is no explicit dependence of $t$ in these equations I dropped the $t + \frac{h}{2}$ term and wrote

$$f(t + \frac{h}{2},\frac{1}{2}(y(t) + y(t+h)) = \dot{y}(t) = \frac{1}{2}(\dot{y}(t) + \dot{y}(t+h))$$

Thank you for taking the time to read my question.

• I think I realised something. By guessing $y_{n+1}$ I get the first component of $\dot{y}_{n+1}$. The second component can be computed since $F$ is a function of $y_{n+1}$ only. – Adam Sturge Jun 29 '17 at 1:10
• You are aware that you are using the implicit trapezoidal method? The midpoint method uses the derivative approximation at the midpoint at $t+h/2$. – LutzL Jun 29 '17 at 5:41
• @Lutzl No actually I didn't know that. I wasn't sure how to put in t+h/2 into the equations. Like what is $f(t+\frac{h}{2},\frac{1}{2}(y_{n}+ y_{n+1}))$ as it's described on the wiki page for the midpoint method?. Given that $y_{n}=y(t)$ and $y_{n+1} = y(t+h)$. Is it $\frac{1}{2} (\dot{y}(t+\frac{h}{2}) + \dot{y}(t+\frac{h}{2} + h))$? – Adam Sturge Jun 29 '17 at 16:37

For the midpoint method you approximate \begin{align} \frac{y(t+h)-y(t)}h&\approx\dot y(t+\tfrac h2) \\ &=f(t+\tfrac h2,y(t+\tfrac h2)) \\ &\approx f(t+\tfrac h2,\tfrac12(y(t+h)+y(t))). \end{align} All approximations with an error $O(h^2)$
This gives the method $$y_{k+1}=y_k+hf(t+\tfrac h2,\tfrac12(y_k+y_{k+1}))$$ This can be decomposed into first one step of the implicit Euler method and then one of the explicit Euler method, both with half the step size \begin{align} \text{implicit Euler: }&&y_{k+\frac12}&=y_k+\tfrac h2 f(t+\tfrac h2,y_{k+\frac12})\\ \text{explicit Euler: }&&y_{k+1}&=y_{k+\frac12}+\tfrac h2 f(t+\tfrac h2,y_{k+\frac12}) \end{align} where $(t+\tfrac h2,y_{k+\frac12})$ is the midpoint of the segment connecting the iteration points with $y_{k+\frac12}=\frac12(y_k+y_{k+1})$
• Name $z=\frac12(y_k+y_{k+1})$, then $y_{k+1}=2z-y_k$, and insert into the midpoint formula, $$2z-y_k=y_k+hf(t+\frac h2,z)$$ which transforms into $$z=y_k+\frac h2f(t+\frac h2,z)$$ which is the formula for an implicit Euler step. So you get to solve $$0=G(z)=z-y_k-\frac h2f(t+\frac h2,z)$$ either via Newton or via fixed-point iteration. – LutzL Jun 30 '17 at 15:46
• Also assuming that I've run Newton's Method to find $z$, and that I know $y_{k}$ from the start of the timestep, can't I just analytically solve for $y_{k+1}$ using $y_{k+1}=2z−y_{k}$? – Adam Sturge Jun 30 '17 at 17:47
• Yes, but you will have evaluated $s=f(t+h/2,z)$ anyway, so that you can just compute $y_{k+1}=y_k+h·s$ avoiding cancellation errors. – LutzL Jun 30 '17 at 18:05