Runge-Kutta fourth order (rk4)

When implementing the Runge-Kutta order 4 scheme(rk4) u1,...,un are computed using the algorithm:

for i = 0,1,...,n-1

k1 = f (ti, ui),

k2 = f (ti + h/2, ui + h/2*k1),

k3 = f (ti + h/2, ui + h/2*k2),

k4 = f (ti + h, ui + h*k3),

ui+1 = ui + h/6 (k1 + 2k2 + 2k3 + k4).

Would part of this code be written in matlab as the following:

for i = 1:n

k1 = h * feval ( f, u(:,i) );

k2 = h * feval ( f, u(:,i) + k1/2 ); 

k3 = h * feval ( f, u(:,i) + k2/2 ); 

k4 = h * feval ( f, u(:,i) + k3 ); 

u(:,i+1) = u(:,i) + h*( k1 + 2*k2 + 2*k3 + k4 ) / 6; 

end;

I'm not too sure if the "h" in the last line (u(:,i+1)) should be there are not.

Any help would be very much appreciated.

Answer 1

[Jan]

Simply compare the pseudo code:

k4 = f (ti + h, ui + h*k3)

with your Matlab code:

k4 = h * feval(f, u(:,i) + k3)

to see, that the h is already inlcuded in your k1, k2, k3, k4.

A good idea is a test: Integrate a COS curve and look if the results are as expected.

[Rebecca]

Hi, thanks for your response Jan.

I've just found out that my initial coding was incorrect. So I have an updated version which is the following:

for i = 1:n

k1 = feval ( f, u(:,i) );

k2 = feval ( f, u(:,i) + h/2 * k1 );

k3 = feval ( f, u(:,i) + h/2 * k2 );

k4 = feval ( f, u(:,i) + h * k3 );

u(:,i+1) = u(:,i) + h*( k1 + 2*k2 + 2*k3 + k4 ) / 6;

end;

Does this seem about right?

This code produces an "egg-shaped" graph.

[Jan]

@Rebecca: Do you get an egg-shaped graph for integrating a COS function? Then you probably found successfully the SIN.

[Rebecca]

I am an absolute novice with matlab, since I have only been using it for a few weeks, so not too sure how to integrate a COS function in matlab. :S

[Matt Tearle]

What was your function f? (That produced the "egg-shaped graph")

But your code looks about right, from inspection. That said, function handles are nicer than feval... :)

Answer 2

[Meysam Mahooti]

https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78?s_tid=srchtitle