RK4

Propagates the state vector forward one time step using the (classic) Runge-Kutta fourth-order method.

Syntax

y_next = RK4(f,t,y,h)


Description

y_next = RK4(f,t,y,h) returns the state vector at the next sample time, y_next, given the current state vector y at time t, the function f(t,y) defining the ODE , and the step size h.

Input/Output Parameters

 Variable Symbol Description Format Input f $\inline&space;\mathbf{f}(t,\mathbf{y})$ multivariate, vector-valued function ($\inline&space;\mathbf{f}:\mathbb{R}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$) defining the ordinary differential equation $\inline&space;d\mathbf{y}/dt=\mathbf{f}(t,\mathbf{y})$ - inputs to f are the current time (t, 1×1 double) and the current state vector (y, p×1 double) - output of f is the state vector derivative (dydt, p×1 double) at the current time/state 1×1function_handle t $\inline&space;t_{n}$ current sample time 1×1double y $\inline&space;\mathbf{y}_{n}=\mathbf{y}(t_{n})$ state vector (i.e. solution) at the current sample time p×1double h $\inline&space;h$ step size 1×1double Output y_next $\mathbf{y}_{n+1}=\mathbf{y}(t_{n+1})$ state vector (i.e. solution) at the next sample time, $t_{n+1}=t_{n}+h$ p×1double

Note

• This documentation is written specifically for the case of vector-valued ODEs. However, this function can also be used for matrix-valued ODEs of the form , where .

Example

Consider the initial value problem

Find the solution until using RK4. Then, compare your result to the solution found by solve_ivp using the classic RK4 method.

First, let's define our ODE () and initial condition in MATLAB.

f = @(t,y) y;
y2 = 3;


Let's define a time vector between and with a spacing of .

h = 0.01;
t = (2:h:10)';


Solving for using RK4 and comparing the result to the result obtained using solve_ivp with the classic RK4 method,

% preallocate vector to store solution
y = zeros(size(t));

% store initial condition
y(1) = y2;

% solving using "RK4"
for i = 1:(length(t)-1)
y(i+1) = RK4(f,t(i),y(i),h);
end

% solving using "solve_ivp"
[t_ivp,y_ivp] = solve_ivp(f,[2,10],y2,h,'RK4');

% maximum absolute error between the two results
max(abs(y_ivp-y))

ans =

1.533408067189157e-09



As expected, the two methods obtain identical results.