project updates and summaries
WEEK ONE
06/11/22  06/17/22
UPDATES: Following a meeting on Monday, I spent the following day learning the basics of Octave/MATLAB. This included downloading Octave onto my computer and understadning arithmetic operations, variables, matricies, loops, scripts, functions, and plotting. While being introduced to the basics of octave, I managed to plot arbitrary values onto a graph which would turn out to be helpful experience for the upcoming projects.
EXERCISE #1  PLOTTING EULER'S METHOD
Exercise #1 included using Euler's Method to build a computational model of a simple hanging harmonic oscillator and comparing this with the analytical method of solving for the exact function of position as a function of time. As seen in image #1 below, the red graph (analytical) and the blue graph (numeric) do not align, thus hinting towards potentially inaccurate in terms of Euler's Method when dealing with energy. This is further corroborated when plotting the energy as a function of time (image #2), which illustrates the numerical energy going astray.
image #1
image #2
Image #1 Code 
Image #2 Code 
% position as a function of time  euler's method
% initial parameters
dt = 0.01; %interval between steps
v1 = 0; %initial velocity
m = 1.0; %object mass
y1 = .1; %initial position
g = 9.8; %gravity
k = 1; %spring constant
tsteps = 1000;
% define and start arrays for numerical solutions
time=zeros(tsteps);
y_Euler = zeros(tsteps);
v_Euler = zeros(tsteps);
% analytical solution
y_exact = zeros(tsteps);
v_exact = zeros(tsteps);
% initial conditions
time(1) = 0;
y_Euler(1)=y1;
y_exact(1)=y1;
v_Euler(1)=v1;
v_exact(1)=v1;
% main loop
for i=2:tsteps
time(i) = (i1)*dt;
v_Euler(i)=v_Euler(i1)k*y_Euler(i1)*dt/m;
y_Euler(i)=y_Euler(i1)+v_Euler(i1)*dt;
v_exact(i)=y1*sqrt(k/m)*sin(sqrt(k/m)*time(i));
y_exact(i)=y1*cos(sqrt(k/m)*time(i));
end
% plotting results
plot(time,y_Euler,'b')
hold on
plot(time,y_exact,'r')
xlabel('time (s)');
ylabel('position (m)');

% including energy comparison
% Euler's Method  Numerical and Analytical Comparison
% parameters
k = 1; % spring constant [N/m].
g = 9.81; % acceleration of gravity [m/s^2].
m = 1; % mass of the object [kg].
y1 = 5; % initial position[m].
v1 = 0; % initial velocity [m/s].
dt = 0.01; % time interval between steps [s].
tstep = 10000; % number of steps.
time(1) = 0; % initial time condition.
y(1) = y1; % initial position condition.
v(1) = v1; % initial velocity condition.
p_energy(1) = .5 * k * y1^2;
k_energy(1) = .5 * m * v1^2;
energy(1) = p_energy(1) + k_energy(1);
% numerical method
for i = 2:tstep
time(i) = time(i1) + dt; % adds time interval to previous time
v(i) = v(i1)  (k./m) .* y(i1) .* dt; % calculates new velocity
y(i) = y(i1) + v(i1) * dt; % calculates new position
k_energy(i) = .5 * m * v(i)^2;
p_energy(i) = .5 * k * y(i)^2;
energy(i) = p_energy(i) + k_energy(i);
end
% analytical method
y_exact = y1 * cos(sqrt(k ./ m) .* time);
v_exact = y1 * sqrt(k ./ m) .* sin(sqrt(k ./ m) .* time);
k_energy_exact = .5 .* m .* v_exact .^ 2;
p_energy_exact = .5 .* k .* y_exact .^ 2;
e_exact = p_energy_exact + k_energy_exact;
% plot results of numerical and analytical computations
plot(time,energy,'c',time,e_exact,'r')
hold on
title('Comparison of Numerical and Analytical  Eulers Method')
xlabel('time');
ylabel('position');
%ploting position
plot(time,y,'b')
hold on
plot(time,y_exact,'r')
legend('Euler Numerical Energy', 'Euler Exact Energy','Numerical Position','Exact Position');

WEEK TWO
06/18/22  06/24/22
UPDATES: This week, we focused our efforts on the Euler Cromer method, and then comparisons between small angle approximations and the eulercromer method. Not only did I gain more insight onto the systems, I also the learned crucial lesson of being mindful of syntax on MATLAB/Octave.
EXERCISE #2  SOLVING WITH EULERCROMER METHOD
With the EulerCromer Method, the final velocity needed to be utilized, instead of the initial velocity. Thus, to account for the change, the velocity formula needed to be altered to "v(i) = v1 + acceleration * time" and the code required the inclusion of acceleration and force. Ultimately, the EulerCromer Method proved to be successful in solving the energy differences between the analytical and numeric solving methods.
CODE 
% euler cromer method comparison
% parameters
k = 1; %spring constant [N/m].
g = 9.81; %acceleration of gravity [m/s^2].
m = 1; %mass of the object [kg].
y1 = 5; %initial position[m].
v1 = 0; %initial velocity [m/s].
dt = 0.01; %time interval between steps [s].
tstep = 10000; %number of steps.
time(1) = 0; %initial time condition.
y(1) = y1; %initial position condition.
v(1) = v1; %initial velocity condition.
force(1) = k*y1; %initial force (K = F/y > Ky = F > K*y)
p_energy(1) = .5 * k * y1^2;
k_energy(1) = .5 * m * v1^2;
energy(1) = p_energy(1) + k_energy(1);
% numerical method
for i = 2:tstep
time(i) = time(i1) + dt; % adds time interval to previous time
force(i) = k.*y(i1);
a(i) = force(i)/m;
v(i) = v(i1) + a(i) * dt; % final velocity = v1 +acceleration(time)
y(i) = y(i1) + v(i) * dt; % calculates new position
p_energy(i) = .5 * k * y(i)^2;
k_energy(i) = .5 * m * v(i)^2;
energy(i) = p_energy(i) + k_energy(i);
end
% analytical method
v_exact = y1 * sqrt(k./ m) .* sin(sqrt(k ./ m) .* time);
y_exact = y1 * cos(sqrt(k./ m) .* time);
p_energy_exact = .5 .* k .* y_exact .^ 2;
k_energy_exact = .5 .* m .* v_exact .^ 2;
e_exact = p_energy_exact + k_energy_exact;
% plot results of numerical and analytical computations
plot(time,energy,'c',time,e_exact,'g')
hold on
title('Euler vs Exact')
xlabel('time');
ylabel('position');
%ploting position
plot(time,y,'b')
hold on
plot(time,y_exact,'r')
legend('Euler Numerical Energy', 'Euler Exact Energy','Numerical Position','Exact Position');

EXERCISE #3  RIGID PENDULUM
The following pendulum system included the comparisons between small angle approximations and a numerical eulercromer usage. Although following similar base structures as the prior exericse, this exercise included new variables and components such as the new value for the coefficient of friction (b) and the inclution of torque (tau * cos(w*time)). The primary goal of the comparison between the two types of solution methods (analytical vs numerical) was to determine which one had more accuracy and/or whether one stopped working after a certain limit.
Between a certain limit for theta, the graphs (of theta and omega for both forms) were pretty expected, yet as the size of theta increased, the small angle approximation turned to be less accurate (and more chaotic) as seen in the non standard looking spiral. The graphs below (image #4, image #5) exist when most initial conditions are 0 and the theta value is 1.
However, if one was to change some of the parameters, such as increasing the theta value and changing most of the initial parameters from 0 to another value, the results would turn more chaotic. To further illustrate that the small angle approximation is not adaquet for all values of theta, image #7 illustrates the loss of energy and general chaotic structure of the spiral.
image #4 image #5
The code for both results is seen in the first column. The second column depicts the slight changes to the code that create the images/results following the code. This second code is virtually the exact same as the prior, but without the specific legends, the 0 values for certain parameters, and the slightly small value for theta.
CODE FOR IMAGE #4 AND #5
% euler cromer method  small angle approximation and numerical
% parameters
k = 1; %spring constant [N/m].
g = 9.8; %acceleration of gravity [m/s^2].
m = 1; %mass of the object [kg].
dt = 0.01; %time interval between steps [s].
theta(1) = 1;
om(1)= 0; % omega
tstep = 1000; % number of steps.
time(1) = 0; % initial time condition.
L = 1; % length
I = (m*L^2)/3
Tau = 0; % value between 0 10
w = 0; % between 0 and 3 rad
b = 0; % value between 0 and 1
a(1) = ((m*g*L)/2)/I * sin(theta(1))(b*om(1))/I + Tau*cos(w*time(1)); % initializing acceleration
for i = 2:tstep
time(i) = time(i1) + dt; % time interval
om(i) = om(i1)+dt*a(i1); % omega interval
theta(i) = theta(i1)+dt *om(i); % theta interval
a(i) = ((m*g*L)/2)/I * sin(theta(i))(b*om(i))/I + Tau*cos(w*time(i)); % acceleration
t(i) = time(i1) +dt;
end
theta_exact = theta(1) .* cos(((3 * g)/(2 * L)) .* time);
om_exact = theta(1).*(3*g)/(2*L).*(sin(((3*g)/(2*L)).*time));
title('Small Angle Approx VS EulerCromer')
xlabel('time');
ylabel('position');
plot (time, theta,'r', time, theta_exact, 'b')
legend ('Small Angle Approx', 'EulerCromer')
figure
plot (theta, om,'g')
hold on
plot (theta_exact, om_exact, 'r')
legend ('Small Angle Approx', 'EulerCromer')

SLIGHT ALTERATIONS
% euler cromer method comparison
% parameters
k = 1; %spring constant [N/m].
g = 9.8; %acceleration of gravity [m/s^2].
m = 1; %mass of the object [kg].
dt = 0.01; %time interval between steps [s].
theta(1) = (120*pi)/180;
om(1)= 0;
tstep = 1000; %number of steps.
time(1) = 0; %initial time condition.
L = 1;
I = (m*L^2)/3
Tau = 2; %value between 0 10
w = 3; % between 0 and 3 rad
b = 0.2; % value between 0 and 1
a(1) = ((m*g*L)/2)/I * sin(theta(1))(b*om(1))/I + Tau*cos(w*time(1));
for i = 2:tstep
time(i) = time(i1) + dt; % adds time interval to previous time
om(i) = om(i1)+dt*a(i1);
theta(i) = theta(i1)+dt *om(i);
a(i) = ((m*g*L)/2)/I * sin(theta(i))(b*om(i))/I + Tau*cos(w*time(i));
t(i) = time(i1) +dt;
end
theta_exact = theta(1) .* cos(((3 * g)/(2 * L)) .* time);
om_exact = theta(1).*(3*g)/(2*L).*(sin(((3*g)/(2*L)).*time));
title('Euler vs Exact')
xlabel('time');
ylabel('position');
plot (time, theta,'r', time, theta_exact, 'b')
figure
plot (theta, om)
hold on
plot (theta_exact, om_exact)

image # 6 image #7 image #8
WEEK THREE
06/24/22  07/01/22
UPDATES: Over the weekend, I was tasked with finishing the exercises for the ridgid pendulum (using the code above). Although we had finished most of it before this weekend, exercise four asked if there was a way to dampen the pendulum to the point it doesn't oscilate. Following this, we focused on a cubic oscilator's period and amplitude, and worked on numeric integration  all which were new concepts relatively new to me.
EXERCISE #4  RIGID PENDULUM CONT'D
For the subexercise, it was found that oscilations will decrease the closer b ([0,1]) gets to 1. Likewise, when experimenting with changing the tau value from 2 to 10, it was visible that as the tau value increased, the figures produced would grow more chaotic. At about 6 or 7 (for tau), the figures turned away from a coiled spiral to a more chaotic result. The code utilized was the same as in the exercise above.
Images #9 and #10 illustrate the results for tau = 2. Images #11 and #12 illustrate the results for tau = 10.
image #9 image #10
image #11 image #12
EXERCISE #5  CUBIC OSCILATOR
While the general parameters remained the same as in prior code, to plot the relationship between the period and amplitude of a cubic oscilator, we needed to include both of those factors. The result illustrated in image #13 showed a general trend between the period and amplitude. We further calculated that the period of the oscillator is proprotional to the inverse of the amplitude through numerical integration, which resulted to about 1.31 when using various different integration methods (including Simpson's 1/3 Rule).
%%%% cubic oscilator  period vs amplitude
pkg load signal
k = 1; %spring constant [N/m].
g = 9.8; %acceleration of gravity [m/s^2].
m = 1; %mass of the object [kg]
dt = 0.1; %time interval between steps [s].
tstep = 1000; % number of steps.
L = 1; % length
% loop for changing amplitude
for n = 1:17
theta(1) = (5*n)*(pi/180);
a(1) = theta(1)^3;
om(1)= 0; % omega
time(1) = 0; % initial time condition.
in_am(n) = 1/theta(1);
for i = 2:tstep
time(i) = time(i1) + dt; % time interval
om(i) = om(i1)+dt*a(i1); % omega interval
theta(i) = theta(i1)+dt *om(i); % theta interval
a(i) = theta(i)^3;
t(i) = time(i1) +dt;
end
[PKS,LOC] = findpeaks(abs(theta'));
LOC;
newtime = time(LOC);
period(n) = 2* (mean(diff(newtime)));
end
plot(in_am, period, 'x')
legend ('EulerCromer')
title('EulerCromer')

numeric integration
image # 13
EXERCISE #6  MODELLING A RIGID PENDULUM THROUGH MATLAB
Slightly altering the code in prior exercises for the rigid pendulum to account for center of mass in x and y direction, we were able to animate and model the pendulum through MATLAB. Linked below is a clip of the animation.
RIGID PENDULUM ANIMATION  YOUTUBE
% code models and animates motion of rigid pendulum
clear all
close all
clc
% parameters
k = 1; %spring constant [N/m].
g = 9.8; %acceleration of gravity [m/s^2].
m = 1; %mass of the object [kg].
dt = 0.01; %time interval between steps [s].
theta(1) = 1;
om(1)= 0; % omega
tstep = 100; % number of steps.
time(1) = 0; % initial time condition.
L = 2; % length
I = (m*L^2)/3
Tau = 0; % value between 0 10
w = 3; % between 0 and 3 rad
b = 0; % value between 0 and 1
a(1) = ((m*g*L)/2)/I * sin(theta(1))(b*om(1))/I + Tau*cos(w*time(1)); % initializing acceleration
for i = 2:tstep
time(i) = time(i1) + dt; % time interval
om(i) = om(i1)+dt*a(i1); % omega interval
theta(i) = theta(i1)+dt *om(i); % theta interval
a(i) = ((m*g*L)/2)/I * sin(theta(i))(b*om(i))/I + Tau*cos(w*time(i)); % acceleration
t(i) = time(i1) +dt;
end
theta_exact = theta(1) .* cos(((3 * g)/(2 * L)) .* time);
om_exact = theta(1).*(3*g)/(2*L).*(sin(((3*g)/(2*L)).*time));
% the following commands are commented out as they were unnecessary for this specific exercise but could be useful for other comparisons that were done in prior exercises and thus are not entirely deleted from the code
%title('Small Angle Approx VS EulerCromer')
%xlabel('time');
%ylabel('position');
%plot (time, theta,'r', time, theta_exact, 'b')
%legend ('Small Angle Approx', 'EulerCromer')
%figure
%plot (theta, om,'g')
%hold on
%plot (theta_exact, om_exact, 'r')
%legend ('Small Angle Approx', 'EulerCromer')
% pendulum animation values
X = .5*L.*sin(theta); % center of mass x direction
Y = .5*L.*cos(theta); % center of mass y direction
figure
% animation for pendulum loop
for n = 1:length(theta)
plot ([0,2*X(n)], [0,2*Y(n)],'linewidth', 2)
axis([2,2,2,0])
hold on
plot(X(n),Y(n),'rx') % center point
pause (.01)
% storing frames
K(n) = getframe(gcf);
hold off
end
videofile = VideoWriter('pendulum.avi', 'Uncompressed AVI')
open(videofile)
writeVideo(videofile,K)
close(videofile)

WEEK FOUR
07/01/22  07/08/22
UPDATE: I had managed to numerically integrate the prior integral in exercise #5 through octave (as prior to today, I had done simpson's 1/3 rule by hand which was unreliable and difficult). I also spent time experimenting with MATLAB and animating both a simple square and later a 3D cuplike shape, of which the codes are seen below. Because this was a shorter week, we had also started on a new exercise, which I have not completed just yet.
EXERCISE #7  EXPERIMENTING WITH ANIMATIONS/VIDEOS
ANIMATED SQUARE 
3D SHAPE 
% rotating square animation
clear all, close all, clc
t=0:pi/2:2*pi;
x=cos(t);
y=sin(t);
axis equal;
hold on;
for i=1:100
clf
hold on;
q=[x;y];
e=pi/10*i;
f=[cos(e) sin(e);sin(e) cos(e)];
l=f*q;
r=l(1,:);
d=l(2,:);
plot(r,d,'g','linewidth',3);
axis square;
pause(0.1);
end 
% 3d shape
h = animatedline('Color','#F5F5DC', 'LineWidth', 2.5);
t = 0:0.01:50;
x = exp(.005*t).*sin(10*t);
y = exp(.005*t).*cos(10*t);
z = t;
axis ([ 1 1 1 1 0 30]);
axis square
grid on
for k = 1:length(t)
addpoints (h,x(k), y(k), z(k));
drawnow limitrate
pause (.001)
end

WEEK FIVE
07/08/22  07/15/22
UPDATE: This week, we worked on various exercises surrounding the differential equation(s) and concept we were introduced with last week. I learned further MATLAB/Ocatve skills while also expanding my physics/mathematical understanding of the topics we were undertaking.
EXERCISE #8  MODELLING WITH FIXED ENDS
Using the midpoint approximation of a second derivative, we managed to derive a formula that we utilized to simulate a wave [(r^2.*y(x+dx,t)+y(xdx,t)) +(2*(1r^2).*2y(x,t)y(x,tdt)]. In this wave, we also had the initial shape of 1/2[1cos(2*pi*x / L)] with the slope and velocity being zero when t =0 and the wave having fixed ends (at 0).
The animation produced is linked below, and following that is the code to this specific animation.
WAVE SIMULATION (fixed ends)
% diff eq code
clear all; close all; clc
% parameters
dx = 0.5;
dt = dx / 10;
v = 5;
l = 10;
r = (dt/dx) * v;
n = (l/dx) + 1;
final_time = 10;
tstep = final_time/dt;
% set up initial conditions
current = 0.5*(1cos((2*pi)/l.*[0:dx:l])); % boundary
past = current;
save = [];
% loop
x_space = [0:dx:l];
for t = 0:dt:final_time
future(1) = 0;
future(2:n1) = r^2.*(current(3:n)+current(1:n2))+(2*(1r^2).*current(2:n1))past(2:n1);
future(n) = 0;
% set up for next tstep
save = [save;current];
past = current;
current = future;
end
final = size(save)
% plotting & animation
for k = 1:final(1)
plot((x_space),(save(k,:)));
axis([0,10,1,1]);
hold on
pause(.05)
Z(k) = getframe(gcf);
hold off
end
% saving video
videofile = VideoWriter('SSI.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,Z)
close(videofile)

Following this concept, we also experimented with different boundary conditions, for which the only things that were changed in the code above were future(1) and future(n) to equal the two new boundary conditions and setting the initial condition to be current = zeros(length(x_space)). Because we were still dealing with fixed end points, not much else had to be changed. Linked below are a few different videos simulating different boundary conditions. This was somewhat meant to resemble an instance with two children spinning a jumprope, in which each kid had their own end of the rope and was spinning it in their own manner. For one of the videos, both ends had the same boundary condition/wave, while in other videos, they had completely differing waves.
Fixed Ends [sin(2t) and cos(15t)] Fixed Ends Wave [sin(2t)] Fixed Ends Simulator [sin(2t), cos(2t)]
EXERCISE #9  MODELLING WITH NOT FIXED ENDS
Using most of the same concepts of the previous exercise, modelling a wave which did not have fixed ends meant altering the initial conditions and the boundary conditions. However, the derivation remained relatively similar with the exception that there were no fixed ends, which made one portion of the boundaries, y(dx,t) problematic as there is not a thing as dx in this concept, as well as something similar with the second boundary condition accounting for slope = 0 at x = L.
The changes to the boundary conditions resulted in these new boundary condition pseudocodes.
%when x = 0, slope = 0
2*r^2.* (current(2))+(2*(1r^2)*current(1))past(1);
%when x = L, slope = 0
2*r^2.* (current(n1))+(2*(1r^2)*current(n))past(n);
Likewise, the code will illustrate a new condition of initial_x, whic is utilized to allow for the initial conditions to create a dip in only one part of the wave, even before the animation begins. This is not necessarily needed to illustrate the not fixed boundaries, as the prior initial conditions of exercise #8 should be able to work, but they would not create the dip we were wishing to see.
The animation produced is linked below, and following that is the code to this specific animation.
Jumprope Simulation
%%%%%%%%%% JUMP ROPE WAVE SIMULATOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all; clc
% parameters
dx = 0.1;
dt = dx/5;
v = 1;
l = 10;
r = (dt/dx) * v;
n = (l/dx) + 1;
final_time = 10;
tstep = final_time/dt;
x_space = [0:dx:l];
initial_x = linspace(0,2*pi,((4/dx)(3/dx)+1));
% set up initial conditions
%current = 0.5*(1cos((2*pi)/l.*[0:dx:l])); % initiating boundary
%current = 0.5+0.5(cos(initial_x));
current = zeros(1,n);
current((3/dx):(4/dx)) = 0.5+0.5*(cos(initial_x));
past = current;
save = [];
% loop
for t = 0:dt:final_time
future(1) = 2*r^2.* (current(2))+(2*(1r^2)*current(1))past(1);
future(2:n1) = r^2.*(current(3:n)+current(1:n2))+(2*(1r^2).*current(2:n1))past(2:n1);
future(n) = 2*r^2.* (current(n1))+(2*(1r^2)*current(n))past(n);
% set up for next tstep
save = [save;current];
past = current;
current = future;
end
final = size(save);
% plotting & animation
for k = 1:final(1)
plot((x_space),(save(k,:)));
axis([0,10,2,2]);
hold on
pause(.1)
Z(k) = getframe(gcf);
hold off
end
% saving video
videofile = VideoWriter('nonfixed.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,Z)
close(videofile)

goals for the weekend and upcoming week:
1. experiment with new boundary and initial conditions
WEEK SIX & SEVEN
Wave Motion of Flexible Strings  Realistic Modeling
UPDATE: These two weeks, we've focused on a more realistic project which includes past research regarding the wave motion of flexible strings using linear elastic theory. What we are focusing on, primarily, is the modelling of this prior research project while including damping and a vibrating end.
The two models we focused on were one with unrealistic assumptions and another that had more realistic assumptions for small amplitudes. The first model had motion only in the y direction, had clamped ends, and utilized d^2y/dt = CT^2 * d^2y (x,t) /dx to model it. The second model had motion in both the x and y direction, as well as a vibrating end. Because it had motion in the y direction, the prior equation was still utilized, along with d^2u/dt^2 = µdu/dt  CL^2 *d^2u/dx^2 = (CL^2  CT^2)y''y'.
To solve and simulate this model, we utilized finite difference methods.
We were able to derive an equation that including damping by combining four individual equations [ d^2 u / dx^2  CL^2 d^2u/dt^2 + KL du/dt = (CL^2  CT^2) dy/dx * d^2y/dx^2] to solve for u(x, t+dt). This ended up coming out to be (2u(x,t)(1R^2)+u(x,tdt)(KL *dt/2 1) +R^2 (U(x+dx,t) + u(xdx,t)) + dt^2(CL^2CT^2)/2dx^2 (y(x+dx,t)  y(xdx,t))(y(x+dx,t)2y(x,t)+y(xdx,t)))/(1+ KL*dt/2)) with R being (dt/dx) * CT.
Implementing this a new code that resmebled the code from the prior week, we were able to adjust the parameters, boundary conditions, and initial conditions to model the experiment done in the research conducted the prior year. We were also able to trace a single point on the line. This is seen in this video: Wave with Stretching
The code to create this model is listed below  with slight modifications to be able to trace three points on the line (one on the left, one on the right, and one in the middle).
%%%%%%%%%% 07/29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all; clc
%%%% parameters %%%%
KT = 20; % coefficient of damping along transverse direction
KL = 20; % coefficient of damping along longitudinal direction
LI = 0.05; % Initial length [m]
L0 = 5*LI; % length
SP = 5; % Spring Constant
e = (L0/LI)1; % strain on the rope
T0 = SP * LI * e; % tension
m = 0.0025; % mass
CT = sqrt(T0 * L0/ m); % velocity in the transverse direction
CL = CT * sqrt((1+1/e)); % velocity in the longitudinal direction
n = 101;
k = pi/L0;
alpha = k*CT;
x_space = linspace(0,L0,n);
dx = L0/(n1);
dt = dx/(2*CT);
A = LI; % amplitude
%%%% set up initial conditions %%%%
current = zeros(1,n); % flat
%current = A * sin(x_space*k); ===>> current = 0.05*sin((pi*x_space)/l);
%current = 0.5*(1cos((2*pi)/l.*[0:dx:l]));
past = current; % rest
%past = current  (A*alpha* sin(x_space*k))*dt;
ucurrent = zeros(1,n); % flat
%ucurrent = A^2*k/16*sin(2*k*x_space)*(CT/CL)^2;
upast = ucurrent; %at rest
usave = [];
save = [];
%%%% loop %%%%
final_time = 20*pi/alpha;
dt = dt;
r = (dt/dx) * CT; %parameter
P = (dt^2 * (CL^2CT^2)/(2 * dx^3)); % parameter
for t = 0:dt:final_time
% solution for displacement along the y direction
future(1) = 0.008 * sin(alpha * t);
%future(1) = 0; % for clamped
future(2:n1) = (r^2.*(current(3:n)+current(1:n2))+(2*(1r^2).*current(2:n1))+ past(2:n1) .* (KT * dt / 2  1)) ./ (1 + (KT * dt / 2));
future(n) = 0;
% solution for displacement along the u direction
ufuture(1) = 0;
ufuture(2:n1) =(2 .* ucurrent(2:n1) .* (1r^2) + upast(2:n1) .* (KL * dt / 2  1) + r^2.*(ucurrent(3:n)+ucurrent(1:n2)))/(1 + (KL * dt / 2))...
+ (P.* (current(3:n)current(1:n2)) .* (current(3:n)  2 * current(2:n1) + current(1:n2)))/(1 + (KL * dt / 2));
% ufuture(2:n1) =
%ufuture(2:n1) = (dt^2/(2*(KL1)).* ((((CL^2CT^2)/dx^2) .* current(3:n)  2*current(2:n1) + current(1:n2) .* (dx (current(3:n)current(2:n1))))  ((1/dx^2) .* ucurrent(3:n)  2* ucurrent(2:n1) + ucurrent(1:n2)) ((2/dt^2) .* ucurrent(2:n1)) + ((1/dt^2) .*upast(2:n1))+ ((KL/dt) .* ucurrent(2:n1))));
ufuture(n) = 0;
%ufuture(2:n1) =((2 .* ucurrent(2:n1) .* (1r^2) + upast(2:n1) .* (KL * dt / 2  1) + r^2.*(ucurrent(3:n)+ucurrent(1:n2))) + P.* (future(3:n)future(1:n2)) .* (future(3:n)  2 * future(2:n1) + future(1:n2)))/(1 + (KL * dt / 2));
%ufuture(2:n1) =(2 .* ucurrent(2:n1) .* (1r^2) + upast(2:n1) .* (KL * dt / 2  1) + r^2.*(ucurrent(3:n)+ucurrent(1:n2)) + (CL^2CT^2)*dt^2*(A^2*k^3/2*(cos(alpha*(t)))^2*sin(2*k*x_space(2:n1))))/(1 + (KL * dt / 2));
%ufuture(2:n1) =(2 .* ucurrent(2:n1) .* (1r^2+KL/2*dt) + upast(2:n1) .* (1) + r^2.*(ucurrent(3:n)+ucurrent(1:n2)) + 2*P.* (current(3:n)current(2:n1)) .* (current(3:n)  2 * current(2:n1) + current(1:n2)))/(1 + (KL * dt));
%
% set up for next tstep
save = [save;current];
past = current;
current = future;
usave = [usave;ucurrent];
upast = ucurrent;
ucurrent = ufuture;
end
final = size(save);
%%%% plotting & animation %%%%
w = alpha;
T = linspace(0, 2*pi/w, 40);
K = pi/L0;
xm = (n  1)/2 + 1; % center
xl = (n  1)/4 + 1; % left
xr = 3*(n  1)/4 + 1; % right
I = round(final(1)/2); % initial point
for k = I:10:final(1)
plot ((x_space), save(k,:),'b');
%axis ([0,L0,A,A]); %%% this will show just one point zoomed in
axis([L0/4L0/10,L0/4+L0/10,A,A]); %% this will show all three points
hold on
plot((x_space + usave(k,:) ),save(k,:),'r' ...
,(x_space(xl) + usave(k,xl) ),save(k,xl),'ro',(x_space(xl) + usave(I:k,xl) ),save(I:k,xl),'r' ...
,(x_space(xm) + usave(k,xm) ), save(k,xm),'ro',(x_space(xm) + usave(I:k,xm) ),save(I:k,xm),'r' ...
,(x_space(xr) + usave(k,xr) ), save(k,xr),'ro',(x_space(xr) + usave(I:k,xr) ),save(I:k,xr),'r');
u_L_vector = A^2*K/16.*sin(2*K*x_space(xl)).*(cos(2*w*T) + 1/(e + 1));
y_L_vector = A * sin(K*x_space(xl)).*cos(w*(T));
plot (u_L_vector +x_space(xl),y_L_vector,'k')
pause(.0001)
%Z(k) = getframe(gcf);
hold off
end
% saving video
%videofile = VideoWriter('wavesim2.avi','Uncompressed AVI');
%open(videofile)
%writeVideo(videofile,Z)
%close(videofile)

Using this, we were also able to discover the difference different amounts of damping can cause. If the coefficient of damping along longitudinal and transverse direction is 20, an 8 shape will be made from tracing a single point. But if the coefficient of damping along longitudinal and transverse direction is 50, a"backwards" cresent will be made. This is visible in the photos below which illustrate the final outcome of the model after the code is run once (with 50 and 20 as the coefficient of damping).
in this case, the black line visible in both photos is what was somewhat "estimated" to happen per the prior research and analytical solution. however, even after using both the backwards and forward approximation, this did not seem to be the case
Based on this week, we concluded that the introduction of damping in the system potentially causes it to lag and then overshoot. And it then tries to "correct" this overshoot, but only continues this cycle, in this case "undershooting," which it then tries to, again, correct by overshooting. This cycle is seen in the right image (potentially being the reason for the "8" shape).
upcoming week: we will be focusing on the second model a bit further with another method that includes focusing on individual atoms.