Applied Numerical Methods Using MATLAB. Won Y. Yang

Note that MATLAB has a built‐in graphic function ‘ ezplot()’, which is much more powerful and convenient to use than ‘ ez_plot()’. You can type ‘ help ezplot’ to see its function and usage.

1.3.7 Adaptive Input Argument List

      A MATLAB function/routine is said to be ‘adaptive’ to users in terms of input arguments if it accepts different number/type of input arguments and makes a reasonable interpretation. For example, let us see the nonlinear equation solver ‘ Newton()’ in Section 4.4. Its input argument list is

       (f,df,x0,TolX,MaxIter)

      where f, df, x0, tol, and MaxIter denote the name of function (to be solved), the name of the derivative function, the initial guess (for solution), the error tolerance, and the maximum number of iterations, respectively. Suppose the user, not having the derivative function, tries to use the routine with just four input argument as follows:

       >Newton(f,x0,tol,MaxIter)

      At first, these four input arguments will be accepted as f, df, x0, and tol, respectively. But when the second line of the program body

       if nargin==4&isnumeric(df), MaxIter=TolX; TolX=x0; x0=df; end

      is executed, the routine will notice something wrong from that df is not any function name but a number, and then interpret the input arguments as f, x0, tol, and kmax to the idea of the user. This allows the user to use the routine in two ways, depending on whether he is going to supply the routine with the derivative function or not. This scheme is conceptually quite similar to function overloading of C++, but C++ requires us to have several functions having the same name, with different argument list.

Problems

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