Anyway to simplify/generalise my RK4 code?

Hi All,

I'm new to using matlab and am struggling to simplify/generalise my code. My code is for 3 coupled equations using the Rk4 method. I need to make it so I can add in more equations and it won't be 100 lines long (if that makes sense).

The three equations I am using in this code are:

 y1' = -y2 - y3;
 y2' = y1 + (k1*y2);
 y3' = k2 + y3*(y1 - k3);
 with k1=k2= 0.2  ; k3=5.7

Here's my code:

 t0 = 0;
 tf = 10;
 N = 1000;
 h = (tf-t0)/N;
 t = [t0:h:tf];

 k1 = 0.2;
 k2 = 0.2;
 k3 = 5.7;
 y1(1) = 1;   
 y2(1) = 1;   
 y3(1) = 1;   

 S_tyxz = @(y_1,y_2,y_3) -y_2 - y_3;
 M_tyxz = @(y_1,y_2,y_3) y_1 + (k1 * y_2);
 E_tyxz = @(y_1,y_2,y_3) k2 + y_3 *(y_1 - k3);

 for i=1:N                          
    s1 = S_tyxz(y1(i),y2(i),y3(i));
    s2 = S_tyxz((y1(i)+0.5*h*s1),(y2(i)+0.5*h*s1),(y3(i)+0.5*h*s1));
    s3 = S_tyxz((y1(i)+0.5*h*s2),(y2(i)+0.5*h*s2),(y3(i)+0.5*h*s2));
    s4 = S_tyxz((y1(i)+s3*h),(y2(i)+s3*h),(y3(i)+s3*h));

    m1 = M_tyxz(y1(i),y2(i),y3(i));
    m2 = M_tyxz((y1(i)+0.5*h*m1),(y2(i)+0.5*h*m1),(y3(i)+0.5*h*m1));
    m3 = M_tyxz((y1(i)+0.5*h*m2),(y2(i)+0.5*h*m2),(y3(i)+0.5*h*m2));
    m4 = M_tyxz((y1(i)+m3*h),(y2(i)+m3*h),(y3(i)+m3*h));

    e1 = E_tyxz(y1(i),y2(i),y3(i));
    e2 = E_tyxz((y1(i)+0.5*h*e1),(y2(i)+0.5*h*e1),(y3(i)+0.5*h*e1));
    e3 = E_tyxz((y1(i)+0.5*h*e2),(y2(i)+0.5*h*e2),(y3(i)+0.5*h*e2));
    e4 = E_tyxz((y1(i)+e3*h),(y2(i)+e3*h),(y3(i)+e3*h));

    y1(i+1) = y1(i) + (1/6)*(s1+2*s2+2*s3+s4)*h;
    y2(i+1) = y2(i) + (1/6)*(m1+2*m2+2*m3+m4)*h;
    y3(i+1) = y3(i) + (1/6)*(e1+2*e2+2*e3+e4)*h;

 end
 plot3(y1,y2,y3,'r');

Any suggestions would be really appreciated,

Thank you

Answer 1

[James Tursa]

1) Define your state as a column vector, the number of elements being the number of equations that you have (since they are all 1st order equations). If there were higher order equations then we would need more states.

2) Create a function file that calculates the derivative of your state vector at a given time. Each element of this vector will correspond to one of your derivative equations. Use the same signature for this function that the MATLAB ode functions use:

% File my_derivative.m
function dy = my_derivative(t,y)
% put code here to calculate dy as a function of t and y
return
end

E.g., for these equations

 y1' = -y2 - y3;
 y2' = y1 + (k1*y2);
 y3' = k2 + y3*(y1 - k3);
 with k1=k2= 0.2  ; k3=5.7

you could have something like this:

% File my_derivative.m
function dy = my_derivative(t,y)
k1 = 0.2;
k2 = 0.2;
k3 = 5.7;
dy = zeros(size(y));
dy(1) = -y(2) - y(3);
dy(2) = y(1) + k1*y(2);
dy(3) = k2 + y(3)*(y(1) - k3);
return
end

Alternatively, if you wanted to pass in the k's as parameters, you could do this:

% File my_derivative.m
function dy = my_derivative(t,y,k1,k2,k3)
dy = zeros(size(y));
dy(1) = -y(2) - y(3);
dy(2) = y(1) + k1*y(2);
dy(3) = k2 + y(3)*(y(1) - k3);
return
end

3) Then in your loop you only have one set of RK4 equations that operates on your state vector, instead of multiple RK4 equations that each operate only on a scalar. Put the iterative results in columns of the result. E.g., a simple outline would be this:

 y(:,1) = your initial column state vector
 t(1) = the initial time
 for i=1:N-1
    % Rewrite your derivative equations using (:,i) for indexing instead of (i)
    % You will only have y(:,i) here, not y1(:,i) and y2(:,i) etc.
    % Be sure to use t(i) in your equations instead of just t (if there is a time dependency)
    % Note that you will have only one set of RK4 equations here
    y(:,i+1) = y(:,i) + (1/6)*(s1+2*s2+2*s3+s4)*h;  % store iterative results in columns of y
    t(i+1) = t(i) + h;
 end