-
Notifications
You must be signed in to change notification settings - Fork 63
/
rk4.m
30 lines (26 loc) · 865 Bytes
/
rk4.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function [t,u] = rk4(dudt,tspan,u0,n)
% RK4 Fourth-order Runge-Kutta for an IVP.
% Input:
% dudt defines f in u'(t)=f(t,u) (function)
% tspan endpoints of time interval (2-vector)
% u0 initial value (vector, length m)
% n number of time steps (integer)
% Output:
% t selected nodes (vector, length n+1)
% u solution values (array, (n+1) by m)
% Define time discretization.
a = tspan(1); b = tspan(2);
h = (b-a)/n;
t = a + (0:n)'*h;
% Initialize solution array.
u = zeros(length(u0),n+1);
u(:,1) = u0(:);
% Time stepping.
for i = 1:n
k1 = h*dudt( t(i), u(:,i) );
k2 = h*dudt( t(i)+h/2, u(:,i)+k1/2 );
k3 = h*dudt( t(i)+h/2, u(:,i)+k2/2 );
k4 = h*dudt( t(i)+h, u(:,i)+k3 );
u(:,i+1) = u(:,i) + (k1 + 2*(k2 + k3) + k4)/6;
end
u = u.'; % conform to MATLAB output convention