Euler, semi-implicit Euler, and RK4

Learning outcomes

Implement three time-stepping methods.
Use stability and convergence tests to explain numerical behavior.
Distinguish local accuracy from long-term energy behavior.

Prerequisites: Linear algebra, ODEs, and earlier lectures in this sequence.

Explicit Euler

For $\dot x=f(x)$ , explicit Euler uses $x_{k+1}=x_k+\Delta t f(x_k)$ . Global error is first order under the usual smoothness and fixed-time assumptions. For $\dot x=\lambda x$ , the amplification factor is $1+\Delta t\lambda$ , so a decaying continuous system remains stable only inside the method's stability region.

For an undamped oscillator with frequency $\omega$ , explicit Euler has eigenvalue magnitudes $\sqrt{1+(\omega\Delta t)^2}>1$ . Discrete energy grows while physical energy stays constant.

Semi-implicit Euler

For a separable mechanical system with position-dependent acceleration,

v_{k+1}=v_k+\Delta t\,a(q_k),\qquad q_{k+1}=q_k+\Delta t\,v_{k+1}.

This symplectic Euler variant is also first order, but can keep oscillator energy bounded with oscillatory error at a suitable timestep. Generic bounded oscillator motion requires

\omega\Delta t<2

. At or beyond the boundary, do not expect a reliable trajectory.

With velocity-dependent forces, an engine may treat selected terms implicitly. MuJoCo's integrator named Euler therefore differs in detail from a hand-written explicit Euler method.

Runge–Kutta 4

RK4 samples four slopes:

k_1=f(x),\quad k_2=f(x+\tfrac{\Delta t}{2}k_1),\quad k_3=f(x+\tfrac{\Delta t}{2}k_2),\quad k_4=f(x+\Delta t k_3),

then combines them with weights

(1,2,2,1)/6

. For smooth problems it has fourth-order global accuracy over a fixed interval. Four evaluations cost more than one. RK4 is not generally symplectic and does not exactly conserve energy.

Integrate the same smooth problem to the same final time with $\Delta t$ , $\Delta t/2$ , and $\Delta t/4$ . In the asymptotic regime, Euler's error roughly halves and RK4's decreases by about sixteen. Compare to an analytic solution when available. Comparing runs at different final times invalidates the ratio.

The practice uses a unit oscillator with $q(0)=1$ , $v(0)=0$ and exact solution $(\cos t,-\sin t)$ . It reports trajectory error and maximum relative energy error.

Practical limits

Contact events can invalidate smooth-order expectations. A stiff spring may force small explicit timesteps even when visible motion is slow. Solver tolerances and floating-point error eventually dominate refinement.

Exercise: Which method is best for a 10-second smooth prediction and which for a very long Hamiltonian simulation? Answer guide: compare accuracy per computation for the first. For the second, investigate geometric structure and long-term energy behavior. Formal order alone does not select a universal winner. Variational integrators appear in the optional post-course extension.