Spring Pendulum
Spring Pendulum (may be updated later)
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:
dvx/dt = -ω2*(1-L/r)*x
dvy/dt = -ω2*(1-L/r)*y
dvz/dt = -ω2*(1-L/r)*z – g
dx/dt = vx
dy/dt = vy
dz/dy = vz
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')