In lib490 implement and test the fourth order RungeKutta ODE

In lib490, implement and test the fourth order Runge-Kutta ODE integrator for an arbitrary number of coupled ordinary differential equations. Create a function to perform one time step. It should have the signature: runge_kutta_4(function, time, y0, dt, *args) It should return a list of updated y values, time is the current time in the integration, y0 is a list of current values for the dynamical variables, and dt is the size of the time step, function is a single callable that returns dy/dt for all dynamical variables. It is expected to have the call signature equations(time, y0, *args) and should return a list of dy/dt values, args is a list of optional arguments that can be passed to function. Create a unit test for your functions to show that it returns, within numerical tolerance, the correct Runge-Kutta 4 updated values for the position and velocity of the simple harmonic oscillator. The choice of initial conditions and time step is up to you but should be non-trivial (y notequalto (0, 0)). To obtain the exact values Runge-Kutta 4 for comparison, integrate the system one time step by hand.

Solution

a)

from math import sqrt

def rungaKutta(function, x0, y0, x1, n):
px = [0] * (n + 1)
py = [0] * (n + 1)
h_val = (x1 - x0) / float(n)
px[0] = x = x0
py[0] = y = y0
for i in range(1, n + 1):
m1 = h_val * function(x, y)
m2 = h_val * function(x + 0.5 * h_val, y + 0.5 * m1)
m3 = h_val * function(x + 0.5 * h_val, y + 0.5 * m2)
m4 = h_val * function(x + h_val, y + m3)
px[i] = x = x0 + i * h_val
py[i] = y = y + (m1 + m2 + m2 + m3 + m3 + m4) / 6
return px, py

def function(x, y):
return x * sqrt(y)

----------------------------------------------------------------------------------------------------------------------------------------------------

(b)

px, py = rungaKutta(function, 0, 1, 15, 200)
for i, j in list(zip(px, py))[::10]:
print(\"%4.1f %10.5f %+12.4e\" % (i, j, j - (4 + i * i)**2 / 16))