Euler Method Simulations
project updates and summaries
WEEK ONE
06/11/22 - 06/17/22
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 Code | Image #2 Code |
|
% position as a function of time -- euler's method % initial parameters dt = 0.01; %interval between steps
|
% including energy comparison % Euler's Method - Numerical and Analytical Comparison
plot(time,energy,'c',time,e_exact,'r')
|
WEEK TWO
06/18/22 - 06/24/22
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].
for i = 2:tstep % analytical method v_exact = -y1 * sqrt(k./ m) .* sin(sqrt(k ./ m) .* time); % plot results of numerical and analytical computations
%ploting position 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]. for i = 2:tstep theta_exact = theta(1) .* cos(((3 * g)/(2 * L)) .* time); title('Small Angle Approx VS Euler-Cromer')
|
SLIGHT ALTERATIONS % parameters k = 1; %spring constant [N/m].
theta_exact = theta(1) .* cos(((3 * g)/(2 * L)) .* time); title('Euler vs Exact')
|
image # 6 image #7 image #8
WEEK THREE
06/24/22 - 07/01/22
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 % loop for changing amplitude for i = 2:tstep [PKS,LOC] = findpeaks(abs(theta')); plot(in_am, period, 'x') 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.
RIGID PENDULUM ANIMATION - YOUTUBE
|
% code models and animates motion of rigid pendulum
k = 1; %spring constant [N/m]. for i = 2:tstep theta_exact = theta(1) .* cos(((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')
% pendulum animation values X = -.5*L.*sin(theta); % center of mass x direction figure % animation for pendulum loop for n = 1:length(theta) % storing frames videofile = VideoWriter('pendulum.avi', 'Uncompressed AVI') |
WEEK FOUR
07/01/22 - 07/08/22
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
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 % parameters % set up initial conditions save = []; % loop x_space = [0:dx:l]; for t = 0:dt:final_time final = size(save) % plotting & animation for k = 1:final(1) % saving video videofile = VideoWriter('SSI.avi','Uncompressed AVI'); |
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.
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.
|
%%%%%%%%%% JUMP ROPE WAVE SIMULATOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; clc % parameters x_space = [0:dx:l]; % set up initial conditions save = [];
final = size(save); % plotting & animation for k = 1:final(1) % saving video videofile = VideoWriter('nonfixed.avi','Uncompressed AVI'); |
1. experiment with new boundary and initial conditions
WEEK SIX & SEVEN
Wave Motion of Flexible Strings - Realistic Modeling
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 LI = 0.05; % Initial length [m]
A = LI; % amplitude %%%% set up initial conditions %%%% current = zeros(1,n); % flat ucurrent = zeros(1,n); % flat usave = [];
% solution for displacement along the y direction future(1) = 0.008 * sin(alpha * t); 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))) + P.* (future(3:n)-future(1:n-2)) .* (future(3:n) - 2 * future(2:n-1) + future(1:n-2)))/(1 + (KL * dt / 2)); final = size(save); %%%% plotting & animation %%%% xm = (n - 1)/2 + 1; % center
u_L_vector = -A^2*K/16.*sin(2*K*x_space(xl)).*(cos(2*w*T) + 1/(e + 1)); % saving video %videofile = VideoWriter('wavesim2.avi','Uncompressed AVI'); |
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).