HW2

Aleksandr Mnatsakanyan HW#2 Problem 1: see rocket_euler.m Problem 2: Here is the plot that I got from the script. We can see that smaller the “step” of the calculation, smaller the error of the calculation but looking at the plot there is a limited benefit of decreasing the size of the delta as the accuracy gained in approximation will not be relative to the computation time. Problem 3: let's assume that the origin (s = 0) is the position of the velociraptor, so that the displacement of the unfortunate human at t = 0 is s h = 40 m . As the acceleration of the velociraptor is constant we can use suvat motion equations to solve for its time and displacement. v 2 = v i 2 + 2 as As v i is 0 , s can be derived to be

s = ¿ v 2 2 a , s = ¿ 25 2 2 ⋅ 4 ¿ 78.125 m Using v = v i + at , we get t = ¿ v a ¿ 25 4 ¿ 6.25 s , that is the time that the velociraptor takes to reach its max speed of 25 m/s. On the other hand the human in that period of time only reaches total displacement of, s = 40 m + 6.25 s ⋅ 6 m / s = 77.5 m which is less than the displacement of the velociraptor before reaching its top speed, therefore to find at which time we can simply equate the total displacement equation for the human with the velociraptor s = v i t + 0.5 at 2 to get. v h t + 40 = 0.5 at 2 Substituting the values and solving the polynomial for t we get t = 6.217 s , therefore the distance that the human traveled before being devoured is s = v h t = 6 ⋅ 6.217 ≈ 37.3 m from their original position. Problem 4: see hw2_p4.m. As we can see in the plot the gator intercepts the human a bit later than 6 seconds which is very close to the analytical solution Now lets have fun with it. I will change the acceleration to 2 m s 2 and the maximum speed of the gator to 11 m s . As after a slight google search these values are a more accurate

2 1 2^-1 0.5 3 0 2^0 0 4 0 2^1 0 5 1 2^2 4 6 1 2^3 8 7 1 2^4 16 8 0 2^5 0 9 1 2^6 64 10 0 2^7 0 11 0 2^8 0 So the number ( 00101110011 ) 2 is ( 64 + 16 + 8 + 4 + 0.5 + 0.25 )=( 92.75 ) 10 Problem 6: see hw2_p6.m I have written a script to automate this process as I was too lazy to do this by hand. The final binary value that I got was 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 1 1 Where it starts with the answer of question 1 and ends with question 20. After searching who could be a person obsessed with math, be an astronomer and have persian origins, there were too many options so I went to my old friend Akinator, and strangely he was asking questions of a familiar style. His name is Muhammad ibn Musa al-Khwarizmi Problem 7: see hw_p7.m As found by the script the value that I got using the algorithm is the same as in the built in function. Problem 8: see hw_p8.m As mentioned in the script, here is how I found the formula for the Taylor series. We know that the formula for the nth order variable is nth variable : f n ( a ) ¿¿ Now after looking at the derivatives of this specific logarithm we get f n ( x )= x − 1 , − x − 2 , 2 x − 3 , 6 x − 4 ,..., ¿ After substituting and canceling the factorial terms as ( n! = n ⋅ ( n − 1 ) ! ) we get

nth variable : − 1 n + 1 ¿¿ ;n ≥ 1 And here is the plot up to 15th degree TSA.

the last value of the acceleration vector, on the other hand backwards approximation is gives smaller values compared to the analytical solution so the mean error is positive(0.858) but cannot be used to approximate the first value of the acceleration vector. “Hold my beer” says centered approximation with mean error of (-0.048) even though the first and last instances are not possible to find with this method, it is very precise as discussed in class. Problem 10: see lorenz.m First plot Second plot As we can see these two plots are slightly different even though we used the values generated from the initial calculations from the first plot, but due to this system of equations having sensitive dependance on initial conditions even a slight change(rounding off mostly in this case) will change the outcome, aka “butterfly effect”. Hw2_p1 function [t,v,m] = rocket_euler(delt,simend,params) %%This function simulates a rocket's velocity during launch using Euler's

%method. %%Inputs: delt - step size for Euler's method (s) % simend - final simulation time (s) % params - launch parameters (structure variable, with fields c, m0, and mdot) %%Outputs: t - time vector for simulation (s) % v - velocity vector of simulation (m/s) % m - mass vector of simulation (kg) %%Author: Aleksandr Mnatsakanyan, University of Florida %Get the variables we need out of the params structure variable c = params.c; %thrust coefficient (m/s) m0 = params.m0; %initial mass (kg) mdot = params.mdot; %fuel burn rate (kg/s) g = 9.81; count = 2; % used as a counter for the vector indexes in the for loop v_p = 0; % intermediate velocity variable used to create the vector with all of the velocities t = zeros(1,(simend/delt)+1); % time variable used as the h in eulers (later found this a better method which proved this one very ineffective, but this worked so did not change it) m = zeros(1,(simend/delt)+1); % also preallocating for mass m(1) = m0; % first index of the mass vector has to be the initial mass accel = zeros(1,(simend/delt)+1); % another preallocation v = zeros(1,(simend/delt)+1); % and another one for ii = delt:delt:simend % not stating with 0 as that is already given by the initial values t(count) = ii; % assigning t values for each iteration of the eulers calculation accel(count) = -g - (c*mdot)/m(count-1); % calculating dv/dt at t sec v(count) = v_p + accel(count)*delt; % eulers formula in action, wow v_p = v(count); % assigning the old value to this intermediate variable m(count) = m0 + mdot*t(count); % calculating new mass count = count + 1; % updating the count for the next index array end end Hw2_p2 clc;clear; %This script uses the rocket_euler function to plot a graph of multiple

%This script uses the initial conditions given to solve the problem using eulers method %The script also plots the displacement of the human and the displacement of the velociraptor dis = 40; % heads-up distance of the human (meters) a = 2; % acceleration of the gator (m/s^2) vmax = 11; % maximum speed of the gator (m/s) vhum = 6; % maximum speed of the human (m/s) delt = 0.1; % delta used (s) tend = 600; % potential end time of the simulation (s) euler = 0:delt:tend; % time vector used for convenience sGator = zeros(1,(tend/delt)); % preallocating for the final variables sHuman = zeros(1,(tend/delt)); sHuman(1) = dis; % giving the human the additional 40 meters count = 2; % counter, used to control indexes throughout the script for t = euler v = a*t; % finding the velocity at time t if v > vmax % if its bigger than 25 then assign it to be 25 therefore the maximum value it can be is 25 v = vmax; end sHuman(count) = sHuman(count-1) + vhum*delt; % as human speed is a constant then its displacement is going to be a line in eulers method as well sGator(count) = sGator(count-1) + v*delt; % using eulers method we can approximate the position of the gator even if it has changing velocity if sGator(count)>= sHuman(count) % here we check if the displacement of the gator is greater than the displacement of the human if so we stop the eulers calculation break ; end count = count + 1; % updating the count variable end sHuman(length(sHuman):-1:(count+1)) = []; % due to the setup of the vectors the last n indixes need to be removed as they are zeros due to preallocation so that we only plot the relevant plot sGator(length(sGator):-1:(count+1)) = []; % we do that by setting n = (size of the vector) - (the index of the first non calculated zero) euler(length(euler):-1:(count+1)) = []; figure plot(euler,sHuman, '-*' ) hold on plot(euler,sGator, '-o' ) hold off legend( 'Human' , 'Oversized gator' ) title( 'Very Vital Plot' ) subtitle( 'Displacement vs Time' ) xlabel( 'Time (sec)' ) ylabel( 'Displacement (m)' ) Hw2_p6

clc; clear; % This script converts any base 10 number to a 20 digit binary number base = 29953.03125; % given base 10 value binary = zeros(1,20); % preallocating for 20 values for ii = 20:-1:1 % reverse for loop to be able to subtract first smaller value from the base current = 2^(ii-6); % last digit corresponds with 2 to the power of 14 as the question states that the first digit is 2^-5 if base - current >= 0 % checking if the current power of 2 is smaller than the reminder of the previous successful calculation binary(ii) = 1; % if so then that index in the binary value shold change from 0 to 1 base = base - current; % new base value end end Hw2_p7 clc; clear; %This script will be used to find the smallest number that can be processed %by the computer and compare it to the built in value epsilon = 1; while true if 1+epsilon <= 1 break else epsilon = epsilon/2; end end epsilon = 2 * epsilon; %2.220446049250313e-16 true_eps = eps(); %2.220446049250313e-16 Hw2_p8

clc; clear; %This script takes the data samples from the file provided and calculates %the forward, backward and centered finite difference approximation and the %error that each method gives. values = load( "hw2_p9.mat" ); % loading the values from the given file t = 0:1:120; % time vector used for calculation params = struct( 'c' , 4500, 'm0' , 3*10^6, 'mdot' , -10000); % constants h = values.t(2,1) - values.t(1,1); % delta t accelf = (values.v(2:7,1)-values.v(1:6,1))./h; % forward - finite difference approximation calculation accelb = (values.v(2:7,1)-values.v(1:6,1))./h; % backward - finite difference approximation calculation accelc = (values.v(3:7,1)-values.v(1:5,1))./(2*h); % centered - finite difference approximation calculation m = params.m0 + params.mdot.*t; % getting the mass vector for all t times a_analytical_plot = -9.81 - params.c.*params.mdot./(params.m0+params.mdot.*t); % getting the acceleration vector for the analytical plot a_analytical_deriv = -9.81 - params.c.*params.mdot./(params.m0+params.mdot.*values.t); % getting the acceleration vetor for the true derivative calculation % I decided to calculate the mean error from all the points in those three % approximations to have a relative understanding of the accuracy that each % of them provides front_error = sum((a_analytical_deriv(1:6) - accelf)./size(accelf,1)); back_error = sum((a_analytical_deriv(2:7) - accelb)./size(accelb,1)); center_error = sum((a_analytical_deriv(2:6) - accelc)./size(accelc,1)); % printing the mean errors fprintf( 'FFD mean Error: %g \nBFD mean Error: %g \nCFD mean Error: %g' ,front_error,back_error,center_error) figure hold on plot(t,a_analytical_plot) plot(values.t(1:6),accelf, "-o" ) plot(values.t(2:7),accelb, "-x" ) plot(values.t(2:6),accelc, "-*" ) hold off title( 'Front, Backward and Centered finite approximation plot' ) subtitle( 'a vs t' ) xlabel( 'Time (sec)' ) ylabel( 'Acceleration (m/s^2)' ) legend( 'Analytical' , 'front' , 'backwards' , 'center' ) Hw2_p10

clc; clear; %This script is used to solve the lorenz differential equation system and %plot the resulting values delt = 0.01; % delta t used in the calculation t_start = 500; % starting t value t = t_start:delt:1000; % time vector used in the calculation and preallocation % preallocation of all the variables used x = zeros(1,size(t,2)); y = zeros(1,size(t,2)); z = zeros(1,size(t,2)); xdot = zeros(1,size(t,2)); ydot = zeros(1,size(t,2)); zdot = zeros(1,size(t,2)); count = 1; % variable used to control the indexes of the vectors created eariler % initial conditions for x,y and z, and calculation of all the derivatives % for the initial conditions x(count) = 4.8072; y(count) = 8.4159; z(count) = 12.6007; xdot(count) = 10*(y(count)-x(count)); ydot(count) = x(count)*(28-z(count))-y(count); zdot(count) = x(count)*y(count)-(8/3)*z(count); for i = t % calculation of the next values for x,y and z using eulers method x(count+1)= x(count) + xdot(count)*delt; y(count+1)= y(count) + ydot(count)*delt; z(count+1)= z(count) + zdot(count)*delt; count = count + 1; % calculation of the next values for all the derivatives using the % equations in the system given xdot(count) = 10*(y(count)-x(count)); ydot(count) = x(count)*(28-z(count))- y(count); zdot(count) = x(count)*y(count)-(8/3)*z(count); end % adding one more instance of time for plotting purposes t(count) = t(count-1) + delt; figure %first plot is the plot of x,y and z vs t from t(995) to t(1000) subplot(1,2,1) plot(t,x,t,y,t,z) xlim([995 1000]) % limiting the x axis with the conditions above title( 'Lorenz differential equations' ) subtitle( 'x, y, z vs t' ) xlabel( 't' ); ylabel( 'x,y and z' ) legend( 'x' , 'y' , 'z' ) subplot(1,2,2) % second plot is the plot of z vs x to show the famous butterfly attractor % from t(900) to t(1000) as I found these limits to work best plot( x(size(t,2)-10000 : size(t,2)), z(size(t,2)-10000 : size(t,2))); xlim([-25 25]) ylim([0 55]) title( 'Lorenz differential equations' ) subtitle( 'strange attractor' ) xlabel( 'x' ); ylabel( 'z' ) legend( 'T(900)-T(1000)' )