Spring Pendulum

This program solves the spring pendulum problem.  RK4 is one of the more accurate methods for solving ODEs, but in this case the equations are such that it doesn’t improve the accuracy over other methods.  This is because there is no dependence in any of these equations of the variable inside the derivative itself.  An example is the x-velocity equation doesn’t have x-velocity on the right hand side.  Here are the six equations:

dv_x/dt = -ω²*(1-L/r)*x

dv_y/dt = -ω²*(1-L/r)*y

dv_z/dt = -ω²*(1-L/r)*z – g

dx/dt = v_x

dy/dt = v_y

dz/dy = v_z

Where “v” indicates velocity and the subscript is the direction, L is the un-stretched spring pendulum length, ω is the square root of the spring constant k over the mass of the pendulum bob m.  See Figure 1 for a rough sketch.

Below is the MATLAB code.  There are plots for the initial conditions given in the code below.

%Goal: Model the Spring Pendulum in 3-D with RK4

clc

clear

format 'long'

m = 1.5; %Mass in kg

g = 9.81; %Gravity in m/s^2

k = 2; %Spring constant in N/m

L = 1.1; %Unstretched length of spring in m

w = sqrt(k/m); %Angular frequency of ossilations in rad/s

x0 = 3; %Initial x-position in m

y0 = 1; %Initial y-position in m

z0 = -2.2; %Initial z-position in m

vx0 = .1; %Initial x-velocity in m/s

vy0 = .2; %Initial x-velocity in m/s

vz0 = .3; %Initial x-velocity in m/s

dt = .001; %Time step in s

t = 25; %Total time to run simulation in s

%Vectors store positions and velocities

xp = []; %Stores x-position

yp = []; %Stores y-position

zp = []; %Stores z-position

xv = []; %Stores x-velocity

yv = []; %Stores y-velocity

zv = []; %Stores z-velocity

rp = []; %Stores radial position

v = []; %Stores total velocity

for i = 1:t/dt

r = sqrt(x0^2 + y0^2 + z0^2); %Radial position

%Solves the x-velocity

dvx1 = dt*(-w^2*(1-L/r)*x0);

dvx2 = dvx1;

dvx3 = dvx2;

dvx4 = dvx3;

vx = vx0 + (dvx1 + 2*dvx2 + 2*dvx3 + dvx4)/6;

%Solves the y-velocity

dvy1 = dt*(-w^2*(1-L/r)*y0);

dvy2 = dvy1;

dvy3 = dvy2;

dvy4 = dvy3;

vy = vy0 + (dvy1 + 2*dvy2 + 2*dvy3 + dvy4)/6;

%Solves the z-velocity

dvz1 = dt*(-w^2*(1-L/r)*z0 - g);

dvz2 = dvz1;

dvz3 = dvz2;

dvz4 = dvz3;

vz = vz0 + (dvz1 + 2*dvz2 + 2*dvz3 + dvz4)/6;

%Solve the x-position

dx1 = dt*vx0;

dx2 = dx1;

dx3 = dx2;

dx4 = dx3;

x = x0 + (dx1 + 2*dx2 + 2*dx3 + dx4)/6;

%Solve the y-position

dy1 = dt*vy0;

dy2 = dy1;

dy3 = dy2;

dy4 = dy3;

y = y0 + (dy1 + 2*dy2 + 2*dy3 + dy4)/6;

%Solve the z-position

dz1 = dt*vz0;

dz2 = dz1;

dz3 = dz2;

dz4 = dz3;

z = z0 + (dz1 + 2*dz2 + 2*dz3 + dz4)/6;

%Appends vectors to save all positions and velocities

xp = [xp,x0];

yp = [yp,y0];

zp = [zp,z0];

xv = [xv,vx0];

yv = [yv,vy0];

zv = [zv,vz0];

rp = [rp,r];

v = [v,sqrt(vx0^2 + vy0^2 + vz0^2)];

%Set the current positions and velocities as the original to iterate.

x0 = x;

y0 = y;

z0 = z;

vx0 = vx;

vy0 = vy;

vz0 = vz;

end

%Plot Time vs. Radial Position

plot(1:t/dt,rp)

title('Number of Time-Steps vs. Radial Position')