help PLease, RK4!

clc
clear all
% constants
m1 = 60;
m2 = 70;
m3 = 80;
c1 = 5;
c3 = 5;
c2 = 10;
k1 = 50;
k3 = 50;
k2 = 100;
%given equations
fx = @(time,x1,x2,x3) m1*x1 +c1*x1+c2*(x1-x2)+k1*x1+k2*(x1-x2);
fy = @(time,x1,x2,x3) m2*x2+c3*(x2-x3)+c2*(x2-x1)+k3*(x2-x3)+k2(x2-x1);
fz = @(time,x1,x2,x3) m3*x3+c3*(x3-x2)+k3(x3-x2);
%initialize
x(1) = 1;
y(1) = 1;
z(1) = 1;
h = 1;
t(1) = 0;

for i=1:500 
t(i+1) = t(i)+h;
k1x = fx(t(i),x(i),y(i),z(i));
k1y = fy(t(i),x(i),y(i),z(i));
k1z = fz(t(i),x(i),y(i),z(i));

k2x = fx(t(i)+h/2, x(i)+h/2*k1x, y(i)+h/2*k1y, z(i)+h/2*k1z);
k2y = fy(t(i)+h/2, x(i)+h/2*k1x, y(i)+h/2*k1y, z(i)+h/2*k1z);
k2z = fz(t(i)+h/2, x(i)+h/2*k1x, y(i)+h/2*k1y, z(i)+h/2*k1z);

k3x = fx(t(i)+h/2, x(i)+h/2*k2x, y(i)+h/2*k2y, z(i)+h/2*k2z);
k3y = fy(t(i)+h/2, x(i)+h/2*k2x, y(i)+h/2*k2y, z(i)+h/2*k2z);
k3z = fz(t(i)+h/2, x(i)+h/2*k2x, y(i)+h/2*k2y, z(i)+h/2*k2z);

k4x = fx(t(i)+h, x(i)+h*k3x, y(i)+h*k3y, z(i)+h*k3z);
k4y = fy(t(i)+h, x(i)+h*k3x, y(i)+h*k3y, z(i)+h*k3z);
k4z = fz(t(i)+h, x(i)+h*k3x, y(i)+h*k3y, z(i)+h*k3z);

x(i+1)=x(i)+h/6*( k1x + 2*k2x + 2*k3x + k4x);
y(i+1)=y(i)+h/6*( k1y + 2*k2y + 2*k3y + k4y);
z(i+1)=z(i)+h/6*( k1z + 2*k2z + 2*k3z + k4z);
end

This is my code for RK4 for the given equations, however, I am completely stuck. It keeps on giving me an error that is as follows: "Array indices must be positive integers or logical values.

Error in FinalCode>@(time,x1,x2,x3)m2*x2+c3*(x2-x3)+c2*(x2-x1)+k3*(x2-x3)+k2(x2-x1)

Error in FinalCode (line 237)
k1y = fy(t(i),x(i),y(i),z(i));"

Please, I would really appreciate the help!!! I have tried everything and am going to lose it.

Thank you

Answer 1

[James Tursa]

You left out some * operators on the k2 and k3. This

fy = @(time,x1,x2,x3) m2*x2+c3*(x2-x3)+c2*(x2-x1)+k3*(x2-x3)+k2(x2-x1);
fz = @(time,x1,x2,x3) m3*x3+c3*(x3-x2)+k3(x3-x2);

should be this

fy = @(time,x1,x2,x3) m2*x2+c3*(x2-x3)+c2*(x2-x1)+k3*(x2-x3)+k2*(x2-x1);
fz = @(time,x1,x2,x3) m3*x3+c3*(x3-x2)+k3*(x3-x2);

[Layla Bitar]

oh my god. I can't believe I missed that this whole time -.- (tears). It has literally caused me so much trouble Thank you! It has now worked ,however for some reason, the graph did not converge. Would it be something wrong with my runge-kutta algorithm? I did try changing with my stepsize

[James Tursa]

What are the original differential equations you are trying to solve? And the original initial conditions? What physical system do they represent?

[Layla Bitar]

Here is the problem that I was given. The only initial condition that was given is t = 0

I would appreciate the help, Thank you!

[James Tursa]

Aha! You have a 6th order system, not a 3rd order system. The total order of a set of differential equations is obtained by adding up the highest order derivative for each variable. In your case, you have x1dotdot, x2dotdot, and x3dotdot as the highest order derivatives. So 2nd order for each of these, making the total order of your system 2 + 2 + 2 = 6. You are going to need six elements to represent this system for integration, not three. The six elements are all of the variables and derivatives up to but not including the highest order derivatives:

x1
x2
x3
x1dot
x2dot
x3dot

I would suggest you start using a 6-element state vector for this instead of using different variables as you are currently doing. You will end up writing 1/6th of the code this way. I would suggest using y for this state vector name to coincide with the nomenclature used in the MATLAB doc. E.g., define the y elements as

y(1) = x1
y(2) = x2
y(3) = x3
y(4) = x1dot
y(5) = x2dot
y(6) = x3dot

Then write a single derivative function for this state vector. Again, using the nomenclature from the doc (e.g., ode45) I would suggest it have the signature (t,y) ... the t is included so that it matches the MATLAB doc even though you don't need t for your equations. So you would write code for this:

function dydt = myderiv(t,y,m1,m2,m3,c1,c2,c3,k1,k2,k3)
dydt = zeros(size(y));
x1 = y(1);
x2 = y(2);
x3 = y(3);
x1dot = y(4);
x2dot = y(5);
x3dot = y(6);
x1dotdot = _____; % you fill this in
x2dotdot = _____; % you fill this in
x3dotdot = _____; % you fill this in
dydt = [ x1dot; x2dot; x3dot; x1dotdot; x2dotdot; x3dotdot];
end

The "fill in" parts above are obtained by solving your original differential equations for the highest order term in each equation. I will leave this up to you.

The last step is to rewrite your hand-written RK4 code to use a state vector instead of multiple scalar variables. This has the advantage of reducing your code and simplifying it. E.g., consider the state vector as a column vector, then the solution will consists of a series of column vectors. We can store all of these in a single matrix, each column of the matrix being a solution at a specific time. Your RK4 code would then reduce to this:

f = @(t,y) myderiv(t,y,m1,m2,m3,c1,c2,c3,k1,k2,k3);
for i=1:500 
t(i+1) = t(i) + h;
k1 = f( t(i)      , y(:,i)         );
k2 = f( t(i) + h/2, y(:,i) + h/2*k1);
k3 = f( t(i) + h/2, y(:,i) + h/2*k2);
k4 = f( t(i) + h  , y(:,i) + h  *k3);
y(:,i+1) = y(:,i) + h/6*( k1 + 2*k2 + 2*k3 + k4 );
end

You simply need to start it out by initializing the first time and state vector:

t(1) = ____; % you fill this in
y = zeros(6,501); % allocate the solution state matrix
y(:,1) = ____; % you fill this in, first column is the starting initial states.

See if you can implement this approach. Another advantage of this is that you can pass all of your time span, initial conditions, and derivative function to ode45( ) directly to check your solution.

[Layla Bitar]

WHAT A CLEAR EXPLANATION! Thank you so much. I have been struggling to understand what the dots signify. I will get cracking on this, now with your explanation, it has become much more clear as to how I should code for the RK4 algorithm.

Thank you again!!!!

[James Tursa]

The "dots" are derivatives with respect to time. I.e.,

x1dot = d(x1)/dt <-- First derivative of x1 with respect to time

x1dotdot = d2(x1)/dt^2 <-- Second derivative of x1 with respect to time

etc.

Your RK4 looping code was actually just fine ... it was your derivatives and state definitions that were messed up.

[Layla Bitar]

Yes makes a lot of sense!

update: I was able to get my graph to converge and it did level off. Thank you so much!!!! Don't know if I can thank you enough!