June 17 2022
For the first part of this internship we compared the difference between Euler's method of integration with Euler-Cromer integration. It was observed that when using Euler's method in a Simple Harmonic Oscillator (SHO) the method would tend to increase the amount of energy in the system which is incorrect.
The energy issue is solved by using the Euler-Cromer method.
As shown on the figures, The Euler-Cromer method holds together much better and energy is conserved.
The code used is as follows.
% This script compares Euler's method of numerical integration to the analytical method.
% This script compares the Euler-cromer method of numerical integration to the analytical method.
for i = 2:n
June 24 2022
During this week we analyzed the Euler-Cromer method of numeric integration in a rigid pendulum. We also added resistance and oscillatory forces into our system.
As you can see the graphs are a lot more chaotic.
The code used is as follows:
g = 9.81;
omega(1) = 0;
July 1 2022
This week we studied the relationship between a cubic oscillator’s period and its amplitude. We calculated that the period of the oscillator is proportional to the inverse of the amplitude.
We also worked with animation and modeling techniques in Matlab.
This code calculates and animates a rigid pendulum.
%%% This code models and animates the motion of a rigid pendulum.
g = 9.81; %gravity
L = 2; % length
m = 5; % mass
b = 0; % resistance
I = (m * L^2) / 3; %moment of inirtia
Tau = 0; % added force
w = 2.5; % added force interval
omega(1) = 0; %initialize omega (initial velocity)
theta(1) = 30 * pi / 180; % initialize theta (initial angle)
alpha(1) = ((3 * g) / (2 * L))^2 * sin(theta(1)); %initialize alpha (initial acceleration)
nsteps = 1000; % sets number of steps
dt = .01; % sets time step
t(1) = 0; %initialize time
% Euler-Cromer method
for n = 2:nsteps
t(n) = t(n-1) + dt;
alpha(n) = -((3 * g) / (2 * L))^2 * sin(theta(n-1)) - (b*omega(n-1))/I + Tau * cos(w*t(n));
omega(n) = omega(n-1) + alpha(n) * dt;
theta(n) = theta(n-1) + omega(n) * dt;
X = -.5 * L .* sin(theta); % locates center of mass in the x direction
Y = -.5 * L .* cos(theta); % locates center of mass in the y direction
% animation for the pendulum
for z = 1:length(theta)
K(z) = getframe(gcf);
videofile = VideoWriter('pendulum.avi','Uncompressed AVI')
I am also working on a simulation of the ISS orbiting Earth, but I have run into an issue when calculating the angle and I have not been able to figure it out yet. Here is my code so far. This is still very rough.
mass_Earth = 5.972e24;
mass_object = 20000;
x(1) = 0;
y(1) = 6786100;
alt(1) = sqrt(x(1)^2 + y(1)^2);
vx(1) = 7778;
vy(1) = 0;
ax(1) = 0;
ay(1) = 0;
nsteps = 100;
dt = 60;
t(1) = 0;
theta(1) = 0;
for n = 2:nsteps
f(n) = -6.67428e-11 * (mass_Earth * mass_object) / alt(n-1)^2;
fx(n) = f(n) * sin(theta(n-1));
fy(n) = f(n) * cos(theta(n-1));
ax(n) = fx(n)/mass_object;
ay(n) = fy(n)/mass_object;
vx(n) = vx(n-1) + ax(n) * dt;
vy(n) = vy(n-1) + ay(n) * dt;
x(n) = x(n-1) + vx(n) * dt;
y(n) = y(n-1) + vy(n) * dt;
alt(n) = sqrt(x(n)^2 + y(n)^2);
v(n) = sqrt(vx(n)^2 + vy(n)^2);
theta(n) = theta(n-1) + (v(n)/(2 * pi * alt(n)))*dt;
t(n) = t(n-1)+dt;
July 15 2022
This week we worked on waves as a function of displacement and time. We were able to come up with a function to predict the future state of a wave.
future(n) = r2(current(n+1)+current(n-1))+2(1-r2)current(n)-past(n)
We applied this function to many different initial conditions including: fixed point “plucked” string, oscillating ends, and “sliding” end string.