III Exploring the world of motions: a refined approach
Since the iteration method used up to now might lead, with convenient step sizes, to physically non-tolerable results, in particular for log times, let us try to improve the method. A possible criticism against approximation (5) is that the rate of change ($v_n$, $a_n$) is taken at the beginning of the next time interval of length $\Delta t$. It is clear that better approximation can be obtained by taking the rate of change in the middle, at half integer time steps. Let us therefore replace (5) by $$ v_{n+1} = v_n + \class{colorRed}{a_{n+1/2}} \Delta t, \quad x_{n+1} = x_n + \class{colorRed}{v_{n+1/2}} \Delta t, \qquad (6) $$ where $\class{colorRed}{a_{n+1/2}}$, $\class{colorRed}{v_{n+1/2}}$ denote the acceleration and the velocity, respectively, at the time instant $(n+1/2)\Delta t$. These are unknown yet. But why not to apply another pair of iterations? We can take the method used up to now, but for only half a time step, $\Delta t/2$, from which $$ \class{colorRed}{v_{n+1/2}} = v_n + a_n \Delta t/2, \quad x_{n+1/2} = x_n + v_n \Delta t/2 \qquad (7) $$ follows. The half-integer acceleration $\class{colorRed}{a_{n+1/2}}$ is then obtained from the acceleration rule $a(x,v)$ given the knowledge of $x_{n+1/2}$ and $v_{n+1/2}$.
Relations (6), (7) provide an improved numerical scheme (called mid-point method) which is more accurate than (5). Its one-step error turns out, in general, to be proportional to the third power of the time step: $\Delta t^3$.
As an illustration, apply the improved scheme to the $x$ coordinate of the vertical projection tackled in the
table of Section I, in the knowledge of $v(t)=10-10t$. Show that the one-step error $5 \Delta t^2$ found with
the earlier scheme disappears here, the exact recursion is recovered, indicating that the new method generates
errors of a higher order.
result.
The midpoint velocity is $v_{n+1/2} = 10-10(n\Delta t+\Delta /2)$. The increment taken with this velocity is $$ x_{n+1} = x_n+(10-10n\Delta t+5 \Delta t) \Delta t = x_n+(10-10n\Delta t) \Delta t + 5 \Delta t^2. $$ Let us determine the position from the exact expression $x(t)=10^t-5 t^2$ at the $n+1$-st and the $n$th instant: $$x((n+1)\Delta t)= 10n\Delta t + 10 \Delta t - 5(n^2 + 2n +1) \Delta t^2,$$ $$x(n\Delta t)= 10n\Delta t - 5 n^2 \Delta t^2.$$ The difference of the two is $10 \Delta t - (10n + 5) \Delta t^2.$ This exact one-step displacement is identical to the one obtained from the iteration. The one-step error is thus zero in this case. In this very special case no $\Delta t^3$ term is found, the method, however, generates order $\Delta t^3$ errors in general cases.
Motion on a spring revisited. Monitor the motion by means of the improved numerical scheme.
Motion on a spring from a new view. The improved accuracy enables us to turn to a novel view of the motion: the velocity vs position view (the so-called phase space view). We plot the $v_n$ values (column I) vs the $x_n$ values (column H). This view provides, in general, a compact pattern since time is eliminated.
Damped spring. Damping of springs is often linearly dependent on velocity. It always slows down the motion, the sign is negative, we choose the proportionality constant to be 0.2. Let us therefore represent this by considering the acceleration rule $$ a(x,v)=-x-0.2v. $$ Monitor the motion both in time and in the $x,v$ view.
Driven damped spring. When an external time-dependent periodic force also acts on the mass fixed to the spring, the motion is called driven. Taking the forcing to be a cosine function of time, the acceleration becomes $$ a(x,v,t)=-x-0.2v+\cos(0.8 t) $$ (the frequency of the driving is 0.8 here, and its period is $2 \pi/0.8=7.85$ time units). This is a case when the acceleration depends also explicitly on time, and therefore in $a_{n+1/2}$ the time should appear with its current value, i.e. as a half integer multiple of the time step.
Repelling-attracting spring. Let us change to a slightly more complicated force in which a cubic term also occurs (without damping and driving in this first case): $$ a(x)=x-x^3. $$ The corresponding force is repulsive and proprotional to x for small distances, while in larger distances attraction occurs: $-x^3$ dominates, for large $x$. Study the motion.
Damped repelling-attracting spring. With a linearly velocity dependent damping (as for the linear spring) the resultant acceleration is $$ a(x,v)=x-x^3-0.2 v.$$ Study the motion.
Revolution of the Earth around the Sun.
As an example for a planar motion monitored by the improved method, let us consider the revolution of the Earth due
to the gravitational force of the Sun whose mass is denoted by $M$. The magnitude of the force is $\gamma Mm/r^2$ with
$r=(x^2+y^2)^{1/2}$ denoting the distance from the Sun. Since the force points in the direction of the Sun, the vector
form of the force is $-\gamma Mm/r^2\ \boldsymbol{r}/r$. For simplicity we take $\gamma M=1$, the two components of the
acceleration are thus
$$ a_x(x,y)=-x/(x^2+y^2)^{3/2}, \quad a_y(x,y)=-y/(x^2+y^2)^{3/2}. $$
Explore the motion with the initial conditions $x_0=0.5,\ y_0=0,\ v_{x,0}=0,\ v_{y,0}=1.63.$