Final Project Problem 2: Numerical Integration

Problem E

Analytical Approach

Solution Steps

First, I took the indefinite integral of equation E1 with respect to time to determine the formula for the distance of the rocket. Then, the values of 5 s, 10 s, and 30 s were inserted into the resulting distance equation to find the height of the rocket.

Analytical Results

Graphical Representation of Problem

Figure 1: The figure above shows the velocity over time function of the rocket.

The integral of the velocity function after 5 seconds is 238.59 m. This is the analytical height of the rocket after 5 seconds.

The integral of the velocity function after 10 seconds is 995.81 m. This is the analytical height of the rocket after 10 seconds.

The integral of the velocity function after 30 seconds is 10880 m. This is the analytical height of the rocket after 30 seconds.

Numerical Approach

Solution Steps

Part A)

The integral of the velocity function of the rocket was approximated using the six segment trapezoidal rule, Simpson's 1/3 Rule, and O(h^8) Romberg Method. To calculate the O(h^8) Romberg Solution, first the O(h^4) Romberg Solution and the O(h^6) Romberg Solution must be calculated from trapezoidal rule approximations. Therefore, I computed, one, two, four, and eight segment trapezoidal rule approximations as well.

Part B)

The acceleration of the rocket was found by approximating the derivative of the velocity function using finite divided differences. Since the times of 5, 10, and 30 seconds were significant in Part A, the derivative at the points from 0 to 30 with a step size of 5, were estimated.

Numerical Results

Part A)

The six segment trapezoidal yielded an integral of 238.73 m, which is the approximate height of the rocket after 5 seconds. After 10 seconds, the approximate height of the rocket is 997.02 m. Lastly, the approximate height of the rocket after 30 seconds using the six segment trapezoidal rule is 10931 m.

Using Simpson's 1/3 Rule, the integral of the velocity function after 5 seconds was 238.59 m. After 10 seconds, the approximate height of the rocket is 995.84 m. Lastly, after 30 seconds, the approximate height of the rocket using Simpson's 1/3 Rule is 10897 m.

Using the O(h^8) Romberg Method, the integral of the velocity function after 5 seconds was 747.76 m. After 10 and 30 seconds, the approximate height of the rocket is 3120.9 m and 34099 m, respectively.

Part B)

The results of the derivative approximations using finite divided differences can be seen in Table 2.

Figure 2: The acceleration of the rocket over time is plotted above.

Error Analysis

Table 1: The percent relative error for method to approximate the integral of the velocity function is shown below.

Table 2: The table below shows the analytical acceleration calculated by taking the derivative of equation E1, the numerical acceleration approximated using finite divided differences, and the percent relative error between the two.

Discussion

Part C)

Simpson's 1/3 Rule produced the lowest error for the height of the rocket. The magnitude of the error was over 100 times larger for the 5 second interval for the trapezoidal rule compared to the Simpson's Rule approximation. I would prefer to use Simpson's rule in this problem due to its higher accuracy and the simplicity of the algorithm.

Finite divided differences proved to not be an effective method for approximating the derivative of the velocity function. The error of the approximation ranged from 77-80% for each time value. Therefore, I would recommend using a built-in coding function like MATLAB's diff function to find the analytical acceleration.

Pseudocode

Set Parameters: g, q, m0, uCreate function handle for v
Analytical SolutionUse the integral function to calculate the
Plot Analytical Solutiont = 0:.1:35plot(t,v(t))
Compute 6 Segment Trapezoidal Rulea = 0b = 5c = 10d = 30n = 6h5 = (b-a)/nApproximate Integral after 5sT5 = (h5/2)*(v(a)+2*(v(a+h5)+v(a+2*h5)+v(a+3*h5)+v(a+4*h5)+v(a+5*h5))+v(b))Error against analytic valuesEt5 = abs(A5-T5)et5 = (Et5/A5)*100
h10 = (c-a)/nApproximate Integral after 10sT10 = (h10/2)*(v(a)+2*(v(a+h10)+v(a+2*h10)+v(a+3*h10)+v(a+4*h10)+v(a+5*h10))+v(c))Error against analytic valuesEt10 = abs(A10-T10)et10 = (Et10/A10)*100
h30 = (d-a)/nApproximate Integral after 30sT30 = (h30/2)*(v(a)+2*(v(a+h30)+v(a+2*h30)+v(a+3*h30)+v(a+4*h30)+v(a+5*h30))+v(d))Error against analytic valuesEt30 = abs(A30-T30)et30 = (Et30/A30)*100
Simpson's 1/3 RuleEstimate height after 5 secondsB = (b-a)/2integral after 5sS5 = (b-a)*((v(a)+4*v(B)+v(b))/6)ErrorEs5 = abs(A5-S5)es5 = (Es5/A5)*100
Estimate height after 10 secondsC = (c-a)/2integral after 10sS10 = (c-a)*((v(a)+4*v(C)+v(c))/6)ErrorEs10 = abs(A10-S10)es10 = (Es10/A10)*100
Estimate height after 30 secondsD = (d-a)/2; % midpoint of intervalintegral after 30sS30 = (d-a)*((v(a)+4*v(D)+v(d))/6)ErrorEs30 = abs(A30-S30)es30 = (Es30/A30)*100
0(h^8) Romberg MethodFirst 0(h^4) and 0(h^6) must be calculated
Part B - Numerical Differentiation% Analytical Solutionta = -5:5:35; % time intervalaccel = diff(v(ta));t=0:5:35;% Plot acceleration over timefigure(2)plot(t,accel,'r-','LineWidth',2)grid on; box on;title('Acceleration of the Rocket over Time','FontSize',16)xlabel('Time (s)','FontSize',14)ylabel('Acceleration (m/s^2)','FontSize',14)
%% Numerical Solution - Finite Divided Differencest1 = 0; t2 = 5; t3 = 10; t4 = 15; t5 = 20; t6 = 25; t7 = 30;
% Forward finite difference approximationf0 = (v(t2)-v(t1))/(t2-t1);
% Centered finite difference approximationf5 = (v(t3)-v(t1))/(t3-t1);f10 = (v(t4)-v(t2))/(t4-t2);f15 = (v(t5)-v(t3))/(t5-t3);f20 = (v(t6)-v(t4))/(t6-t4);f25 = (v(t7)-v(t5))/(t7-t5);
% Backward finite difference approximationf30 = (v(t7)-v(t6))/(t7-t6);
num_diff = [f0 f5 f10 f15 f20 f25 f30];
%% Error between Numerical and AnalyticalE0 = abs((accel(1)-f0)/accel(1))*100;E5 = abs((accel(2)-f5)/accel(2))*100;E10 = abs((accel(3)-f10)/accel(3))*100;E15 = abs((accel(4)-f15)/accel(4))*100;E20 = abs((accel(5)-f20)/accel(5))*100;E25 = abs((accel(6)-f25)/accel(6))*100;E30 = abs((accel(6)-f30)/accel(7))*100;E = [E0 E5 E10 E15 E20 E25 E30];

MATLAB Code

%%% EGR 312%%% Cameron King%%% 5/8/20%%% Final Project - Problem E
% Set Parametersg = 9.8; % gravity in m/s^2q = 2500; % fuel consumption rate in kg/sm0 = 160000; % initial mass of rocket in kgu = 1800; % velocity at which fuel is expelled in m/sv = @(t) u.*log(m0./(m0-q.*t))-g.*t;
%% Analytical SolutionA5 = integral(v,0,5);A10 = integral(v,0,10);A30 = integral(v,0,30);
%% Plot t = 0:.1:35; % time intervalfigure(1)plot(t, v(t),'LineWidth', 2)grid on; box on;title('Velocity of the Rocket over Time','FontSize',16)xlabel('Time (s)','FontSize',14)ylabel('Velocity (m/s)','FontSize',14)
%% Six Segment Trapezoidal Rulen = 6; % number of segments
% Estimate height after 5 secondsa = 0; % initial time in sb = 5; % final time in sh5 = (b-a)/n; % width of trapezoid% Integral after 5s (approximate height)T5 = (h5/2)*(v(a)+2*(v(a+h5)+v(a+2*h5)+v(a+3*h5)+v(a+4*h5)+v(a+5*h5))+v(b)); % ErrorEt5 = abs(A5-T5); % true erroret5 = (Et5/A5)*100; % percent relative error
% Estimate height after 10 secondsc = 10; % final time in sh10 = (c-a)/n; % width of trapezoid% Integral after 10s (approximate height)T10 = (h10/2)*(v(a)+2*(v(a+h10)+v(a+2*h10)+v(a+3*h10)+v(a+4*h10)+v(a+5*h10))+v(c));% ErrorEt10 = abs(A10-T10); % true erroret10 = (Et10/A10)*100; % percent relative error
% Estimate height after 30 secondsd = 30; % final time in sh30 = (d-a)/n; % width of trapezoid% Integral after 30s (approximate height)T30 = (h30/2)*(v(a)+2*(v(a+h30)+v(a+2*h30)+v(a+3*h30)+v(a+4*h30)+v(a+5*h30))+v(d)); % ErrorEt30 = abs(A30-T30); % true erroret30 = (Et30/A30)*100; % percent relative error
%% Simpson's 1/3 Rule% Estimate height after 5 secondsB = (b-a)/2; % midpoint of intervalS5 = (b-a)*((v(a)+4*v(B)+v(b))/6); % integral after 5s (approximate height)% ErrorEs5 = abs(A5-S5); % true errores5 = (Es5/A5)*100; % percent relative error
% Estimate height after 10 secondsC = (c-a)/2; % midpoint of intervalS10 = (c-a)*((v(a)+4*v(C)+v(c))/6); % integral after 10s (approximate height)% ErrorEs10 = abs(A10-S10); % true errores10 = (Es10/A10)*100; % percent relative error
% Estimate height after 30 secondsD = (d-a)/2; % midpoint of intervalS30 = (d-a)*((v(a)+4*v(D)+v(d))/6); % integral after 30s (approximate height)% ErrorEs30 = abs(A30-S30); % true errores30 = (Es30/A30)*100; % percent relative error
%% O(h^8) Romberg Method% To obtain the O(h^8) estimate, 4 different approximations of the% trapezoidal rule must be computed. I will compute 1, 2, 4, and 8 segment% trapezoidal rule approximations for 5s, 10s, and 30s so the interval, h,% is halved in each.
%% One Segment Trapezoidal Rule% Integral after 5s (approximate height)T5_1 = (b-a)*((v(a)-v(b))/2);% ErrorEt5_1 = abs(A5-T5_1); % true erroret5_1 = (Et5_1/A5)*100; % percent relative error
% Integral after 10s (approximate height)T10_1 = (c-a)*((v(a)-v(c))/2);% ErrorEt10_1 = abs(A10-T10_1); % true erroret10_1 = (Et10_1/A10)*100; % percent relative error
% Integral after 30s (approximate height)T30_1 = (d-a)*((v(a)-v(d))/2);% ErrorEt30_1 = abs(A30-T30_1); % true erroret30_1 = (Et30_1/A30)*100; % percent relative error
%% Two Segment Trapezoidal Rulen = 2; % number of segments
% Estimate height after 5 secondsh5_2 = (b-a)/n; % width of trapezoid% Integral after 5s (approximate height)T5_2 = (h5_2/2)*(v(a)+2*v(a+h5_2)+v(b)); % ErrorEt5_2 = abs(A5-T5_2); % true erroret5_2 = (Et5_2/A5)*100; % percent relative error
% Estimate height after 10 secondsh10_2 = (c-a)/n; % width of trapezoid% Integral after 10s (approximate height)T10_2 = (h10_2/2)*(v(a)+2*v(a+h10_2)+v(c));% ErrorEt10_2 = abs(A10-T10_2); % true erroret10_2 = (Et10_2/A10)*100; % percent relative error
% Estimate height after 30 secondsh30_2 = (d-a)/n; % width of trapezoid% Integral after 30s (approximate height)T30_2 = (h30_2/2)*(v(a)+2*v(a+h30_2)+v(d));% ErrorEt30_2 = abs(A30-T30_2); % true erroret30_2 = (Et30_2/A30)*100; % percent relative error
%% Four Segment Trapezoidal Rulen = 4; % number of segments
% Estimate height after 5 secondsh5_4 = (b-a)/n; % width of trapezoid% Integral after 5s (approximate height)T5_4 = (h5_4/2)*(v(a)+2*(v(a+h5_4)+v(a+2*h5_4)+v(a+3*h5_4))+v(b)); % ErrorEt5_4 = abs(A5-T5_4); % true erroret5_4 = (Et5_4/A5)*100; % percent relative error
% Estimate height after 10 secondsh10_4 = (c-a)/n; % width of trapezoid% Integral after 10s (approximate height)T10_4 = (h10_4/2)*(v(a)+2*(v(a+h10_4)+v(a+2*h10_4)+v(a+3*h10_4))+v(c)); % ErrorEt10_4 = abs(A10-T10_4); % true erroret10_4 = (Et10_4/A10)*100; % percent relative error
% Estimate height after 30 secondsh30_4 = (d-a)/n; % width of trapezoid% Integral after 30s (approximate height)T30_4 = (h30_4/2)*(v(a)+2*(v(a+h30_4)+v(a+2*h30_4)+v(a+3*h30_4))+v(d)); % ErrorEt30_4 = abs(A30-T30_4); % true erroret30_4 = (Et30_4/A30)*100; % percent relative error
%% Eight Segment Trapezoidal Rulen = 8; % number of segments
% Estimate height after 5 secondsh5_8 = (b-a)/n; % width of trapezoid% Integral after 5s (approximate height)T5_8 = (h5_8/2)*(v(a)+2*(v(a+h5_8)+v(a+2*h5_8)+v(a+3*h5_8)+v(a+4*h5_8)+v(a+5*h5_8)+v(a+6*h5_8)+v(a+7*h5_8))+v(b)); % ErrorEt5_8 = abs(A5-T5_8); % true erroret5_8 = (Et5_8/A5)*100; % percent relative error
% Estimate height after 10 secondsh10_8 = (c-a)/n; % width of trapezoid% Integral after 10s (approximate height)T10_8 = (h10_8/2)*(v(a)+2*(v(a+h10_8)+v(a+2*h10_8)+v(a+3*h10_8)+v(a+4*h10_8)+v(a+5*h10_8)+v(a+6*h10_8)+v(a+7*h10_8))+v(c)); % ErrorEt10_8 = abs(A10-T10_8); % true erroret10_8 = (Et10_8/A10)*100; % percent relative error
% Estimate height after 30 secondsh30_8 = (d-a)/n; % width of trapezoid% Integral after 30s (approximate height)T30_8 = (h30_8/2)*(v(a)+2*(v(a+h30_8)+v(a+2*h30_8)+v(a+3*h30_8)+v(a+4*h30_8)+v(a+5*h30_8)+v(a+6*h30_8)+v(a+7*h30_8))+v(d)); % ErrorEt30_8 = abs(A30-T30_8); % true erroret30_8 = (Et30_8/A30)*100; % percent relative error
%% O(h^4) Romberg Methods% 5sI4_5a = (4/3)*T5_2-(1/3)*T5_1;% ErrorE5h4a = abs(A5-I4_5a); % true errore5h4a = (E5h4a/A5)*100; % percent relative error
I4_5b = (4/3)*T5_4-(1/3)*T5_2;% ErrorE5h4b = abs(A5-I4_5b); % true errore5h4b = (E5h4b/A5)*100; % percent relative error
I4_5c = (4/3)*T5_8-(1/3)*T5_4;% ErrorE5h4c = abs(A5-I4_5c); % true errore5h4c = (E5h4c/A5)*100; % percent relative error
% 10sI4_10a = (4/3)*T10_2-(1/3)*T10_1;% ErrorE10h4a = abs(A10-I4_10a); % true errore10h4a = (E10h4a/A10)*100; % percent relative error
I4_10b = (4/3)*T10_4-(1/3)*T10_2;% ErrorE10h4b = abs(A10-I4_10b); % true errore10h4b = (E10h4b/A10)*100; % percent relative error
I4_10c = (4/3)*T10_8-(1/3)*T10_4;% ErrorE10h4c = abs(A10-I4_10c); % true errore10h4c = (E10h4c/A10)*100; % percent relative error
% 30sI4_30a = (4/3)*T30_2-(1/3)*T30_1;% ErrorE30h4a = abs(A30-I4_30a); % true errore30h4a = (E30h4a/A30)*100; % percent relative error
I4_30b = (4/3)*T30_4-(1/3)*T30_2;% ErrorE30h4b = abs(A30-I4_30b); % true errore30h4b = (E30h4b/A30)*100; % percent relative error
I4_30c = (4/3)*T30_8-(1/3)*T30_4;% ErrorE30h4c = abs(A30-I4_30c); % true errore30h4c = (E30h4c/A30)*100; % percent relative error
%% O(h^6) Romberg Methods% 5sI6_5a = (16/5)*I4_5b-(1/15)*I4_5a;% ErrorE5h6a = abs(A5-I6_5a); % true errore5h6a = (E5h6a/A5)*100; % percent relative error
I6_5b = (16/5)*I4_5c-(1/15)*I4_5b;% ErrorE5h6b = abs(A5-I6_5b); % true errore5h6b = (E5h6b/A5)*100; % percent relative error
% 10sI6_10a = (16/5)*I4_10b-(1/15)*I4_10a;% ErrorE10h6a = abs(A10-I6_10a); % true errore10h6a = (E10h6a/A10)*100; % percent relative error
I6_10b = (16/5)*I4_10c-(1/15)*I4_10b;% ErrorE10h6b = abs(A10-I6_10b); % true errore10h6b = (E10h6b/A10)*100; % percent relative error
% 30sI6_30a = (16/5)*I4_30b-(1/15)*I4_30a;% ErrorE30h6a = abs(A30-I6_30a); % true errore30h6a = (E30h6a/A30)*100; % percent relative error
I6_30b = (16/5)*I4_30c-(1/15)*I4_30b;% ErrorE30h4b = abs(A30-I6_30b); % true errore30h4b = (E30h4b/A30)*100; % percent relative error
%% O(h^8) Romberg Methods% 5sI8_5 = (64/63)*I6_5b-(1/63)*I6_5a;% ErrorE5h8 = abs(A5-I8_5); % true errore5h8 = (E5h8/A5)*100; % percent relative error
% 10sI8_10 = (64/63)*I6_10b-(1/63)*I6_10a;% ErrorE10h8 = abs(A10-I8_10); % true errore10h8 = (E10h8/A10)*100; % percent relative error
% 30sI8_30 = (64/63)*I6_30b-(1/63)*I6_30a;% ErrorE30h8 = abs(A30-I8_30); % true errore30h8 = (E30h8/A30)*100; % percent relative error
%% Part B - Numerical Differentiation% Analytical Solutionta = -5:5:35; % time intervalaccel = diff(v(ta));t=0:5:35;% Plot acceleration over timefigure(2)plot(t,accel,'r-','LineWidth',2)grid on; box on;title('Acceleration of the Rocket over Time','FontSize',16)xlabel('Time (s)','FontSize',14)ylabel('Acceleration (m/s^2)','FontSize',14)
%% Numerical Solution - Finite Divided Differencest1 = 0; t2 = 5; t3 = 10; t4 = 15; t5 = 20; t6 = 25; t7 = 30;
% Forward finite difference approximationf0 = (v(t2)-v(t1))/(t2-t1);
% Centered finite difference approximationf5 = (v(t3)-v(t1))/(t3-t1);f10 = (v(t4)-v(t2))/(t4-t2);f15 = (v(t5)-v(t3))/(t5-t3);f20 = (v(t6)-v(t4))/(t6-t4);f25 = (v(t7)-v(t5))/(t7-t5);
% Backward finite difference approximationf30 = (v(t7)-v(t6))/(t7-t6);
num_diff = [f0 f5 f10 f15 f20 f25 f30];
%% Error between Numerical and AnalyticalE0 = abs((accel(1)-f0)/accel(1))*100;E5 = abs((accel(2)-f5)/accel(2))*100;E10 = abs((accel(3)-f10)/accel(3))*100;E15 = abs((accel(4)-f15)/accel(4))*100;E20 = abs((accel(5)-f20)/accel(5))*100;E25 = abs((accel(6)-f25)/accel(6))*100;E30 = abs((accel(6)-f30)/accel(7))*100;E = [E0 E5 E10 E15 E20 E25 E30];

Page updated

Report abuse