Calculate the numeric solution at each step to the ODE f(x, y) using RK4 https://en.wikipedia.org/wiki/Runge-Kutta_methods Arguments: f -- The ode as a function of x and y y0 -- the initial value for y x0 -- the initial value for x h -- the stepsize x_end -- the en
(f, y0, x0, h, x_end)
| 2 | |
| 3 | |
| 4 | def runge_kutta(f, y0, x0, h, x_end): |
| 5 | """ |
| 6 | Calculate the numeric solution at each step to the ODE f(x, y) using RK4 |
| 7 | |
| 8 | https://en.wikipedia.org/wiki/Runge-Kutta_methods |
| 9 | |
| 10 | Arguments: |
| 11 | f -- The ode as a function of x and y |
| 12 | y0 -- the initial value for y |
| 13 | x0 -- the initial value for x |
| 14 | h -- the stepsize |
| 15 | x_end -- the end value for x |
| 16 | |
| 17 | >>> # the exact solution is math.exp(x) |
| 18 | >>> def f(x, y): |
| 19 | ... return y |
| 20 | >>> y0 = 1 |
| 21 | >>> y = runge_kutta(f, y0, 0.0, 0.01, 5) |
| 22 | >>> float(y[-1]) |
| 23 | 148.41315904125113 |
| 24 | """ |
| 25 | n = int(np.ceil((x_end - x0) / h)) |
| 26 | y = np.zeros((n + 1,)) |
| 27 | y[0] = y0 |
| 28 | x = x0 |
| 29 | |
| 30 | for k in range(n): |
| 31 | k1 = f(x, y[k]) |
| 32 | k2 = f(x + 0.5 * h, y[k] + 0.5 * h * k1) |
| 33 | k3 = f(x + 0.5 * h, y[k] + 0.5 * h * k2) |
| 34 | k4 = f(x + h, y[k] + h * k3) |
| 35 | y[k + 1] = y[k] + (1 / 6) * h * (k1 + 2 * k2 + 2 * k3 + k4) |
| 36 | x += h |
| 37 | |
| 38 | return y |
| 39 | |
| 40 | |
| 41 | if __name__ == "__main__": |