Prev Next Index-> contents reference index search external Up-> pycppad runge_kutta_4 runge_kutta_4_correct.py runge_kutta_4_correct.py

runge_kutta_4 A Correctness Example and Test

Discussion
Source Code

Discussion
Define $y : \B{R} \rightarrow \B{R}^n$ by $$y_j (t) = t^{j+1}$$ It follows that the derivative of $y(t)$ satisfies the runge_kutta_4 ODE equation where $y(0) = 0$ and $f(t, y)$ is given by $$f(t , y)_j = y_j '(t) = \left\{ \begin{array}{ll} 1 & {\; \rm if \;} j = 0 \\ (j+1) y_{j-1} (t) & {\; \rm otherwise } \end{array} \right.$$

Source Code

def fun(t , y) :
n        = y.size
f        = numpy.zeros(n)
f[0]     = 1.
index    = numpy.array( range(n-1) ) + 1
f[index] = (index + 1) * y[index-1]
return f
n  = 5              # size of y(t) (order of method plus 1)
ti = 0.             # initial time
dt = 2.             # a very large time step size to test correctness
yi = numpy.zeros(n) # initial value for y(t); i.e., y(0)

# take one 4-th order Runge-Kutta integration step of size dt
yf = runge_kutta_4(fun, ti, yi, dt)

# check the results
t_jp = 1.                                # t^0 at t = dt
for j in range(n-1) :
t_jp = t_jp * dt                    # t^(j+1) at t = dt
assert abs( yf[j] - t_jp ) < 1e-10  # check yf[j] = t^(j+1)

Input File: example/runge_kutta_4_correct.py