Читать книгу Applied Numerical Methods Using MATLAB - Won Y. Yang - Страница 42

1 1.1 Creating a Data File and Retrieving/Plotting Data Saved in a Data FileUsing the MATLAB editor, make a script “nm01p01a.m” that lets its user input data pairs of heights (ft) and weights (lb) of as many persons as he wants till he presses <Enter> and save the whole data in the form of an N × 2 matrix into an ASCII data file ‘***.txt’ named by the user. If you have no idea how to compose such a script, you can permutate the statements in the box beneath to make your script. Store the script in the file named “nm01p01a.m” and run it to save the following data into the data file named ‘hw.txt’:%nm01p01a.m: input data pairs and save them into an ASCII data file k=0; while 1 end k=k+1; x(k,1)=h; h=input('Enter height:') x(k,2)=input('Enter weight:') if isempty(h), break; end cd('c:\MATLAB\nma') % to change current working directory filename=input('Enter filename(.txt):','s'); filename=[filename '.txt']; % String concatenation save(filename,'x','-ascii')Make a MATLAB script “nm01p01b.m” that reads (loads) the data file ‘hw.txt’ made in (a), plots the data as in Figure 1.1a in the upper‐left region of the screen divided into four regions like Figure 1.5 and plots the data in the form of PWL graph describing the relationship between the height%nm01p01b.m: to read the data file and plot the data cd('c:\MATLAB\nma') % to change current working directory weight=hw(I,2); load hw.txt clf subplot(221) plot(hw) subplot(222) axis([5 7 160 200]) plot(height,weight,'-+') [height,I]=sort(hw(:,1)); and the weight in the upper‐right region of the screen. Let each data pair be denoted by the symbol ‘+’ on the graph. Also let the ranges of height and weight be [5,7] and [160,200], respectively. If you have no idea, you can permutate the statements in the below box. Additionally, run the script to check if it works fine.

2 1.2 Text Printout of Alphanumeric DataMake a MATLAB function ‘ max_array(A)’ that uses the ‘ max()’ command to find one of the maximum elements of a matrix A given as its input argument and uses the ‘ fprintf()’ command to print it onto the screen together with its row/column indices in the following format. '\n Max(A) is A(%2d,%2d)=%5.2f\n',row_index,col_index,maxAAdditionally, try it to have the maximum element of an arbitrary matrix (generated by the following two consecutive commands) printed in this format onto the screen. >rand('state',sum(100*clock)), rand(3)

3 1.3 Plotting the Mesh Graph of a Two‐dimensional FunctionConsider the MATLAB script “nm01p03a.m”, whose objective is to draw a cone.The statement on the sixth line seems to be dispensable. Run the script with and without this line and see what happens.If you want to plot the function ‘ fcone(x,y)’ defined in another M‐file “fcone.m”, as an inline function, or as a function handle, how will you modify this script?If you replace the fifth line by ‘ Z=1‐abs(X)‐abs(Y);’, what difference does it make?%nm01p03a.m: to plot a cone x=-1:0.02:1; y=-1:0.02:1; [X,Y]=meshgrid(x,y); Z=1-sqrt(X.̂2+Y.̂2); Z=max(Z,zeros(size(Z))); mesh(X,Y,Z)function z=<b>fcone</b><![CDATA[(x,y) z=1-sqrt(x.̂2+y.̂2);

4 1.4 Plotting the Mesh Graph of Stratigraphic StructureTable P1.4 The depth of the rock layer.y‐Coordinatex‐Coordinate0.11.22.53.64.80.54103903804204501.43953754104354552.23654054304554703.53704004204454354.6385395410395410Consider the incomplete MATLAB script “nm01p04.m”, whose objective is to draw a stratigraphic structure of the area around Pennsylvania State University from the several perspective point of view. The data about the depth of the rock layer at 5 × 5 sites are listed in Table P1.4. Supplement the incomplete parts of the script so that it serves the purpose and run the script to answer the following questions. If you complete it properly and run it, MATLAB will show you the four similar graphs at the four corners of the screen and be waiting for you to press any key.At what value of k does MATLAB show you the mesh/surface‐type graphs that are the most similar to the first graphs? From this result, what do you guess are the default values of the azimuth or horizontal rotation angle and the vertical elevation angle (in degrees) of the perspective view point?As the first input argument Az of the command ‘ view(Az,E1)’ decreases, in which direction does the perspective view point revolve round the z‐axis, clockwise or counterclockwise (seen from the above)?As the second input argument El of the command ‘ view(Az,E1)’ increases, does the perspective view point move up or down along the z‐axis?What is the difference between the plotting commands ‘ mesh()’ and ‘ meshc()’?What is the difference between the usages of the command ‘ view()’ with two input arguments Az, El, and with a three‐dimensional vector argument [x,y,z]?%nm01p04.m: to plot a stratigraphic structure clear, clf x=[0.1 .. .. . ]; y=[0.5 .. .. . ]; Z=[410 390 .. .. .. .. ]; [X,Y]=meshgrid(x,y); subplot(221), mesh(X,Y,500-Z) subplot(222), surf(X,Y,500-Z) subplot(223), meshc(X,Y,500-Z) subplot(224), meshz(X,Y,500-Z) pause for k=0:7 Az=-12.5*k; El=10*k; Azr=Az*pi/180; Elr=El*pi/180; subplot(221), view(Az,El) subplot(222), k, view([sin(Azr),-cos(Azr),tan(Elr)]), pause %pause(1) end

5 1.5 Plotting the Graph of a Function over an Interval Containing Its Singular PointNoting that the tangent function f(x) = tan(x) is singular at x = π/2, 3π/2, …, let us plot its graph over [0, 2π] as follows:Define the domain vector x consisting of sufficiently many intermediate point xis along the x‐axis and the corresponding vector y consisting of the function values at xis and plot the vector y over the vector x. You may use the following statements.>x=[0:0.01:2*pi]; y=tan(x); >subplot(221), plot(x,y); xlim([0 2*pi])Which one is the most similar to what you have got, among the graphs shown in Figure P1.5? Is it far from your expectation?Expecting to get the better graph, we scale it up along the y‐axis by using the following command.>axis([0 2*pi -10 10])Which one is the most similar to what you have got, among the graphs shown in Figure P1.5? Is it closer to your expectation than what you got in (a)?Most probably, you must be nervous about the straight‐lines at the singular points x = π/2 and 3π/2. The severer you get disturbed by the lines that must not be there, the better you are at the numerical stuffs. As an alternative to avoid such a singular happening, you can try dividing the interval into three sections excluding the two singular points as follows:x1=[0:0.01:pi/2-0.01]; x2=[pi/2+0.01:0.01:3*pi/2-0.01]; x3=[3*pi/2+0.01:0.01:2*pi]; y1=tan(x1); y2=tan(x2); y3=tan(x3); subplot(222), plot(x1,y1,x2,y2,x3,y3), axis([0 2*pi -10 10])How about trying the easy plotting command ‘ ezplot()’ or ‘ fplot()’? Do they work properly?>ezplot('tan(x)',[0 2*pi]) >fplot(@(x)tan(x),[0 2*pi]) % plot an anonymous functionFigure P1.5 Plotting the graph of f(x) = tan x.

6 1.6 Plotting the Graph of a Sinc FunctionThe sinc function is defined as(P1.6.1) whose value at x = 0 is(P1.6.2) We are going to plot the graph of this function over [ ‐4π,4π].(a) Casually, you may try as follows:>x=[-100:100]*pi/25; y=sin(x)./x; >plot(x,y), axis([-15 15 -0.4 1.2]) Even if no warning message shows up, is there anything odd about the graph?(b) How about trying with a different domain vector?>x=[-4*pi:0.1:+4*pi]; y=sin(x)./x; >plot(x,y), axis([-15 15 -0.4 1.2]) Surprisingly, MATLAB gives us the function values without any complaint and presents a nice graph of the sinc function. What is the difference between (a) and (b)?(cf) Actually, we would have no problem if we used the MATLAB built‐in function ‘ sinc()’.

7 1.7 Term‐Wise (Element‐by‐Element) Operation in MATLAB User FunctionsLet the function f1(x) be defined without one or both of the dot ( .) operators in Section 1.1.6. Could we still get the output vector consisting of the function values for the several values in the input vector? You can run the following statements and see the results.>f1=@(x)1./(1+8*x̂2); f1([0 1]) >f1=@(x)1/(1+8*x.̂2); f1([0 1]) Let the function f1(x) be defined with both of the dot ( .) operators as in Section 1.1.6. What would we get by running the following statements?>f1=@(x)1./(1+8*x.̂2); f1([0 1]')

8 1.8 In‐line Function and M‐file Function with the Integral Function ‘ quad()’As will be seen in Section 5.8, one of the MATLAB built‐in functions for computing numerical integrals is ‘ quad()’, the usual usage of which is(P1.8.1) wheref: the name of the integrand function (M‐file name should be categorized by ' ' )a,b: the lower/upper bound of the integration intervaltol: the error tolerance (10−6 by default []), trace: 1(on)/0(off) (0 by default [])p1,p2, ..: additional parameters to be passed directly to function fLet us use this ‘ quad()’ function with an in‐line function and an M‐file function to obtain(P1.8.2a) and(P1.8.2b) where(P1.8.3) Complete and run the following script “nm01p08.m” to compute the two integrals (P1.8.2a and P1.8.2b).%nm01p08.m m=1; sigma=2; x0=1; Gpdf='exp(-(x-?).̂2/2/?????̂2)/sqrt(2*pi)/sigma'; xGpdf=inline(['(x-x0).*' Gpdf],'x','m','?????','x0'); x2Gpdf=inline(['(x-x0).̂2.*' Gpdf],'x','?','sigma','x0'); int_xGpdf=quad(xGpdf,m-10,m+10,[],0,?,sigma,x0) int_x2Gpdf=quad(x2Gpdf,m-10,m+10,[],0,m,?????,x0)

9 1.9 μ‐Law Function Defined in an M‐FileThe so‐called μ‐law function and μ−1‐law function used for uniform quantization is defined as(P1.9.1a) (P1.9.1b) Below are the incomplete μ‐law function ‘ mulaw()’, μ−1‐law function ‘ mulaw_inv()’, and a script “nm01p09.m” that are all supposed to be saved as M‐files. Complete them and run the script to do the following jobs with μ = 10, 50, and 255:– Find the values y of the μ‐law function for x=[‐1:0.005:1] and plot the graph of y vs. x.– Find the values x0 of the μ−1‐law function for y.– Compute the discrepancy between x and x0.function [y,xmax]=<b>mulaw</b><![CDATA[(x,mu,ymax) if nargin<3, ymax=1; end xmax=max(abs(x)); y=ymax*log(1+mu*???(x/xmax))./log(1+mu).*????(x); % Eq.(P1.9a)function x=<b>mulaw_inv</b><![CDATA[(y,mu,xmax) % Inverse of mu-law if nargin<3, xmax=1; ymax=1; else ymax=max(abs(y)); end if nargin<2, mu=255; end x=xmax.*(((1+??).̂(abs(?)/y???)-1)/??).*sign(y); % Eq.(P1.9b)%nm01p09.m: to plot the mulaw curve x=[-1:.005:1]; mu=[10 50 255]; for i=1:3 [y,xmax]=mulaw(x,mu(i),1); x0=mulaw_inv(y,mu(i),xmax); discrepancy=norm(x-x0) plot(x,y,'b-', x,x0,'r-'), hold on end

10 1.10 Analog‐digital converter (ADC)Below are two ADC routines ‘ adc1(a,b,c)’ and ‘ adc2(a,b,c)’, which assign the corresponding digital value c(i) to each one of the analog data belonging to the quantization interval [ b(i), b(i+1)]. Let the boundary vector and the centroid (level) vector be, respectively, b=[-3 -2 -1 0 1 2 3]; c=[-2.5 -1.5 -0.5 0.5 1.5 2.5];Note that the digital value corresponding to an analog data smaller than b(1) (the smallest boundary) should be c(1) (the smallest centroid) and that corresponding to an analog data larger than b(N) (the largest boundary) should be c(N) (the largest centroid) where the boundary and centroid vectors are assumed to be arranged in ascending order and N is the number of the centroids (quantization levels) or quantization intervals.Explain the jobs of lines 8, 9, and 10 of the following ADC function ‘ adc1()’.Explain the jobs of lines 3, 5, and 7 of the following ADC function ‘ adc2()’.Run the script “nm01p10.m” (see below box) to get two graphs as shown in Figure P1.10. Explain about the graphs.Figure P1.10 The characteristic of an analog‐to‐digital converter (ADC).function d=<b>adc1</b><![CDATA[(a,b,c) %Analog-to-Digital Converter %Input a=analog signal, b(1:N+1)=boundary vector % c(1:N)=centroid vector %Output: d=digital samples N=length(c); for n=1:length(a) I=find(a(n)<b(2:N)); if ∼isempty(I), d(n)=c(I(1)); else d(n)=c(N); end endfunction d=<b>adc2</b><![CDATA[(a,b,c) N=length(c); d(find(a<b(2)))=c(1); for i=2:N-1 index=find(b(i)<=a&a<b(i+1)); d(index)=c(i); end d(find(b(N)<=a))=c(N);%nm01p10.m b=[-3 -2 -1 0 1 2 3]; % Boundary vector c=[-2.5 -1.5 -0.5 0.5 1.5 2.5]; % Centroid vector % Plot the input-output relationship xa=[-300:300]/100; % Analog data in the range [-3∼+3] xd1=adc1(xa,b,c); xd2=adc2(xa,b,c); subplot(221) plot(xa,xd1, xa,xd2,'r:') % Output of ADC to a sinusoidal input t=[0:200]/100*pi; xa=3*sin(t); xd1=adc1(xa,b,c); %xd2=adc2(xa,b,c); subplot(222) plot(t,xa, t,xd1,'r')

11 1.11 Decimal‐to‐Binary/Octal ConversionsReferring to the decimal‐to‐binary conversion process illustrated in Figure 1.6, complete the following function ‘ dec2bin_my()’ so that it can perform the conversion process. Then, using the function, convert two decimal numbers 3 and 14 to their corresponding binary numbers and compare the results with those obtained using the MATLAB built‐in function ‘ dec2bin()’. Can you modify the function so that the binary numbers corresponding to each of x can be placed columnwise, that is, from top to bottom, rather than from left to right when a vector consisting of multiple decimal numbers is given as the first input argument x?Complete the following function ‘ bin2dec_my()’ so that it can perform the binary‐to‐digital conversion process:function y=<b>dec2bin_my</b><![CDATA[(x) % converts given decimal numbers into binary numbers of N bits y=[]; xmax=max(x); N=1; while xmax>1 xmax=floor(xmax/2); N=N+1; end for n=1:length(x) xn=x(n); for i=?:-1:1; xn_1= floor(xn/?); yn(i)= xn-?*xn_1; xn=xn_1; end y= [y yn]; endfunction y=<b>bin2dec_my</b><![CDATA[(x) [M,N]=size(x); % Number of binary numbers and Number of bits for m=1:M y(m,:)=x(m,:)*?.̂[?-1:-1:0]'; endThen, using the function, convert two binary numbers [0 0 1 1] and [1 1 1 0] to their corresponding digital numbers and compare the results with those obtained using the MATLAB built‐in function ‘ bin2dec()’, by typing the following statements at MATLAB prompt:bin2dec_my([0 0 1 1; 1 1 1 0]) bin2dec(['0011'; '1110'])(c) Modify the decimal‐to‐binary converting function ‘ dec2bin_my()’ into another function ‘ dec2oct_my()’ so that it can perform the decimal‐to‐octa conversion process. Then, using the function, convert two decimal numbers 14 and 3 to their corresponding octal numbers and compare the results with those obtained using the MATLAB built‐in function ‘ deci2oct()’.

12 1.12 Truth Table for a Logical (Boolean) ExpressionNoting that the logical AND, OR, and NOT operators are &, |, and ∼ in MATLAB, complete and run the following script “make_truth_table.m” to construct the truth table for the following logical (Boolean) expression:(P1.12.1) that lists the (output) value of the expression for all possible values of the input variables.%make_truth_table.m % makes the truth table for a logical (Boolean) expression N=3; % Number of Boolean variables Inputs=dec2bin_my([0:2̂N-1]); %All possible values of input variables A=Inputs(:,1); B=Inputs(:,2); C=Inputs(:,3); truth_table=[Inputs (?A&?&C)?(B?∼C)]

13 1.13 Playing with PolynomialsPolynomial evaluation: polyval() Write a MATLAB statement to compute(P1.13.1) (b) Polynomial addition/subtraction by using compatible vector addition/subtractionWrite a MATLAB statement to add the following two polynomial coefficient vectors:(P1.13.2) (c) Polynomial multiplication: conv()Write a MATLAB statement to get the following product of polynomials:(P1.13.3) (d) Polynomial division – deconv()Write a MATLAB statement to get the quotient and remainder of the following polynomial division:(P1.13.4) (e) Routines for differentiation/integration of a polynomialWhat you see in the below box is the function ‘ poly_der(p)’, which gets a polynomial coefficient vector p (in the descending order) and outputs the coefficient vector pd of its derivative polynomial. Likewise, compose a function ‘ poly_int(p)’, which outputs the coefficient vector of the integral polynomial for a given polynomial coefficient vector.(cf) MATLAB has the built‐in functions ‘ polyder()’/‘ polyint()’ that can be used to find the derivative/integral of a polynomial.function pd=<b>poly_der</b><![CDATA[(p) %p: Vector of polynomial coefficients in descending order N=length(p); if N<=1, pd=0; % constant else for i=1:N-1, pd(i)=p(i)*(N-i); end end function pi=<b>poly_int</b><![CDATA[(p) % p: Vector of polynomial coefficients in descending order N=length(p); pi(N+1)=0; for i=1:N, pi(i)=p(?)/(N-?+1); end(f) Roots of a polynomial equation: roots()Write a MATLAB statement to get the roots of the following polynomial equation:(P1.13.5) You can check if the result is right, by using the MATLAB command ‘ poly()’, which generates a polynomial having a given set of roots.(g) Partial fraction expansion (PFE) of a rational polynomial function – residue()/ residuez()(i) The MATLAB function ‘ [r,p,k]=residue(B,A)’ finds the PFE for a ratio of given polynomials B(s)/A(s) as(P1.13.6a) which is good for taking the inverse Laplace transform. Use the function to find the PFE for(P1.13.7a) (ii) The MATLAB function ‘ [r,p,k]=residuez(B,A)’ finds the PFE for a ratio of given polynomials B(z)/A(z) as(P1.13.6b) which is good for taking the inverse z‐transform. Use the function to find the PFE for(P1.13.7b) (h) Piecewise polynomial – mkpp()/ ppval()Suppose we have an M × N matrix P, the rows of which denote M (piecewise) polynomials of degree (N − 1) for different (nonoverlapping) intervals with (M + 1) boundary points bb=[b(1) .. b(M+1)], where the polynomial coefficients in each row are supposed to be generated with the interval starting from x = 0. Then we can use the MATLAB command ‘ pp=mkpp(bb,P)’ to construct a structure of piecewise polynomials, which can be evaluated by using ‘ ppval(pp)’.Figure P1.13 shows a set of piecewise polynomials {p1(x+3), p2(x+1), p3(x −2)} for the intervals [ −3, −1], [ −1, 2], and [2, 4], respectively, where(P1.13.8) (i) Complete and run the upper part of the above MATLAB script “nm01p13h.m” to use ‘ mkpp()’/‘ ppval()’ for making this graph.(ii) Complete and run the lower part of the above MATLAB script “nm01p13h.m” to use ‘ polyval()’ for making this graph.(iii) (cf) You can type ‘ help mkpp’ to see a couple of examples showing the usage of ‘ mkpp()’.%nm01p13h.m % to plot the graph of a piecewise polynomial bb=[-3 -? 2 4]; % Boundary point vector % Matrix having three polynomial coefficient vectors % in its rows P=[1 0 0;-1 2 -?;1 0 -2]; pp=mkpp(??,P) xx=bb(1)+[0:1000]/1000*(bb(end)-bb(1)); % Whole range plot(xx,ppval(??,xx)), hold on% Alternative without using ppval() % to plot the polynomial curves one by one for i=1:size(P,1) xx=[0:100]/100*(bb(i+1)-bb(i)); plot(xx+??(i),polyval(?(i,:),xx),'r:') endFigure P1.13 The graph of piece‐wise polynomial functions.

14 1.14 Routine for Matrix MultiplicationAssuming that MATLAB cannot perform direct multiplication on vectors/matrices, supplement the following incomplete function l.l ‘ multiply_matrix(A,B)’ so that it can multiply two matrices given as its input arguments only if their dimensions are compatible, but display an error message if their dimensions are not compatible. Try it to get the product of two arbitrary 3 × 3 matrices generated by the command ‘ rand(3)’ and compare the result with that obtained by using the direct multiplicative operator *. Note that the matrix multiplication can be described as(P1.14.1) function C=<b>multiply_matrix</b><![CDATA[(A,B) [M,K]=size(A); [K1,N]=size(B); if K1∼=K error('The # of columns of A is not equal to the # of rows of B') else for m=1:? for n=1:? C(m,n)=A(m,1)*B(1,n); for k=2:? C(m,n)=C(m,n)+A(m,k)*B(k,n); end end end end

15 1.15 Function for Finding Vector NormAssuming that MATLAB does not have the ‘ norm()’ command finding us the norm of a given vector/matrix, make a routine ‘ norm_vector(v,p)’ which computes the norm of a given vector as(P1.15.1) for any positive integer p, finds the maximum absolute value of the elements for p=inf, and computes the norm as if p=2, even if the second input argument p is not given. If you have no idea, permutate the statements in the below box and save it in a file named “norm_vector.m”. Additionally, try it to get the norm with p=1, 2,∞ ( inf) for an arbitray vector generated by the command ‘ rand(2,1)’. Compare the result with that obtained by using the ‘ norm()’ command.function nv=<b>norm_vector</b><![CDATA[(v,p) if nargin<2, p=2; end nv= sum(abs(v).̂p)̂(1/p); nv= max(abs(v)); if p>0&p∼=inf elseif p==inf end

16 1.16 Backslash (\) OperatorLet us play with the backslash ( \) operator.Use the backslash ( \) command, the minimum‐norm solution (2.1.7), and the ‘ pinv()’ command to solve the following equations, find the residual error ||Aix ‐bi||'s and the rank of the coefficient matrix Ai, and fill in Table P1.16 with the results.(P1.16.1) (P1.16.2) (P1.16.3) (b) Use the backslash( \) command, the least‐squares (LSs) solution (2.1.10), and the ‘ pinv()’ command to solve the following equations and find the residual error ||Aix ‐bi||'s and the rank of the coefficient matrix Ai, and fill in Table P1.16 with the results.(P1.16.4) (P1.16.5) (P1.16.6) (cf) If some or all of the rows of the coefficient matrix A in a system of linear equations can be expressed as a linear combination of other row(s), the corresponding equations are dependent, which can be revealed by the rank deficiency, i.e. rank(A) < min(M,N) where M and N are the row and column dimensions, respectively. If some equations are dependent, they may have either inconsistency (no exact solution) or redundancy (infinitely many solutions), which can be distinguished by checking if augmenting the RHS vector b to the coefficient matrix A increases the rank or not, that is, rank([Ab]) > rank(A) or not [M-2].(c) Based on the results obtained in (a) and (b) and listed in Table P1.16, answer the following questions:Based on the results obtained in (a‐(i)), which one yielded the nonminimum‐norm solution among the three methods, i.e. the backslash ( \) operator, the minimum‐norm solution (2.1.7), and the ‘ pinv()’ command? Note that the minimum‐norm solution means the solution whose norm (||x||) is the minimum over the many solutions.Based on the results obtained in (a), which one is most reliable as a means of finding the minimum‐norm solution among the three methods?Based on the results obtained in (b), choose two reliable methods as a means of finding the LS solution among the three methods, i.e. the backslash ( \) operator, the LS solution (2.1.10), and the ‘ pinv()’ command. Note that the LS solution means the solution for which the residual error (||Ax−b||) is the minimum over the many solutions.Table P1.16 Results of operations with backslash( \) operator and ‘ pinv()’ command.backslash ( \)Minimum‐norm or least‐squares (LS) pinv()Remarkx||Aix − bi||x||Aix − bi||x||Aix − bi||rank(Ai) redundant/inconsistentA1x − b11.5000 0 1.50006.4047 × 10−15A2x − b20.3143 0.6286 0.94291.7889A3x − b3∞ ∞ ∞∞A4x − b42.5000 0.00001.2247A5x − b5A6x − b6

17 1.17 Operations on VectorsFind the mathematical expression for the computation to be done by the following MATLAB statements: >n=0:100; S=sum(2.̂-n)Write a MATLAB statement which performs the following computation:Write a MATLAB statement that uses the commands ‘ prod()’ and ‘ sum()’ to compute the product of the sums of each row of a 3 × 3 random matrix.How does the following function ‘ repetition(x,M,N)’ convert a scalar or a matrix given as the first input argument x to make a new sequence y? How does it compare with the MATLAB built‐in function ‘ repmat()’?function y=<b>repetition</b><![CDATA[(x,M,N) y1=x; for n=2:N, y1 = [y1 x]; end y=y1; for m=2:M, y = [y; y1]; endComplete the following function ‘ zero_insertion(x,M,m)’ so that it can insert m zeros just after every Mth element of a given row vector sequence x to make a new sequence. Write a MATLAB statement to use the function for inserting two zeros just after every third element of x = [1 3 7 2 4 9] to gety = [1 3 7 0 0 2 4 9 0 0]function y=<b>zero_insertion</b><![CDATA[(x,M,m) Nx=length(x); N=floor(Nx/M); y=[reshape(x(1:M*N),?,N); ????s(m,N)]; y=[y(:)' x(M*N+1:Nx)];How does the following function ‘ zeroing(x,M,m)’ convert a given row vector sequence x to make a new sequence y?function y=<b>zeroing</b><![CDATA[(x,M,m) m=mod(m,M); Nx=length(x); N=floor(Nx/M); y=x; y(M*[1:N]-m)=0; Complete the following function ‘ sampling(x,M,m)’ so that it 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 ‘ besselj(k,beta)’ computes the value of the Bessel function of the first kind, Jk(β), of order k even when k is not an integer and β is a complex number.

24 1.24 Find the four functions in Chapters 5 and 7 that are fabricated in a recursive (self‐calling) structure.(cf) Does not those algorithms, which are the souls of the functions, seem to have been born to be in a nested structure?

25 1.25 Avoiding Runtime Error in Case of Deficient/Non‐admissible Input ArgumentsConsider the MATLAB function ‘ rotate_r(x,M)’, that you made in Problem 1.17(h). Does it work somehow when the user gives a negative integer as the second input argument M? If not, modify it so that it can perform the rotation left by – M samples for M<0, say, making rotate_r([1 2 3 4 5],-2)=[3 4 5 1 2]Consider the function ‘ trpzds(f,a,b,N)’ in Section 5.6, which computes the integral of function f over [ a, b] by dividing the integration interval into N sections and applying the trapezoidal rule. If the user tries to use it without the fourth input argument N, will it work? If not, make it work with N = 1000 by default even without the fourth input argument N.function INTf=trpzds(f,a,b,N) %integral of f(x) over [a,b] by trapezoidal rule with N segments if abs(b-a)<eps|N<=0, INTf=0; return; end h=(b-a)/N; x=a+[0:N]*h; fx=feval(f,x); values of f for all nodes INTf= h*((fx(1)+fx(N+1))/2+sum(fx(2:N))); % Eq.(5.6.1)

26 1.26 Parameter Passing Through vararginConsider the integration function ‘ trpzds(f,a,b,N)’ in Section 5.6. Can you use it to compute the integral of a function with some parameter(s), like the ‘ Bessel_integrand(x,beta,k)’ that you defined in Problem 1.23? If not, modify it so that it works for a function with some parameter(s) (see Section 1.3.6) and save it in the M‐file named “trpzds_par.m”. Then replace the ‘ quad()’ statement in the script “nm01p23b.m” (presented in Problem 1.23) by an appropriate ‘ trpzds_par()’ statement (with N= 1000) and run the script. What is the discrepancy between the integration results obtained by this function and the nested computing based on Eq. (P1.23.4b)? Is it comparable with that obtained with ‘ quad()’? How do you compare the running time of this function with that of ‘ quad()’? Why do you think it takes so much time to execute the ‘ quad()’ function?

27 1.27 Adaptive Input Argument to Avoid Runtime Error in Case of Different Input ArgumentsConsider the integration function ‘ trpzds(f,a,b,N)’ in Section 5.6. If some user tries to use this function with the following statement, will it work? trpzds(f,[a b],N) or trpzds(f,[a b])function INTf=<b>trpzds_bnd</b><![CDATA[(f,a,b,N) if numel(a)==2 if nargin>2, N=b; else N=1000; end b=a(?); a=a(?); else if nargin<4, N=1000; end end % . . . . . . . . . . . . . . . . . . . . .If not, modify it so that it works for such a usage (with a bound vector as the second input argument) as well as for the standard usage and save it in the M‐file named “trpzds_bnd.m”. Then try it to find the integral of e−t for [0,100] by running the following statements. What did you get?

28 1.28Continuous‐time Fourier transform (CtFT) of a SignalConsider the following definitions of CtFT and inverse continuous‐time Fourier transform (ICtFT) [W-8].(P1.28.1a) (P1.28.1b) Complete the following two MATLAB functions, ‘ CtFT1(x,Dt,w)’ computing the CtFT (P1.28.1a) of x(t) over [ ‐Dt,Dt] for w and ‘ ICtFT1(X,Bw,t)’ computing the ICtFT (P1.28.1b) of X(w) over [‐Bw, Bw] for t. You can choose whatever integral function including ‘ trpzds_par()’ (Problem 1.25) and ‘ quad()’, considering the running time.The following script “nm01p28.m” finds the CtFT of a rectangular pulse (with duration [ −1,1]) defined by ‘ rDt()’ for ω = [ −6π,+6π] and the ICtFT of a sinc spectrum (with bandwidth 2π) defined by ‘ sincBw()’ for t = [ −5,+5]. After having saved the functions and script into M‐files with the appropriate names, run the script to see the rectangular pulse, its CtFT spectrum, a sinc spectrum, and its ICtFT as shown in Figure P1.28. If it does not work, modify/supplement the functions so that you can rerun it to see the signals and their CtFT spectra.Figure P1.28 Graphs for Problem 1.28. (a1) A rectangular pulse function rD(t); (a2) the ICtFT of XB(ω); (b1) the CtFT spectrum of rD(t); and (b2) a sinc spectrum XB(ω).function Xw=<b>CtFT1</b><![CDATA[(x,Dt,w) % CtFT (Continuous-time Fourier Transform) x_ejkwt=inline([x '(t).*exp(-j*w*t)'],'t','w'); Xw=trpzds_par(x_ejkwt,-Dt,??,1000,?);function xt=<b>ICtFT1</b><![CDATA[(X,Bw,t) % ICtFT (Inverse Continuous-time Fourier Transform) Xejkwt=inline([X '(w).*exp(j*w*t)'],'w','t'); xt=trpzds_par(Xejkwt,-??,Bw,1000,?)/2/pi;%nm01p28.m : CtFT and ICtFT global B D % CtFT of A Rectangular Pulse Function t=[-50:50]/10; % Time vector w=[-60:60]/10*pi; % Frequency vector D=1; % Duration of a rectangular pulse rD(t) for k=1:length(w) Xw(k)=CtFT1('rDt',D*5,w(k)); end subplot(221), plot(t,rDt(t)); subplot(222), plot(w,abs(Xw)) % ICtFT of a Sinc Spectrum B=2*pi; % Bandwidth of a sinc spectrum sncB(w) for n=1:length(t) xt(n)=ICtFT1('sincBw',B*5,t(n)); end subplot(223), plot(t,real(xt)); subplot(224), plot(w,sincBw(w))function x=<b>rDt</b><![CDATA[(t) % Rectangular pulse function global D x=(-D/2<=t?t<=D/2);function X=<b>sincBw</b><![CDATA[(w) % Sinc function global B X=2*pi/?*sinc(?/B);