Spring Pendulum

HOME PAGE Web Site Contents Mars Report Contents Mars Report Abstract CV for Dr. David Roffman Diplomas PhD Thesis PhD Thesis Powerpoint Mars PowerPoint MSL Weather Reports MSL Yr 3 Winter Weather MSL Fall Yr 3 Weather MSL Yr. 3 Summer Weather MSL Yr. 3 Spring Weather Martian plume March 25 2017 MSL Ultraviolet Desai, EDL, Parachutes & ExoMars Mars winter vs. summer temps Sea at Utopia Planitia, Mars Tree Stump at MSL? Spherical life on Mars? Mars Report Abstract, 1-1.2 Mars Report Sec.2-2.1 Report 2.2-2.4 Report 2.5-2.5.2 Report 2.5.3-2.7 Report 3-4.1.2 Report 5 to 6 Report  7-7.2.1 Report 8 Report 9 Report 10-11 Report  12-12.2 Report 12.3-12.5 Report 12.6 Report 13-14 Report 14.1 Report 14.2-14.3 Report 14.4-14.6.2 Report 14.6.3-14.7 Report 15-19 Report References Report Afterword Rebuttal of REMS Report Running water on Mars MSL Year 0 Weather MSL Yr 2 Winter-Spring Weather MSL Yr 2 Summer Weather MSL Yr 2 Fall Weather MSL Yr 2-3 Winter Weather Adiabatics MSL Hi Temps MSL Low Temps Organic Chem found by MSL Oxygen in Mars Air MSL Day length & Temp Warm winter ground temps 155-Mile High Mars Plume Radiation Diurnal Air Temp Variation Mars Temps Fahrenheit Beagle found JPL/NASA Pressure Mistakes Enter MarsCorrect Sol 370, 1160 & 1161 Histories Mars-Radio-Show JPL Fudges Pressure Curves MSL Temp. ∆ Mast to Ground High & Low Pressures Normalized Mars soil 2% water Moving rock Mars MAVEN MSL Relative Humidity Claim Ashima Concedes Original MSL Weather Record Old MSL Weather Record MSL Summer Weather Pressure Estimate REMS Wind MSL Pressures REMS Reports Curiosity Geology CERN-2013-pics Daylight Math MSL Errors P1 MSL Errors P2 MSL-Chute-Flap MSL daylight Ashima Sols 15 to 111 Ashima Sol 112 to 226 Ashima Sol 227 on New Ashima Sols 270+ MSL Summer to Sol 316 Updated Secrets of Mars Weather Forecast Wind Booms MSL Credibility MSL Temp. Swings MSL Temperatures Sample Analysis at Mars (SAM) VL2 - MSL Ls Comparson Ashima MIT Mars GCM Dust Storm Nonsense Mars Slideshow Moving Sand & Martian Wind 3 DEC12 Press Conf. MSL Press Conf. 15NOV2012 Sol Numbering MSL Pressure Graph to Ls 218.8 MSL Sky Color Mars Sky Color DATA DEBATE! Zubrin's Letter Phoenix Vaisala Vaisala Pressure Sensors Phoenix &MSL Flawed MSL REMS Viking pressure sensors failed MSL landing site Mars Landings Phobos Grunt Martian Air Supersaturation Mars & CH4 Mars and MSL Time Viking Pressure Audit Links Mars Society 2008 Quant Finance Frontiers Home Front. Preface Frontiers Ch. 1 Frontiers Ch. 2 Antimatter Lightning Frontiers Ch. 3 Frontiers Ch. 4 Frontiers Ch. 5 Frontiers Ch. 6 Frontiers Ch. 7 Frontiers Ch. 8 Frontiers Ch. 9 Frontiers Ch 10 Frontiers Ch 11 Frontiers Ch 12 Frontiers Ch 13 Frontiers Ch 14 Frontiers Ch 15 Frontiers Ch 16 Frontiers Ch 17 Frontiers Ch 18 Frontiers Ch 19 Frontiers Ch 20 Frontiers Ch 21 Frontiers Ch 22 World Tour Spring-Break -13 Other Travels Asteroid Impact? ExoMars data Unit Issues Viking Pressures Tavis CADs Landing Long Scale Heights LS of Max/Min Pressures Tavis Report Tavis Failures Lander Altitude Martian Trees? Code Experiment Gedanken Report Mars Nuke? Martian Flares Mach Numbers MOLA (altitude) Original Mars Report Mariner 9 & Pressure Mars  Temps MSL Time MPF Pressure Blog Debates Spring Pendulum Plasma Model Reporting Errors Orbital Parameters Anderson Localization P. 1 Anderson Localization P. 2 Moving rock old Navigating Mars Mars Report Section Links Mars Report Figure Link Gillespie Lake rock outcrop MSL Sol 200 Anomaly Sol 1300&1301 Anomalies Gilbert Levin & Labeled Release Brine on Mars Ceres Lights Yr 1 Table 1 amfivan Missing data Mitchell Report Old Mars Report All MPF Temps ExoMars fails Did Spirit find past life? MSL ground temps go haywire Seasonal Pressure Altitude Calculations OPACITY AT MSL

Spring Pendulum (may be updated later)

I will start solving physical systems numerically with MATLAB.  This is the first of these problems. 

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')

 

Figure 1

Figure 2

With the initial conditions given in the code, the motion in somewhat periodic, but the pattern of radial position shifts up and down with time.