MATLAB: Runge-Kutta fourth order (rk4)

differential equationsinitial conditionsMATLABode

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.

Best Answer

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.

Related Solutions

MATLAB: Solving nonlinear scalar ode’s

Given that you have to use AB2 (with, apparently, a fixed stepsize) I going to take a guess that this is a numerical analysis homework problem? (If not, use one of the built-in MATLAB linear multistep methods, such as ode113.)

For fixed stepsize, you can calculate all the t values immediately. Look at either the linspace function or the range operator (:).

You can and should make room for the u values before you start to calculate them. Look at zeros.

Set the first value of u from the initial condition.

Loop to fill in the rest of the values, based on the previous ones. Given that you know how many steps you're going to take, use a for-loop.

You can hard-code the rate equation inside your loop if you like. But if you want to be a bit fancier, use an anonymous function handle to define it. Something like f = @(x,y) -2*y + ... Then you can evaluate it anywhere later by simply calling f(t(k),y(k)) (or whatever).

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

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.

Best Answer

Related Solutions

MATLAB: Solving nonlinear scalar ode’s

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

Related Question