M-files and Newton's Method in MATLAB: Understanding Scripts and Functions, Study Guides, Projects, Research of Matlab skills

Download M-files and Newton's Method in MATLAB: Understanding Scripts and Functions and more Study Guides, Projects, Research Matlab skills in PDF only on Docsity! MATLAB: M-files; Newton’s Method Last revised : March, 2003 Introduction to M-files In this session we learn the basics of working with M-files in MATLAB, so called because they must use “.m” for their filename extension. There are two types of M-files: • A script is a collection of MATLAB commands which are executed just as if they had been typed in at the command line. When the script is finished running, any variables defined in the script are still available in the current MATLAB command-line environment. • A function accepts input arguments and returns output arguments, similar to functions in structured programming languages such as JAVA, C, and FORTRAN. Variables defined in a function are local to the function. The variables created in the function are not directly accessible from the command line once the function completes execution. The Windows version of MATLAB provides a text editor for creating and modifying M-files. Throughout this document it is assumed that the reader is running the student version of MATLAB6 on the Windows platform. In this case the default working directory is typically C:\matlab_sv12\work and the M-file editor will save the M-files to this directory unless the user specifies otherwise.1 A Sample Script Executing MATLAB commands from within a script, even for relatively simple tasks, makes it much easier to reproduce results later, or to rerun the same calculations with different parameter settings. At the same time the script serves as a record of your activity and can serve as a useful starting point for future projects. Start up MATLAB and from the command line enter the command pwd to confirm the name of the current working directory. >> pwd ans = C:\matlab_sv12\work % print name of working directory % M-files will be saved here Note: When copying commands from this document into your own M-files or to the command line, the percent sign (%) marks the rest of the line as a comment. It is not necessary to copy the comments. They are included here for the sake of clarity. To create a new M-file, select File→ New→M-file from the menu bar. The M-file editor starts up with an empty window. Begin copying the following MATLAB commands into the new M-file. p = [5 -8 -1 2]; x = [-5.0 : 0.01 : 5.0]; y = p(4)*x.ˆ3 + p(3)*x.ˆ2 + p(2)*x + p(1); h1 = plot(x, y); % coefficients for cubic polyn % x=[-5,-4.99,-4.98,...,4.99,5] % evaluate y=f(x) % plot y=f(x) and save the handle Question : Why do we include the period with the exponentiation operator (.ˆ) when calculating y ? 1Use the pwd command to confirm the current working directory. To change the current directory, enter cd dir_name from the command line, or choose MATLAB→ Current Directory from the MATLAB Launch Pad. 1 To save the new M-file to the disk, choose File→ Save from the editor’s pull-down menu. The pop-up file browser should put you in the current working directory. We will name this first script polyn1.m. Enter the filename and click on the Save button. To test the new script select Debug→ Run from the editor’s menu bar, or enter the name of the script from the Command Window: >> polyn1 % execute the script named polyn1.m If MATLAB reports any errors, return to the editor, make the necessary corrections as indicated from the error message, save the revisions, and rerun polyn1, either from the Command Window or using Debug→ Run. For −5 ≤ x ≤ 5 the extreme values of this cubic polynomial are approximately ±200. MATLAB chooses the plotting range in the y direction based on these maximum and minimum values. Return to the M-file editor and use MATLAB’s axis function to choose a more reasonable range in both x and y. Turning on the grid lines makes it easier to read values from the graph. These additions to polyn1.m are highlighted below using a bold red font. Edit the script, save the changes, and rerun the script. p = [5 -8 -1 2]; x = [-5.0 : 0.01 : 5.0]; y = p(4)*x.ˆ3 + p(3)*x.ˆ2 + p(2)*x + p(1); h1 = plot(x, y); axis([-3 3 -15 15]); grid on % coefficients for cubic polyn % x=[-5,-4.99,-4.98,...,4.99,5] % evaluate y=f(x) % plot y=f(x) and save the handle % set axis [xLo, xHi, yLo, yHi] % turn on grid lines To dress up the graph a bit, we overlay a plot of the x-axis in red and add a title. Recall the use of the hold command to enable overlaying of plots. The clf command is used to clear the current figure before starting a new set of plots. Add the MATLAB code as shown below (highlighted in bold red font), save the changes, and rerun the revised script. p = [5 -8 -1 2]; x = [-5.0 : 0.01 : 5.0]; y = p(4)*x.ˆ3 + p(3)*x.ˆ2 + p(2)*x + p(1); clf h1 = plot(x, y); axis([-3 3 -15 15]); grid on hold on h2 = plot([x(1) x(end)], [0 0], ’r-’); title(’Cubic Polynomial’) % coefficients for cubic polyn % x=[-5,-4.99,-4.98,...,4.99,5] % evaluate y=f(x) % clear current figure % plot y=f(x) and save the handle % set axis [xLo, xHi, yLo, yHi] % add grid lines % overlay additional plots % plot the x-axis (y=0) % add a title 2 x 0 x 1 f(x) Figure 1: Illustration of Newton’s method; one iteration This iteration method results in a sequence of approximations {x0, x1, . . . , xn}, with the next iteration deter- mined by the formula xn+1 = xn − f (xn) f ′(xn) . (3) The distance between successive approximations |xn+1 − xn| serves as an estimate for the error |α − xn| where α represents the exact solution f (α) = 0. This iteration method continues until the error estimate is smaller than some prescribed error tolerance. The program should also be written to terminate if the number of iterations exceeds a certain value. This protects against the possibility that the sequence of approximations is not converging. A script for Newton’s method From completing the first part of the tutorial there should be three M-files, a script polyn1.m and two functions, cubic.m and dcubic.m. The plot generated by polyn1 shows that the cubic polynomial 2x3 − x2 − 8x + 5 has three distinct roots, with one root in the interval (0, 1.5). We try to find this root using Newton’s method with the initial guess x0 = 1.2. The script polyn1.m will be used as a starting point to create a new script for executing Newton’s method. Begin by opening a new file in the MATLAB editor (File→ New→M-file). Copy and paste the contents of polyn1.m into this new window and save the new file under the name newton.m. The MATLAB code in newton.m will need the following revisions: 1. Delete or comment out the lines for calculating and plotting the derivative function. 2. Parameter initialization; set x0, the error tolerance, and maximum number of iterations. 3. A loop that repeats the Newton iteration step until a stopping condition is satisfied. 4. Display the results at the end of the program. The first part of the script newton.m is listed below with the new MATLAB code highlighted in bold font. This part of the script takes us as far as reproducing a plot of f (x), plotting the initial guess, and pausing until the program receives input from the keyboard3 After entering the code, verify that the script is working. 3The code for pausing the the program pause_for_user.m and animating the Newton iteration newton_plot.m should already be downloaded and saved in the same directory as the script newton.m. The URL is given on page 3. The calls to these functions are optional. The script will still run to completion if these lines are commented out, but without the visualization. 5 x0 = 1.2; err_tol = 1.0e-08; max_iter = 20; p = [5 -8 -1 2]; x = [-5.0 : 0.01 : 5.0]; y = cubic(x, p); clf h1 = plot(x, y); axis([-3 3 -15 15]); grid on hold on h2 = plot([x(1) x(end)], [0 0], ’r-’); title(’Cubic Polynomial’) plot(x0, 0, ’r.’, ’markersize’, 12); if ( pause_for_user(-1) < 0 ) return end % initial guess for root of f(x) % set error tolerance % maximum number of iterations % coefficients for cubic polynomial % x=[-5,-4.99,-4.98,...,4.99,5] % evaluate y=f(x) % clear current figure % plot y=f(x) and save the handle % set axis [xLo, xHi, yLo, yHi] % add grid lines % overlay additional plots % plot the x-axis (y=0) % add a title % plot the initial guess x0 % wait for input from the keyboard % the user may choose to quit here Now construct the main iteration loop using the while command. Each time through the loop the program calculates the next Newton iteration using the formula in equation (3) and prints 4 the new iteration xn+1 and the error estimate xn+1 − xn . Next the script calls the function newton_plot to provide a simple animation of Newton’s method. This program plots the tangent line, marks the next Newton iteration as the x-intercept of the tangent line, and pauses the program until there is input from the keyboard. The user can choose to quit after each iteration in which case, newton_plot returns cf = 0. After returning from newton_plot the program checks the stopping conditions. The program will terminate if the error estimate is below the required error tolerance err_tol, or if the number of iterations has exceeded max_iter. The new approximation is copied onto xold in preparation for the next iteration and the script continues looping until the flag cf is no longer equal to 1. Some additional variables need to be initialized before entering the main iteration loop. All of the code shown below is to be inserted at the end of the existing script. After entering this portion of the code, be sure to save the results and run the code to check for errors. If the script is working properly, the approximations will quickly converge to the root at x ≈ 0.64 and the script should terminate after about six iterations. n = 0; xold = x0; cf = 1; % initialize the iteration counter % xold is the current guess % cf=1 means continue iterating 4MATLAB’s fprintf function is used to display formatted output. The syntax and conversion specifications are very similar to the fprintf command in C. The first argument is the file identifier; fprintf(1,...) writes to standard output (the screen). The most common conversion specifications include %d for decimal integers, %f for floating point real numbers, %e for exponential format, and %s for character strings. A specification such as %.2e says to display floating point numbers in exponential format using exactly two digits to the right of the decimal point. The specification \n produces a linefeed. 6 while ( cf == 1 ) n = n + 1; f = cubic(xold, p); df = dcubic(xold, p); xnew = xold - (f/df); err_est = abs(xnew - xold); fprintf(1, ’ xnew= %.8e, err_est= %.2e\n’, xnew, err_est); cf = newton_plot(x, xold, xnew, f, df); if ( err_est < err_tol ) cf = 0; elseif ( n >= max_iter ) cf = -1; end xold = xnew; end % continue as long as cf=1 % increment loop counter % evaluate f(x_n) % evaluate f’(x_n) % next Newton approximation % error estimate % animation % stop if error is small enough % stop if too many iterations % copy new iteration to xold % end of while loop Before exiting, the script should print the most important results to the screen, including the final approxi- mation, the initial guess, the error estimate, and the final number of iterations. See the footnote on page 6 for a description of the fprintf command. hf = plot(xnew, 0, ’b.’, ’markersize’, 12); if ( cf == -1 ) fprintf(1, ’ Newton method failed to converge.\n’); end fprintf(1, ’\n’); fprintf(1, ’ n= %d,’, n); fprintf(1, ’ x0= %.2e,’, x0); fprintf(1, ’ x= %.8e\n’, xnew); fprintf(1, ’ err_est= %.2e,’, err_est); fprintf(1, ’ f(x)= %.2e\n’, cubic(xnew,p)); % plot final guess % cf=-1 means no convergence % blank line % number of iterations % initial guess % final approx; newline % error estimate % f(x_n); newline After completing the script newton.m run it from the Command Window or by selecting Debug→ Run from the menu in the MATLAB editor. If the script and functions are working correctly, you should get output very similar to the following: >> newton ? "Enter" to continue, "q" to quit << ... ? "Enter" to continue, "q" to quit << n= 6, x0= 1.200000e+000, x= 6.39221640e-001 err_est= 2.73e-011, f(x)= 0.00e+000 >> 7