Applied Numerical Methods Using MATLAB. Won Y. Yang

      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];