Applied Numerical Methods Using MATLAB. Won Y. Yang
can sample every ( kM − m)th element of a given row vector sequence x to make a new sequence y. Write a MATLAB statement to use the function for sampling every (3k − 2)th element of x = [1 3 7 2 4 9] to getfunction y=<b>sampling</b><![CDATA[(x,M,m) m=mod(m,M); Nx=numel(x); N=floor(Nx/M); y=x([1:N]*?-?); if Nx-N*M>=M-m, y=[y x(N*M+M-m)]; end Complete the following function ‘ rotate_r(x,M)’ so that it can rotate a given row vector sequence x right by M samples to make a new sequence y. Write a MATLAB statement to use the function for rotating x = [1 2 3 4 5] to getfunction xl=<b>rotate_r</b><![CDATA[(x,M) N=size(x,2); M=mod(M,N); xl=[x(:,end-?+1:end) x(:,1:end-?)];
18 1.18 Distribution of a Random Variable – HistogramComplete the following function ‘ randu(N,a,b)’ that uses the MATLAB function ‘ rand()’ to generate an N‐dimensional random vector having the uniform distribution over [ a, b] and draws the histogram (with 20 bins) for the distribution of the elements of the generated vector as Figure 1.5. Then, find the average height of the histogram that you can get by running the following statement: >randu(1000,-2,2)function x=<b>randu</b><![CDATA[(N,a,b) % generates an N-dimensional random vector with U(a,b) x=a+(?-?)*rand(1,N); if nargout==0, hist(x(:),20); end
19 1.19 Number RepresentationIn Section 1.2.1, we looked over how a number is represented in 64 bits. For example, the IEEE 64‐bit floating‐point number system represents the number 3(21 ≤ 3 < 22) belonging to the range R1 = [21,22) with E = 1 as010000000000100000000000...........00000000000000000000400800........00000where the exponent and the mantissa areThis can be confirmed by typing the following statement into MATLAB command window: >fprintf('3=%bx\n',3) or >format hex, 3, format short which will print out onto the screen 4008000000000000This is exactly the hexadecimal representation of the number 3 as we expected. Find the IEEE 64‐bit floating‐point number representation of the number 14 and use the command ‘fprintf()’ to check if the result is right.(cf) Since the INTEL system stores numbers in little‐endian (byte order) format with more significant bytes in the memory of higher address number, you might see 0000000000000840in the old versions of MATLAB, which is reversely ordered in the unit of byte (8 bits = 2 hexadecimal digits) so that the number is represented with the most/least significant byte on the right/left side.
20 1.20 Resolution of Number Representation and Quantization ErrorIn Section 1.2.1, we have seen that adding 2−22 to 230 makes some difference, while adding 2−23 to 230 makes no difference due to the bit shift by over 52 bits for alignment before addition. How about subtracting 2−23 from 230? In contrast with the addition of 2−23 to 230, it makes difference as you can see by running the following MATLAB statement: >x=2̂30; x+2̂-23==x, x-2̂-23==x which will give you the logical answer 1 (true) and 0 (false). Justify this result based on the difference of resolution of two ranges [230,231) and [229,230) to which the true values of computational results (230 + 2−23) and (230 − 2−23) belong, respectively. Note from Eq. (1.2.5) that the resolutions, i.e. the maximum quantization errors are ΔE = 2E−52 = 2−52+30 = 2−22 and 2−52+29 = 2−23, respectively. For details, refer to Figure P1.20, which illustrates the process of addition/subtraction with 4 mantissa bits, 1 hidden bit, and 1 guard bit.Figure P1.20 Process of addition/subtraction with four mantissa bits.
21 1.21 Resolution of Number Representation and Quantization ErrorWhat is the result of running the following statement? >7/100*100-7How do you compare the absolute value of this answer with the resolution Δ of the range to which 7 belongs?Find how many numbers are susceptible to this kind of quantization error caused by division/ multiplication by 100, among the numbers from 1 to 31.What will be the result of running the following script? Why?%nm01p21.m: Quantization Error x=2-2̂-50; for n=1:2̂3 x=x+2̂-52; fprintf('%20.18E\n',x) end
22 1.22 Avoiding Large Errors/Overflow/UnderflowFor x = 9.8201 and y = 10.2199, evaluate the following two expressions that are mathematically equivalent and tell which is better in terms of the power of resisting the overflow.(P1.22.1a) (P1.22.1b) Also for x = 9.8−201 and y = 10.2−199, evaluate the above two expressions and tell which is better in terms of the power of resisting the underflow.With a = c = 1 and for 100 values of b over the interval [107.4, 108.5] generated by the MATLAB command ‘ logspace(7.4,8.5,100)’, evaluate the following two formulas (for the roots of a quadratic equation) that are mathematically equivalent and plot the values of the second root of each pair. Noting that the true values are not available and so the shape of solution graph is only one practical basis on which we can assess the quality of numerical solutions, tell which is better in resisting the loss of significance.(P1.22.2a) (P1.22.2b) For 100 values of x over the interval [10−9, 10−7.4], evaluate the following two expressions that are mathematically equivalent, plot them, and based on the graphs, tell which is better in terms of resisting the loss of significance.(P1.22.3a) (P1.22.3b) For 100 values of x over the interval [1014, 1016], evaluate the following two expressions that are mathematically equivalent, plot them, and based on the graphs, tell which is better in terms of resisting the loss of significance.(P1.22.4a) (P1.22.4b) On purpose to find the value of (300125/125!)e−300, run the following statement:>300̂125/prod([1:125])*exp(-300)What is the result? Is it of any help to change the order of multiplication/division? As an alternative, referring to the script “nm131_2_1.m” (Section 1.3.1), complete the following function ‘ Poisson_pdf()’ so that it can compute(P1.22.5) with nested computing (in a recursive way), say, like p(k + 1) = p(k) * λ/k and then, use the function to find the value of (300125/125!)e−300.function p=<b>Poisson_pdf</b><![CDATA[(K,lambda) p=exp(-??????); for k=1:?, p= ?*lambda/?; endMake a function that computes the sum(P1.22.6) and then uses the function to find the value of S(155).function S=<b>Poisson_pdf_sum</b><![CDATA[(K,lambda) p=exp(-lambda); S=0; for k=1:K, p= ?*lambda/k; S=?+p; end
23 1.23 Nested Computing for Computational Efficiency(a) The Hermite polynomial [K−2, W−3]Consider the Hermite polynomial defined as(P1.23.1) (i) Show that the derivative of this polynomial function can be written as(P1.23.2) whence the (N + 1)th‐degree Hermite polynomial can be obtained recursively from the Nth‐degree Hermite polynomial as(P1.23.3) (ii) Complete the following MATLAB function ‘ Hermitp(N)’ so that it can use Eq. (P1.23.3) to generate the Nth‐degree Hermite polynomial HN(x).function p=<b>Hermitp</b><![CDATA[(N) % Hn+1(x)=2xHn(x)-Hn'(x) if N<=0, p=1; else p=[2 ?]; for n=2:N, p= 2*[p ?]-[0 ? polyder(p)]; end end(b) The Bessel function of the first kind [K‐2, W‐1]Consider the Bessel function of the first kind of order k defined as(P1.23.4a) (P1.23.4b) function [J,JJ]=<b>Jkb</b><![CDATA[(K,beta) %first kind of kth-order Bessel ftn tmpk= ones(size(beta)); for k=0:K tmp= tmpk; JJ(k+1,:)=tmp; for m=1:100 tmp= -tmp.*????.̂2/4/m/(m+?); JJ(k+1,:)= JJ(k+1,:)?tmp; if norm(tmp)<.001, break; end end tmpk=tmpk.*beta/2/(k+1); end J=JJ(K+1,:);function y=<b>Bessel_integrand</b><![CDATA[(x,beta,k) y=cos(k*x-beta*sin(x));%nm01p23b.m: Bessel_ftn beta=0:.05:15; K=15; tic for i=1:length(beta) % Integration J15_qd(i)=quad('Bessel_integrand',0,pi,[],0,beta(i),K)/pi; end time_qd=toc tic, J15_Jkb=Jkb(K,beta); time_Jkb=toc/K % Nested Computing tic, J15_bj=besselj(K,beta); time_bj=toc discrepancy1 = norm(J15_qd-J15_bj) discrepancy2 = norm(J15_Jkb-J15_bj)(i) Define the integrand of (P1.23.4a) in the name of ‘ Bessel_integrand (x,beta,k)’ and save it in an M‐file named “Bessel_integrand.m”.(ii) Complete the above function ‘ Jkb(K,beta)’ so that it can use (P1.23.4b) in a recursive way to compute Jk(β) of order k = 1:K for given K and β ( beta).(iii) Run the above script “nm01p21b3.m”, which computes J15(β) for β = 0 : 0.05 : 15 in three ways, that is, using Eq. (P1.23.4a), Eq. (P1.23.4b) (cast into ‘ Jkb()’), and the MATLAB built‐in function ‘ besselj()’. Do they conform with each other?(cf) Note that ‘ Jkb(K,beta)’ computes Jk(β) of order k = 1:K, while the integration does for only k = K.(cf) Note also that the MATLAB built‐in function