Applied Numerical Methods Using MATLAB. Won Y. Yang

Чтение книги онлайн.

Читать онлайн книгу Applied Numerical Methods Using MATLAB - Won Y. Yang страница 31

Applied Numerical Methods Using MATLAB - Won Y. Yang

Скачать книгу

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

Скачать книгу