I'm trying to solve a spring-mass-damper system using 4th order Runge-Kutta method. I came up with the following function, but I'm not able to get the outputs as needed.

function yout = smd_rk()

%% steps to define const, DEs

Fo=1;

m=1;

k=10;

b=100;

w=1;

f1=@(t,y1,y2)(y2);

f2=@(t,y1,y2)(Fo*sin(w*t)-b*y2-k*y1)/m;

%% Step size and initial conditions

h=0.01;

t(1)=0;

y1(1)=0;

y2(1)=0;

%% For loop for each timestep

for i=1:1000

t(i+1)=t(i)+h;

k1y1 = h*f1(t(i),y1(i),y2(i));

k1y2 = h*f2(t(i),y1(i),y2(i));

k2y1 = h*f1(t(i)+h/2,y1(i)+k1y1/2,y2(i)+k1y2/2);

k2y2 = h*f2(t(i)+h/2,y1(i)+k1y1/2,y2(i)+k1y2/2);

k3y1 = h*f1(t(i)+h/2,y1(i)+k2y1/2,y2(i)+k2y2/2);

k3y2 = h*f2(t(i)+h/2,y1(i)+k2y1/2,y2(i)+k2y2/2);

k4y1 = h*f1(t(i)+h,y1(i)+k3y1,y2(i)+k3y2);

k4y2 = h*f2(t(i)+h,y1(i)+k3y1,y2(i)+k3y2);

y1(i+1)=y1(i) + (k1y1 + 2*k2y1 + 2*k3y1 + k4y1)/6;

y2(i+1)=y2(i) + (k1y2 + 2*k2y2 + 2*k3y2 + k4y2)/6;

end

%% Output

yout = [t y1 y2];

%% plot command

plot(t,y1)

%% I'm trying to get the time, displacement and velocity as a table but am unable to. Also I feel like the equation may be going wrong somewhere after looking at the plot, but I can't figure out where.

Answer 1

[Ayush Modi]

Hi Ishit,

As per my understanding, you would like to store the calculated values “time”, “displacement”, “velocity” into a table. And you would like to change the equation as the “time-displacement” plot does not conform to the behaviour of the spring-mass-damper system.

Here is an example showing how you can implement the spring-mass-damper system.

    % Calculate the derivatives
    dxdt = v(i);  
    dvdt = -(k/m) * x(i) - (c/m) * v(i);
    % Update the state variables using the fourth-order Runge-Kutta method
    k1x = dxdt;
    k1v = dvdt;
    k2x = v(i) + 0.5 * dt * k1v;
    k2v = -(k/m) * (x(i) + 0.5 * dt * k1x) - (c/m) * k2x;
    k3x = v(i) + 0.5 * dt * k2v;
    k3v = -(k/m) * (x(i) + 0.5 * dt * k2x) - (c/m) * k3x;
    k4x = v(i) + dt * k3v;
    k4v = -(k/m) * (x(i) + dt * k3x) - (c/m) * k4x;
    % Update the state variables
    t(i+1) = t(i) + dt;
    x(i+1) = x(i) + (dt/6) * (k1x + 2*k2x + 2*k3x + k4x);
    v(i+1) = v(i) + (dt/6) * (k1v + 2*k2v + 2*k3v + k4v);

You can store the calculated values by using the “array2table” function. Refer to below code for the same:

    arrr = [t x v]  
    output_table = array2table(arrr);

Please refer to the following documentation to know more about “array2table” function:

https://www.mathworks.com/help/matlab/ref/array2table.html

I hope this resolves the issue you were facing.