Second Order Differential Equation Solving Algorithm

Second order differential equation problems can be written as the following initial value problem:

• y' = f(t, x, y), 0 < t <= T
• y = x'
• x0 = x(0) [initial value]
• y0 = x'(0) [initial derivative]

Let x[i] = x(i * dt) and y[i] = y(i * dt) for i = 0 to i = N - 1, where dt = T / (N - 1). As with the first order ode problem, for dt small enough, we can use the forward difference approximation to obtain

y(t) = y(t - dt) + f(t - dt, x(t - dt), y(t - dt)) * dt

Again, we can use the forward difference approximation y(t) = x'(t) = [x(t + dt) - x(t)] / dt to get

x(t + dt) = x(t) + y(t) * dt

So, we can compute all values of x[i] and y[i] with the following:

x[i] = { x[i - 1] + y[i - 1] * dt   0 < i <= N - 1
       { x0                         i = 0

y[i] = { y[i - 1] + f(i * dt, x[i - 1], y[i - 1]) * dt    0 < i <= N - 1
       { y0                                               i = 0

Thus, the algorithm simplifies to the following:

x, y = np.zeros([2, N])
x[0], y[0] = x0, y0
for i in range(1, N):
  x[i] = x[i - 1] + y[i - 1] * dt
  y[i] = y[i - 1] + f((i - 1) * dt, x[i - 1], y[i - 1]) * dt

[full code implementation]