 # Euler Method Simulations

June 13, 2022, 8:29 p.m.
abeyret

### 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) = (i-1)*dt; v_Euler(i)=v_Euler(i-1)-k*y_Euler(i-1)*dt/m; y_Euler(i)=y_Euler(i-1)+v_Euler(i-1)*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(i-1) + dt; % adds time interval to previous time   v(i)      = v(i-1) - (k./m) .* y(i-1) .* dt; % calculates new velocity   y(i)      = y(i-1) + v(i-1) * 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');

### 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 euler-cromer 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 EULER-CROMER METHOD

With the Euler-Cromer 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 Euler-Cromer 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(i-1) + dt; % adds time interval to previous time   force(i)        = -k.*y(i-1);   a(i)            = force(i)/m;   v(i)            = v(i-1) + a(i) * dt; % final velocity = v1 +acceleration(time)   y(i)            = y(i-1) + 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 euler-cromer 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(i-1) + dt; % time interval   om(i)           = om(i-1)+dt*a(i-1); % omega interval   theta(i)        = theta(i-1)+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(i-1) +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 Euler-Cromer') xlabel('time'); ylabel('position'); plot (time, theta,'r', time, theta_exact, 'b') legend ('Small Angle Approx', 'Euler-Cromer') figure plot (theta, om,'g') hold on plot (theta_exact, om_exact, 'r') legend ('Small Angle Approx', 'Euler-Cromer') 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(i-1) + dt; % adds time interval to previous time   om(i)           = om(i-1)+dt*a(i-1);   theta(i)        = theta(i-1)+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(i-1) +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

### 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(i-1) + dt; % time interval   om(i)           = om(i-1)+dt*a(i-1); % omega interval   theta(i)        = theta(i-1)+dt *om(i); % theta interval   a(i)  =   -theta(i)^3;   t(i) = time(i-1) +dt; end [PKS,LOC] = findpeaks(abs(theta')); LOC; newtime = time(LOC); period(n) = 2* (mean(diff(newtime))); end plot(in_am, period, 'x') legend ('Euler-Cromer') title('Euler-Cromer')  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.

 % 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(i-1) + dt; % time interval   om(i)           = om(i-1)+dt*a(i-1); % omega interval   theta(i)        = theta(i-1)+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(i-1) +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 Euler-Cromer') %xlabel('time'); %ylabel('position'); %plot (time, theta,'r', time, theta_exact, 'b') %legend ('Small Angle Approx', 'Euler-Cromer') %figure %plot (theta, om,'g') %hold on %plot (theta_exact, om_exact, 'r') %legend ('Small Angle Approx', 'Euler-Cromer') % 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)

### 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 cup-like 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

### 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(x-dx,t)) +(2*(1-r^2).*2y(x,t)-y(x,t-dt)]. In this wave, we also had the initial shape of 1/2[1-cos(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.

 % 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*(1-cos((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:n-1) = r^2.*(current(3:n)+current(1:n-2))+(2*(1-r^2).*current(2:n-1))-past(2:n-1);     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 pseudo-codes.

%when x = 0, slope = 0
2*r^2.* (current(2))+(2*(1-r^2)*current(1))-past(1);
%when x = L, slope = 0
2*r^2.* (current(n-1))+(2*(1-r^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*(1-cos((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*(1-r^2)*current(1))-past(1);     future(2:n-1) = r^2.*(current(3:n)+current(1:n-2))+(2*(1-r^2).*current(2:n-1))-past(2:n-1);     future(n) = 2*r^2.* (current(n-1))+(2*(1-r^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

### 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)(1-R^2)+u(x,t-dt)(KL *dt/2 -1) +R^2 (U(x+dx,t) + u(x-dx,t)) + dt^2(CL^2-CT^2)/2dx^2 (y(x+dx,t) - y(x-dx,t))(y(x+dx,t)-2y(x,t)+y(x-dx,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/(n-1); 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*(1-cos((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^2-CT^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:n-1) = (r^2.*(current(3:n)+current(1:n-2))+(2*(1-r^2).*current(2:n-1))+ past(2:n-1) .* (KT * dt / 2 - 1)) ./ (1 + (KT * dt / 2));     future(n) = 0;          % solution for displacement along the u direction     ufuture(1) = 0;     ufuture(2:n-1) =(2 .* ucurrent(2:n-1) .* (1-r^2) + upast(2:n-1) .* (KL * dt / 2 - 1) + r^2.*(ucurrent(3:n)+ucurrent(1:n-2)))/(1 + (KL * dt / 2))...    + (P.* (current(3:n)-current(1:n-2)) .* (current(3:n) - 2 * current(2:n-1) + current(1:n-2)))/(1 + (KL * dt / 2));    % ufuture(2:n-1) =     %ufuture(2:n-1) = (dt^2/(2*(KL-1)).* ((((CL^2-CT^2)/dx^2) .* current(3:n) - 2*current(2:n-1) + current(1:n-2) .* (dx (current(3:n)-current(2:n-1)))) - ((1/dx^2) .* ucurrent(3:n) - 2* ucurrent(2:n-1) + ucurrent(1:n-2))- ((2/dt^2) .* ucurrent(2:n-1)) + ((1/dt^2) .*upast(2:n-1))+ ((KL/dt) .* ucurrent(2:n-1))));     ufuture(n) = 0;    %ufuture(2:n-1) =((2 .* ucurrent(2:n-1) .* (1-r^2) + upast(2:n-1) .* (KL * dt / 2 - 1) + r^2.*(ucurrent(3:n)+ucurrent(1:n-2))) + P.* (future(3:n)-future(1:n-2)) .* (future(3:n) - 2 * future(2:n-1) + future(1:n-2)))/(1 + (KL * dt / 2));    %ufuture(2:n-1) =(2 .* ucurrent(2:n-1) .* (1-r^2) + upast(2:n-1) .* (KL * dt / 2 - 1) + r^2.*(ucurrent(3:n)+ucurrent(1:n-2)) + (CL^2-CT^2)*dt^2*(-A^2*k^3/2*(cos(alpha*(t)))^2*sin(2*k*x_space(2:n-1))))/(1 + (KL * dt / 2));    %ufuture(2:n-1) =(2 .* ucurrent(2:n-1) .* (1-r^2+KL/2*dt) + upast(2:n-1) .* (-1) + r^2.*(ucurrent(3:n)+ucurrent(1:n-2)) + 2*P.* (current(3:n)-current(2:n-1)) .* (current(3:n) - 2 * current(2:n-1) + current(1:n-2)))/(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/4-L0/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.