help PLease, RK4!
2 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
James Tursa
on 5 May 2020
Edited: James Tursa
on 5 May 2020
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);
7 Comments
James Tursa
on 6 May 2020
Edited: James Tursa
on 6 May 2020
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!