Runge-Kutta methods are one of the most widely used methods for numerically solving differential equations. In this post, I will explain how to implement the vectorised code for solving coupled ODEs using fourth order RK methods in GNU-Octave. The same code can be easily translated to matlab with minimal changes.
The explanation for the code is provided below. Note that this is not a 100% vectorised implementation. RK methods are, of course, iterative algorithms. The aim is to eliminate unnecessary loops through vectorisation. If this code can be improved further, please let me know in the comments section below.
RK4 Method
Consider a set of coupled, first order, Ordinary Differential Equations of the general form,
dtdxβ=fxβ(t,x,y,z)
dtdyβ=fyβ(t,x,y,z)
dtdzβ=fzβ(t,x,y,z)
whose initial conditions t0β,x0β,y0β,z0β are known. Then, for a step size h=tn+1ββtnβ, the RK4 method is given by
xn+1β=xnβ+61β(k1βxβ+2k2βxβ+2k3βxβ+k4βxβ)
where [w1ββw2ββw3ββw4ββ]=[61ββ62ββ62ββ61ββ].
Single Iteration
In this code, we will implement the idea to generate [kxββkyββkzββ] simultaneously. This will eliminate need for multiple for loops. Consider a single iteration of RK4. We will create a function f=[fxββfyββfzββ]. Then,
k=h.f(v)
For k1β, this takes in a vector v=[tnββxnββynββznββ] and generates an output. This output, when multiplied with h, will provide k1β=[k1βxββk1βyββk1βzββ].
For k2β, this takes in the vector [tnβ+2hββxnβ+2k1βxβββynβ+2k1βyβββznβ+2k1βzβββ]. Output of this iteration multiplied with h to generate k2β=[k2βxββk2βyββk2βzββ].
Similarly, to generate k3β, we will use [tnβ+2hββxnβ+2k2βxβββynβ+2k2βyβββznβ+2k2βzβββ]
For k4β, we will use [tnβ+hβxnβ+k3βxββynβ+k3βyββznβ+k3βzββ]
Also note that, to implement each step we need the output of previous step. Hence, we cannot vectorise this process. Now, consider the step #2 from above. It can be rewritten as,
v2β=[tnββxnββynββznββ]+c[hβk1βxββk1βyββk1βzββ]
where, c=21β.
This can be generalised as,
$$
\displaystyle{v_j = v + c_j \begin{bmatrix} h & {k_{j-1}}x & {k{j-1}}y & {k{j-1}}_z \end{bmatrix}}
$$
where jβ{2,3,4}. This has been implemented in the code as
function[values, K] =multi_rk4(@fns, in_val, t, h)% The inputs are% fns(tn, xn, yn, zn ...) is a function which outputs [fx_n fy_n fz_n ...]% in_val = [x0 y0 z0 ... ] % Initial values.% No of values 'm = length(in_val)'% t = [t0 tn] %tn being the last value of t% h is the time step h = t_n+1 - t_n% The outputs are% values = [n x m] matrix% EXAMPLE% values = [x0 y0 z0;% x1 y1 z1;% .% .% xn yn zn]% K = [4 x m x n-1] matrix. Its [4 x m] matrix stacked for n-1 iterations.%-----weight vectors-----%weights=[1/6,2/6,2/6,1/6];c=[00.50.51];%-----eo vectors-----%%-----initialize-----%t=t(1):h:t(2);n=length(t);m=length(in_val);values=zeros(n,m);values(1,:)=in_val;K=zeros(4,m,n-1);%-----eo initialize-----%fori=1:n-1forj=1:4ifj==1inputs=[t(i)values(i,:)];elseinputs=[t(i)+c(j)*h(values(i,:)+c(j)*K(j-1,:,i))];endifK(j,:,i)=h*fns(inputs);endforvalues(i+1,:)=values(i,:)+weights*K(:,:,i);endforendfunction
I have written this code as a part of my learning experience. If there can be any improvements or clarifications, please let me know in the comments section below.
Please note that, in the above implementation, I have made k=hf which is different from the general implementations as seen in Wikipedia