Solution 2 27/08/02
Solution 2: Numerical dierentiation and integration
Problem 1
As the space shuttle reenters the atmosphere, the astronauts undergo extreme accelerations as the shuttle decelerates from its orbital velocity of 7620 km/s to its landing velocity of 352 km/hr over the course of roughly 100 minutes. You are given a data set of position points equally spaced in time of a test shuttle as it reenters the atmosphere, as shown in Figure 1, and are assigned the task of determining whether or not the test shuttle decelerates at
Figure 1: Position as a function of time for the test shuttle reentry prole more than the unconsciousness limit of gu = 6.3 Gs, or 6.3 times gravitational acceleration of 9.8m/s2 . (source: [Link] Accelerations higher than the unconsciousness limit will cause the astronauts to faint, and the test shuttle will have to be redesigned, wasting billions of Rand. The 100-point data set can be found at [Link] and consists of two columns, the rst being the time in seconds, and the second being the position in meters. If s is position, then the velocity is given by v = ds/dt and the acceleration is given by a = dv/dt. You decide to compute the velocity at points in time halfway between the data points n with the second order accurate approximation vn+1/2 sn+1 sn t n {1, 2, . . . , N 1} , (1)
Solution 2 27/08/02 and the acceleration at the data points with the second order accurate approximation an vn+1/2 vn1/2 t n {2, . . . , N 1} ,
(2)
where N is the number of data points. 1. Plot the deceleration of the shuttle normalized by the unconsciousness limit (an /gu ) as a function of time. Does the shuttle need to be redesigned? Figure 2 depicts the normalized shuttle deceleration as a function of time. Because the solid line does not cross the critical consciousness line, the astronauts will not faint, and the shuttle does not need to be redesigned. 2. If partial loss of vision occurs at gv = 3.4 Gs, how long would the astronauts estimate to experience vision problems if they were in the test shuttle (round this answer to the nearest minute)? Figure 2 shows that the astronauts may experience vision problems, since the acceleration exceeds the critical vision limit for about 15 minutes between 22 minutes and 37 minutes.
Figure 2: Normalized shuttle deceleration, showing how the astronauts will not faint and thus the shuttle does not need to be redesigned. The astronauts may experience vision problems, however, because the shuttle deceleration exceeds the critical deceleration for vision problems.
Solution 2 27/08/02
Problem 2
The error function is a good example of an integral that cannot be obtained analytically, and is given by y 2 2 Erf(y) = ex dx . (3) 0 1. Write two functions, trapz and simp13, each of which compute the integral of an arbitrary function f on an equispaced grid x using the trapezoidal rule and Simpsons 1/3 rule, so that in Octave, the functions are used as >> I1 = trapz(x,f); >> I2 = simp13(x,f); The two functions can be downloaded from [Link] [Link] and are given below Trapezoidal rule code (trapz.m) % % Introduction to numerical methods for finance students % Assignment 2 Problem 2: trapz function % 27/08/02 % Oliver Fringer % % This m-file computes the integral of a function f over the % equally spaced points x using the trapezoidal rule. % function I = trapz(x,f) h = N = I = for I end I = x(2)-x(1); length(f); 0; n=2:N-1, = I + h*f(n); I + h/2*(f(1) + f(N));
endfunction
Solution 2 27/08/02 Simpsons 1/3 rule code (simp13.m) % % Introduction to numerical methods for finance students % Assignment 2 Problem 2: simp13 function % 27/08/02 % Oliver Fringer % % This m-file computes the integral of a function f over the % equally spaced points x using Simpsons 1/3 rule. The number % of intervals must be even! % function I = simp13(x,f) h = x(2)-x(1); I = 0; N = length(f); if(rem(length(f)-1,2)~=0) error(Simpsons 1/3 rule must use an even number of elements); end for n=1:length(f)/2 I = I + 4/3*h*f(2*n); end for n=1:length(f)/2-1 I = I + 2/3*h*f(2*n+1); end I = I + h/3*(f(1)+f(N)); endfunction
2. Using these functions, compute the value of Erf(1) and show that the trapezoidal rule is second order accurate with respect to the grid spacing h and that Simpsons 1/3 rule is fourth order accurate. Use N = 20, 40, 80, 160, 320, and 640 intervals to discretize x with xi = (i 1)h, i {1, . . . , N + 1} and h = 1/N , and plot the error as a function of h for each method. The error is given by Error = ErfN (1) Erfexact (1) , Erfexact (1) (4)
where ErfN (1) is the estimate you get using one of the methods with N intervals, and Erfexact (1) is the exact value given to you by the erf function in Octave. The code for this problem can be downloaded from [Link]
Solution 2 27/08/02
and the convergence results are shown in Figure 3. Because the slope of the trapezoidal rule error follows the h2 line, this shows that the trapezoidal rule is second order accurate, and because the slope of the Simpsons 1/3 rule error follows the h4 line, Simpsons 1/3 rule is fourth order accurate.
Figure 3: Convergence results for the Trapezoidal and Simpsons 1/3 rule, showing that the trapezoidal rule is second order accurate and Simpsons 1/3 rule is fourth order accurate