Applied Numerical Methods Using MATLAB. Won Y. Yang
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)