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
