I have created a midpoint algorithm to solve 2nd order ODEs in Matlab. And now Im comparing my solver with built in - ode23s. I have used a harmonic motion described as a 2nd order ODE, the .m file with the ODE looks like this

function dy=oscillator(t,y)


And after plotting through ode23s and my Midpoint I got this graph. Now my question is how is it possible that using the modified euler method(midpoint) to solve the oscillator problem does not give me the same result as using ode23s?

  • 1
    $\begingroup$ Looks like a stability issue. $\endgroup$ – Eff May 15 '17 at 16:35

Your ODE seems to be the system $\dot{a} = b$, $\dot{b} = -a$.

Looking at $a^2 + b^2$, we have $(a^2 + b^2)^\bullet = 2 a \dot{a} + 2 b \dot{b} = 0$.

For simplicity, let's take the explicit Euler method (the underlying problem is the same with the explicit(?) midpoint rule you are using), and check whether it respects the invariant.

For the explicit Euler method we have $a_{n+1} = a_n + h b_n$, $b_{n+1} = b_n - h a_n$, whence $$ a_{n+1}^2 + b_{n+1}^2 = (1 + h^2) (a_n^2 + b_n^2) . $$

Thus the magnitude of $a_n$ and $b_n$ increases over the long run, which is what you observe.

Suitable numerical integrators for this kind of system should be the symplectic ones.

Note though that the numerical method is stable: small perturbations in the input result in relatively small perturbations in the output (even if they are amplified by a factor that depends on the time $t$ where you look).

  • $\begingroup$ Thank you kindly for your clean explanation! $\endgroup$ – Denis Benka May 15 '17 at 16:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.