Physics 303 Mechanics

%Fourier Series Square Wave

t = [0:0.01:2*pi];

b1=4/(1*pi);
b3=4/(3*pi);
b5=4/(5*pi);
b7=4/(7*pi);

F=zeros(1,length(t));

for i=1:length(t)-1
if(i<length(t)/2)
F(i)=+1;
else
F(i)=-1;
end
end

plot(t,F,t,b1*sin(t),"g",t,b3*sin(3*t),"r",t,b5*sin(5*t),"y",t,b7*sin(7*t),"k",t,b1*sin(t)+b3*sin(3*t)+b5*sin(5*t)+b7*sin(7*t));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%RK4_Newton

% Calculates the motion of a damped, driven oscillator using RK4 and

% inline function definition of the acting Force.

clc; clf;

clear all;

dt=0.1;

t = 0:dt:100;

a = zeros(1,length(t));

v = zeros(1,length(t));

x = zeros(1,length(t));

m = 1;k=1;c=1;

w0=sqrt(k/m);

wD=sqrt( (k/m)-(c/(2*m))^2);

F0=1;

wF=2;

x(1) = 1;

v(1) = 0;

F = @(t,x,v) -k*x - c*v + F0*cos(wF*t); % change the inline function as you desire

for i=1:(length(t)-1)

kv1 = F(t(i),x(i),v(i))/m;

kv2 = F(t(i),x(i)+0.5*dt,v(i)+0.5*dt*kv1)/m;

kv3 = F(t(i),x(i)+0.5*dt,v(i)+0.5*dt*kv2)/m;

kv4 = F(t(i),x(i)+dt,v(i)+dt*kv3)/m;

v(i+1) = v(i) + dt*(1/6)*(kv1+2*kv2+2*kv3+kv4); % main equation

x(i+1) = x(i) + dt*(v(i)+v(i+1))/2;

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

end

subplot(2,1,1);

plot(t,x,'r');

xlabel('t');ylabel('x');

title([char(F),' ,m=',num2str(m),' ,k=',num2str(k),' ,c=',num2str(c) ,' ,w0=',num2str(w0) ,' ,wD=',num2str(wD),' ,wF=',num2str(wF)]);

subplot(2,1,2);

plot(t,v,'b',t,F0*cos(wF*t),'g');

xlabel('t');ylabel('v in blue, Fext in green');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Simple harmonic motion - comparison of Euler, Euler Cromer,

% 2nd order Runge Kutta and built in MATLAB Runge Kutta

% function ODE45 to solve ODEs.

% Equation is d2x/dt2 = -w0^2 x

% Calculate the numerical solution using Euler method in red

% Needs Inputs (x0,v0,w0,N,dt)

[time,x] = SHM_Euler (10,0,1,1000,0.1);

subplot(2,2,1);

plot(time,x,'r' );

axis([0 100 -100 100]);

xlabel('Time');

ylabel('Displacement');

legend ('Euler method');

% Calculate the numerical solution using Euler Cromer method in blue

[time,x] = SHM_Euler_Cromer (10,0,1,1000,0.1);

subplot(2,2,2);

plot(time,x,'b' );

axis([0 100 -20 20]);

xlabel('Time');

ylabel('Displacement');

legend ('Euler Cromer method');

% Calculate the numerical solution using Second order Runge-Kutta method in green

[time,x] = SHM_Runge_Kutta2 (10,0,1,1000,0.1);

subplot(2,2,3);

plot(time,x,'g' );

axis([0 100 -20 20]);

xlabel('Time');

ylabel('Displacement');

legend ('RK2 method');

% Use the built in MATLAB ODE45 solver to solve the ODE

% The time values range from 0 to 100

% The initial velocity is 0, the initial displacement is 10

[t,x]=ode45(@SHM_ODE45_function,[0,100],[0,10]);

% We need to plot the second column of the y matrix, containing the

% displacement against time in black

subplot(2,2,4);

plot(t,x(:,2),'k');

axis([0 100 -20 20]);

xlabel('Time');

ylabel('Displacement');

legend ('ODE45 Solver');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [time,x] = SHM_Euler (x0,v0,w0,N,dt)

%%Needs Inputs x0,v0,w0,N,dt

a = zeros(N,1);

v = zeros(N,1); % initializes v, a col vector of dimension N X 1,to being all zeros

x = zeros(N,1); % initializes x, a col vector of dimension N X 1,to being all zeros

time = zeros(N,1); % this initializes the col vector time to being all zeros

x(1) = x0; % initial displacement

v(1) = v0; % initial velocity

% Euler solution

for step = 1:N-1 % loop over the timesteps

a(step) = -w0^2*x(step);

v(step+1) = v(step) + a(step)*dt;

x(step+1) = x(step) + v(step)*dt;

time(step+1) = time(step) + dt;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [time,x] = SHM_Euler_Cromer (x0,v0,w0,N,dt)

a = zeros(N,1);

v = zeros(N,1);

x = zeros(N,1);

time = zeros(N,1);

x(1) = x0;

v(1) = v0;

% Euler Cromer solution

for step = 1:N-1 % loop over the timesteps

a(step) = - w0^2 * x(step);

v(step+1) = v(step) + a(step)*dt;

x(step+1) = x(step) + v(step+1)*dt; %Notice v(step+1) insterad of v(step) in Euler

time(step+1) = time(step) + dt;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [time,x] = SHM_Runge_Kutta2 (x0,v0,w0,N,dt)

% 2nd order Runge Kutta solution

a = zeros(N,1);

v = zeros(N,1);

x = zeros(N,1);

time = zeros(N,1);

a_dash = zeros(N,1);

v_dash = zeros(N,1);

x_dash = zeros(N,1);

x(1) = x0;

v(1) = v0;

for step = 1:N-1 % loop over the timesteps

a(step)=- w0^2 * x(step);

v_dash(step)=v(step)+0.5*a(step)*dt;

x_dash(step)=x(step)+0.5*v(step)*dt;

a_dash(step) = - w0^2 * x_dash(step);

v(step+1) = v(step)+ a_dash(step)*dt;

x(step+1) = x(step)+ v_dash(step)*dt;

time(step+1) = time(step)+dt;

end;

%%%%%%%%%%%%%%%%%%%%%%%%

function dy = SHM_ODE45_function(~,y)

% y is the state vector

y1 = y(1); % y1 is displacement

v = y(2); % y2 is velocity

% write down the state equations

dy1=v;

dy2=-y1;

% collect the equations into a column vector, velocity in column 1,

% displacement in column 2 Kevin Berwick Page 15

dy = [dy1;dy2];

Problem

No Transient Question

Brachistochrone

Catenary

Experimental Catenary Program

-