RK4/AB4, need help with correct code for 2 second order equations in Matlab

I have examples of code to follow for one second order equation but struggling to write correct code for 2 second order equations. grateful for any help!

[James Tursa]

One 2nd order equation will result in a 2-element state vector. Two 2nd order equations will result in a 4-element state vector. Just write your derivative function for a 4-element state vector. What have you tried so far?

[Anne Ellaway]

hi James, thanks for getting back to me so fast! I've started writing the code for 4 different arrays so far.

[James Tursa]

Let us know if you need more help. We can't give much advice until we see the differential equations and the code you have written.

[Anne Ellaway]

here's the code I have so far

function tryout()
close all 

h=0.01;
t0=0;
tf=10;

N=(tf-t0)/h;
t=(t0:h:tf);

%create array 
w=zeros(4,N+1); %x1
x=zeros(4,N+1); %x2
y=zeros(4,N+1); %v1
z=zeros(4,N+1); %v2


%initial conditions
w0(1)=1; %x1(0)=1
x0(2)=6; %x2(0)=6
y0(1)=0; %dx1/dt=v1=0
z0(2)=3; %dx2/dt=v2=3

%inputting functions 
f_x1='f_x1';
f_x2='f_x2';
f_V1='f_V2';
f_V2='f_V2';

%runge kutta 

for i=1:N
    k1 = h * feval(f_x1, w(i), x(i), y(i), z(i)); 
    k2 = h * feval(f_x1, w(i)+(k1/2), x(i)+(k1/2), y(i)+(k1/2), z(i)+(k1/2));
    k3 = h * feval(f_x1, w(i)+(k2/2), x(i)+(k2/2), y(i)+(k2/2), z(i)+(k2/2));
    k4 = h * feval(f_x1, w(i)+k3, x(i)+k3, y(i)+ k3, (z(i)+k3));
       
   w(i+1) = w(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;
   x(i+1) = x(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;
   y(i+1) = y(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;
   z(i+1) = z(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;

Answer 1

[James Tursa]

So, first define a 4-element state vector. To keep the nomenclature the same as the MATLAB docs, I will use the variable name y. The elements of y are defined as

y(1) = x1

y(2) = x2

y(3) = x1dot

y(4) = x2dot

Next step is to solve your matrix differential equations for x1dotdot and x2dotdot so that you have two expressions. You can do this easily on paper. I.e., do the matrix multiply on paper and then solve the two equations for the highest order derivatives. You will get:

x1dotdot = some expression involving x1, x2, x1dot, and x2dot

x2dotdot = some expression involving x1, x2, x1dot, and x2dot

Then rewrite these expressions using the definitions above, so that you will get

x1dotdot = some expression involving y(1), y(2), y(3), and y(4)

x2dotdot = some expression involving y(1), y(2), y(3), and y(4)

Now you are ready to write code. The derivative function will be

dy = @(t,y) [y(3); y(4); the expression for x1dotdot; the expression for x2dotdot]

This can be used in your RK4 code, or even passed into ode45( ).

[Anne Ellaway]

thank you so much James, I've pasted my code above hope I've done it right!

[James Tursa]

So, your latest code has too many state variables for my taste. You have w, x, y, z. While you can go this route, you will wind up writing 4 times the amount of code. We don't want that. We want simple code. So stick with the state variable approach. You only need this one line pre-allocation code:

y = zeros(4,N+1);

The definitions are:

y(1,i) is x1 at time t(i)

y(2,i) is x2 at time t(i)

y(3,i) is x1dot at time t(i)

y(4,i) is x2dot at time t(i)

I.e., y(:,i) is the full 4-element state vector at time t(i). Everywhere you want a state vector in your code, it will be of the form y(:,i)

Then your initialization code is simply

y0 = [1;6;0;3];
y(:,1) = y0; % The state vector at time t(1) is y0

The derivative function is what I have showed you (you fill in the expressions):

dy = @(t,y) [y(3); y(4); the expression for x1dotdot; the expression for x2dotdot];

And your RK4 loop would look something like this (use spaces to make the code more readable):

for i=1:N
    k1 = dy( t(i)      , y(:,i)          );
    k2 = dy( t(i) + h/2, y(:,i) + k1*h/2 );
    k3 = dy( t(i) + h/2, y(:,i) + k2*h/2 );
    k4 = dy( t(i) + h  , y(:,i) + k3*h   );
    y(:,i+1) = y(:,i) + h * ( k1 + 2*k2 + 2*k3 + k4 ) / 6;
end

I am just following the equations straight from this link (where my dy is f in the link):

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method

See if you can implement this approach, which is much simpler than carrying four separate state variables.

[Anne Ellaway]

Thank you so much James, I will try this tomorrow and see how it goes. You have been incredibly helpful and I really appreciate it.

[Anne Ellaway]

hi James, I'm still having a problem, as when I enter the dy function, the xdotdot equations have an unknown variable (v=dx/dt) but dx/dt can't be solved so I don't know how to state what v is so that MATLAB understands. many thanks again.

[James Tursa]

Please post your dy code. Also, I don't see the original DE pic you posted earlier. The matrix DE you had needs to be solved for the highest order derivatives x1dotdot and x2dotdot and that becomes part of your dy. You can either do this by hand and write the code explicitly, or use the backslash operator \ in MATLAB to do it for you. In general the variables are defined as

dx1/dt = x1dot = v1 = y(3,i)

dx2/dt = x2dot = v2 = y(4,i)

[Anne Ellaway]

dy = @(t,y) [V1; V2; -3.*x1+1.5.*x2-0.1.*V1+0.05.*V2; 1.5.*x1-3.*x2+0.05.*V1-0.1.*V2 ];

[James Tursa]

All of the calculations in the dy need to be in terms of t, y, or known constants. So replace the x1, x2, V1, V2, etc. with the definitions in terms of y:

dy = @(t,y) [y(3); y(4); -3.*y(1)+1.5.*y(2)-0.1.*y(3)+0.05.*y(4); 1.5.*y(1)-3.*y(2)+0.05.*y(3)-0.1.*y(4) ];

That being said, I would have left the equations in terms of m1, m2, k1, k2, k3, c1, c2, c3 for clarity and ease of changing things later on. E.g.,

m1 = whatever
m2 = whatever
k1 = whatever
k2 = whatever
k3 = whatever
c1 = whatever
c2 = whatever
c3 = whatever
M = [m1,0;
      0,m2];
K = [k1+k2,-k2;
    -k2   ,k2+k3];
C = [c1+c2,-c2;
    -c2   ,c2+c3];
dy = @(t,y) [y(3); y(4); M\(-K*y(1:2) - C*y(3:4))];

[Anne Ellaway]

Ah okay I understand what you mean. I will try to do this! Could you explain how the backslash operator works with this?

Thank you so much again!

[James Tursa]

The backslash operator was a bit of overkill here since the equations can be easily solved by hand. But basically you have this matrix equation

M*Xdotdot + K*X + C*Xdot = 0

where M, K, and C are as above and you have the following state vector definitions

X = [x1;x2]

Xdot = [x1dot;x2dot]

Xdotdot = [x1dotdot;x2dotdot]

You are trying to solve for the highest order derivatives in the equation, which is Xdotdot. So just doing the math:

M*Xdotdot = -K*X - C*Xdot

Now you have an equation in the form Ax=b. MATLAB solves this with the backslash operator, x = A\b. So in your case we get

Xdotdot = M \ (-K*X - C*Xdot)

Put in terms of your variable names, X = [x1;x2] = y(1:2), and Xdot = [x1dot;x2dot] = y(3:4).

[Anne Ellaway]

Thank you so much again James, I really appreciate all your help in this. All the best to you in these challenging times.