How to write general function of 4th Order Runge-Kutta Method?

I am trying to develop a Matlab function for the 4th Order Runge-Kutta Method. It needs to be able to work with any function for given initial conditions, step size, etc. and then plot the results afterwards. Here is the code that I have so far. There are no error in the function itself. However, when I try to use it, I get a couple of error messages that I have also shown below the mfile. I am not sure what the error messages mean or how I would go about correcting them. Any help would be much appreciated!

FUNCTION:

function [tp,yp] = rk4(dydt,tspan,y0,h)

if nargin < 4, error('at least 4 input arguments required'), end
if any(diff(tspan)<=0), error('tspan not ascending order'), end

n = length(tspan);
ti = tspan(1); tf = tspan(n);
tp = ti:h:tf;
yp = zeros(ti,length(tspan));

for n = ti:tf
    tp = ti;
    k1 = dydt(tp(n),y0(n));
    k2 = dydt(tp(n) + (1/2)*h, y0(n) + (1/2)*k1*h);
    k3 = dydt(tp(n) + (1/2)*h, y0(n) + (1/2)*k2*h);
    k4 = dydt(tp(n) + h, y0(n) + k3*h);
    yp = y0(n) + ((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h);
    ti = ti + h;
    y0 = yp;
end
plot(tp,yp)

COMMAND WINDOW:

>> [tp,yp] = rk4(@(CA) ((1/5)*(20-CA))-(0.12*CA),15,0,1)
Attempted to access tp(15); index out of bounds because
numel(tp)=1.

Error in rk4 (line 13)
    k1 = dydt(tp(n),y0(n));

[Stephanie Valerio]

How would you input a given equation?

Answer 1

[Star Strider]

I didn’t run your code, but see if changing the initial ‘tp’ assignment to:

tv = ti:h:tf;

and then in your loop, changing:

tp = ti;

to:

tp = tv(n);

and changing all the subsequent references to ‘tp(n)’ to simply ‘tp’.

Also, you will need change this line (adding a subscript to ‘yp’) to save ‘yp’ at each iteration:

yp(n) = y0(n) + ((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h);

See if that improves things. You may have to make some other changes as well to accommodate these, and there may be other problems in your code, but this should keep you progressing towards a successful conclusion.

[Alicia]

I have still been getting the same error messages. Thanks for trying though!

[Star Strider]

My pleasure!

I assume ‘dydt’ is a function handle. What is it returning at each iteration?

How are you calling your ‘rk4’ function? What are its arguments?

[Alicia]

I am not sure what the dydt is returning, but for the rk4 I am putting in the function, the span of independent values, initial value for first independent, and then h is the step size.

[Star Strider]

I made some changes, and it works!

Updated code:

function [tp,yp] = rk4(dydt,tspan,y0,h)

if nargin < 4, error('at least 4 input arguments required'), end
if any(diff(tspan)<=0), error('tspan not ascending order'), end

n = diff(tspan)/h;
ti = tspan(1); tf = tspan(2);
tp = ti:h:tf;
yp = zeros(ti,length(tspan));

for n = 1:length(tp)
    k1 = dydt(tp(n),y0);
    k2 = dydt(tp(n) + (1/2)*h, y0 + (1/2)*k1*h);
    k3 = dydt(tp(n) + (1/2)*h, y0 + (1/2)*k2*h);
    k4 = dydt(tp(n) + h, y0 + k3*h);
    yp(n,:) = y0 + ((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h)
    ti = ti + h;
    y0 = yp(n);
end
plot(tp,yp)
end

[tp,yp] = rk4(@(CA,T) ((1/5)*(20-CA))-(0.12*CA), [0 15],0,1)

[Alicia]

Thank you so much!

[Star Strider]

My pleasure!

[Steve Chuks]

Hi all... @StarStrider & Alicia.

I have a similar work as to the Runge-Kutta method to solve for ODE.

Can you give me a hint on how the plot will look like and the error estimation of that? Just thinking of how to implement that on matlab.

I appreciate your assistance. Thanks.

Answer 2

[SHIVANI TIWARI]

hi everyone.. can u guys help me to generate the code for solving 1st order differential equations using improved runge kutta method.

[ahti Maître]

this worked for me at least, with y1 the derivative of your function,(yo,ti) the starting point, tf the final abscissa, and h the step size:

function [tp,yp] = rk4_1(y1,ti,tf,h,y0)

n=fix((tf-ti)/h)+1;

tv = ti:h:tf;

yp = zeros(1,n);

for j = 1:n

tp = tv(j);

k1 = y1(tp,y0);

k2 = y1(tp + (1/2)*h, y0 + (1/2)*k1*h);

k3 = y1(tp + (1/2)*h, y0 + (1/2)*k2*h);

k4 = y1(tp + h, y0 + k3*h);

yp(j) = y0+((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h);

ti = ti + h;

y0 = yp(j);

end

plot(tv,yp)

[ahti Maître]

largely inspired by alicias code