Unformatted text preview: 9/24/2010 Function inputs may be functions. If fplot was not part of Matlab, we could easily create something very similar. function [ ] = DIYfplot( f, range ) %DIYfplot Do it yourself version of fplot. Plots a function over a range. % Inputs: f = function to be plotted (need not be able to handle vectors) % range = two element vector containing start and end of range % No outputs (like a C++ void function) n = 100; % number of points x = linspace (range(1), range(2), n); % generate x values % generate y values y = zeros ([1 n]); % preallocate for efficiency for k = 1 : n y(k) = f(x(k)); end plot (x, y); end Just like the built‐in fplot, our version can be used to plot file functions, built‐in functions, and anonymous functions. In the first two cases an [email protected] is required (see previous set of slides). >> DIYfplot( @PR, [0 7]) >> grid on >> xlabel ('Mach number'); >> ylabel ('Pressure Ratio'); >> DIYfplot( @sin, [0 4*pi]) >> f = @(x) x^3 + 2* x^2 + 7; >> DIYfplot( f, [‐5 5]) Note: We could also create our own version of fzero and will soon be doing exactly this. 1 9/24/2010 Another potentially useful function is shown below. It produces a table of function values. function [ ] = ftable( f, vector ) %ftable Produces table of function values. % Inputs: f = function to be tabulated % vector = vector of values for the independent variable % No Outputs (like a C++ void function) fprintf (' Independent Dependent\n'); % Note: this loop can be simplified for k = 1 : length(vector) x = vector(k); y = f(x); fprintf ('%12f%16f\n', x, y); end end Sample usage: >> ftable (@PR, 1:0.5:6)
Independent 1.000000 1.500000 2.000000 2.500000 3.000000 3.500000 4.000000 4.500000 5.000000 5.500000 6.000000 Dependent 1.892929 3.413275 5.640441 8.526136 12.060965 12 16.242001 21.068081 26.538665 32.653474 39.412352 46.815206 Improved Loop: The loop in the function can be reduced to: for x = vector fprintf ('%12f%16f\n', x, f(x)); end 2 9/24/2010 Problem (text 5.14): You buy a $25,000 piece of equipment for nothing down at $5,500 per year for 6 years. What interest rate are you paying? The formula below relates the payment amount (A) to the present worth of the item (P), the number of years (n), and the interest rate (i, expressed as a fraction, 1% = 0.01) A P i (1 i ) n (1 i ) n 1 If we knew P, n, and i we could easily calculate A, but this isn’t the problem. We know A, P, and n and want to calculate i. Ideally we’d manipulate the equation to obtain an expression for i in terms of A, P, and n (remember that analytic solutions are best when we can obtain them). i = somefunction (A, P, n) Since we can’t do this easily we will convert the problem into root finding form
P i (1 i ) n A0 (1 i ) n 1 and use numerical methods (e.g. fzero) to find the value of i that satisfies the equation. Step 1: Define the basic function. >> P = 25000; n = 6; >> calcA = @(i) P * (i .* (i + 1).^n) ./ ((i + 1).^n ‐ 1); Step 2: Plot the function to get some idea of its behaviour and the approximate location of the solution. >> fplot (calcA, [0.01 0.20]) % plot for interest rates from 1% to 20% Step 3: Define the function for root finding (the f in f(x) = 0). >> A = 5500; >> f = @(i) calcA(i) ‐ A; Step 4: Use fzero to locate the desired root. In this case the plot tells us that the answer is somewhere between 6% and 10%. >> i = fzero(f, [0.06 0.10]) i = 0.0856 >> calcA(i) % check that answer is correct (always a good idea) ans = 5.5000e+003 3 9/24/2010 The basic idea can be used to produce a useful function m‐file that calculate the interest rate for any given P, A, and n. function [ i ] = P5_14( P, A, n ) %P5_14 Calculate interest rate. % NOTE: Assumes that interest rate is >= 0.005. % Inputs: P = present value of item % A = annual payment % n ‐ number of years % Outputs: i = interest rate (as a fraction) f = @(i) (P * (i .* (i + 1).^n) ./ ((i + 1).^n ‐ 1)) ‐ A; maxI = A/P; i = fzero(f, [ 0.005 maxI ] ); end The search for the root cannot be made to start at zero because f is undefined at i = 0 (as the formula for calculating A does not work for i = 0). This difficulty can be overcome by producing a calcA function that properly deals with the special case of i = 0. function [ i ] = P5_14v2( P, A, n ) %P5_14v2 Calculates interest rate. % Inputs and outputs as before… f = @(i) calcA (P, n, i) ‐ A; maxI = A/P; i = fzero(f, [ 0 maxI ] ); end function [A] = calcA (P, n, i) % calcA Calculates payment amount. if i == 0 A = P / n; else A = P * (i .* (i + 1).^n) ./ ((i + 1).^n ‐ 1); end end 4 9/24/2010 If function calcA will only ever be used by function P5_14v2, it can be placed in the same m‐file. This makes it a subfunction (see text section 3.1.3). function [ i ] = P5_14v2( P, A, n ) %P5_14v2 Calculates interest rate. % Inputs and outputs as before… f = @(i) calcA (P, n, i) ‐ A; maxI = A/P; i = fzero(f, [ 0 maxI ] ); end function [A] = calcA (P, n, i) % A subfunction if i == 0 A = P / n; else A = P * (i .* (i + 1).^n) ./ ((i + 1).^n ‐ 1); end end Subfunctions: A function m‐file may contain multiple functions (one after the other). The first of these functions is the primary function • Function name is the same as the name of the file • Accessible from the command window and other function m‐files All of the other functions are subfunctions • Can be used within the primary function • Can be used within subfunctions in the same m‐file • NOT accessible from the command window and other function m‐files 5 9/24/2010 Problem: Find all of the points that satisfy both y = x2‐17x+60 and y = 50sin(x/2) (i.e. find all of the intersections of the curves defined by these equations). Step 1: Plot the two curves. >> f1 = @(x) 50 * sin(0.5 * x); >> f2 = @(x) x.^2 ‐ 17 * x + 60; >> x = linspace (0, 20, 100); >> y1 = f1(x); y2 = f2(x); >> plot (x, y1, 'r', x, y2, 'b') >> grid on The plot function will accept more than one pair of x and y vectors. This allows multiple plots to be placed on the same graph. Each pair of vectors may be followed by a string containing plotting options. Some of the permitted characters are: colour: ‘r’ = red, ‘b’ = blue, ‘k’ = black, ‘g’ = green line style: ‘‐’ = solid, ‘:’ = dotted, ‘‐‐’ = dashed, ‘‐.’ = dash dot data point markers: ‘x’ = crosses, ‘o’ = circles Options may be combined (e.g. ‘r:x’ gives a dotted red plot with a crosses) 6 9/24/2010 Note: The hold command is another way of placing multiple plots on the same graph. >> plot (x, y1, 'r‘) >> hold on >> plot (x, y2, 'b') >> hold off % be careful ‐ hold for a figure stays in effect until it is turned off >> grid on Step 2: Define the function for root finding (and perhaps plot it) >> f = @(x) f1(x) ‐ f2(x); % we want the points where the functions are equal >> figure (2) >> fplot (f, [0 20]) >> grid on; Step 3: Find the roots: >> x1 = fzero(f, [0 2]); >> x2 = fzero(f, [6 8]); >> x3 = fzero(f, [12 14]); >> x4 = fzero(f, [16 18]); Problem: The distance (C) between ends of a model leaf spring is to be 5.1875” when the spring is loaded. Under these conditions the height (H) of the spring is 0.25”. C = 5.1875 H = 0.25” Assuming that the spring is an arc of a circle, what is the length (L) of the spring and what is its radius of curvature (R)? Once one has derived the formula (the interesting part) the actual calculation is just “plug and chug”.
R H 2 (C / 2) 2 2H L 2 R sin 1 (C / 2 R ) The formula for L assumes that sin‐1(C/2R) is in radians. As in C++, all Matlab trigonometric functions work in radians. 7 9/24/2010 Applying the formula gives R = 13.58” and L = 5.220”. >> C = 5.1875; >> H = 0.25; >> R = (H^2 + (C/2)^2) / ( 2 * H) R = 13.5801 >> L = 2 * R * asin (C / (2 * R)) L = 5.2196 Problem (continued): When the spring is unloaded its height increases to 0.5”. Assuming that the spring still forms an arc of a circle, how far apart are the ends and what is the radius of curvature? C = ?? H = 0.5” The applicable formulae are: θ R H = 0.5” L =5.220 C 2 H / L (1 cos( )) / R L /(2 ) C 2 R sin The value of θ cannot be found directly. Instead root finding is required. Rearranging the equation gives: 2 H / L (1 cos( )) / 0
From the diagram it is evident that θ must lie between 0 and π/2 radians. It is not uncommon for a root to be bounded like this (consider the floating ball problem). 8 9/24/2010 Zero cannot actually be used as the lower bound in using fzero as the function is undefined at θ = 0 and fzero needs to be able to evaluate the function at both bounds. This problem can be overcome by using a very small value instead. >> H = 0.5; >> f = @(theta) (2 * H) / L ‐ (1 ‐ cos(theta))/theta; >> theta = fzero(f, [1e‐6, pi/2]); >> R = L / (2 * theta) R = 6.7259 >> C = 2 * R * sin (theta) C = 5.0896 Notes: scientific format constants are as for C++ (1e‐6 is 1 x 10‐6) “pi” is a built‐in constant The commands required to solve the problem can be packaged into a script file. C = 5.1875; H = 0.25; R = (H^2 + (C/2)^2) / ( 2 * H); L = 2 * R * asin (C / (2 * R)); fprintf ('The loaded radius is %f\n', R); fprintf ('The length is %f\n', L); H = 0.5; f = @(theta) (2 * H) / L ‐ (1 ‐ cos(theta))/theta; theta = fzero(f, [1e‐6, pi/2]); R = L / (2 * theta); C = 2 * R * sin (theta); fprintf ('The unloaded radius is %f\n', R); fprintf ('The ends of the spring are %f apart\n', C); 9 9/24/2010 It is also possible to package them as a function with no inputs and no outputs. function = spring2 () C = 5.1875; H = 0.25; R = (H^2 + (C/2)^2) / ( 2 * H); L = 2 * R * asin (C / (2 * R)); fprintf ('The loaded radius is %f\n', R); fprintf ('The length is %f\n', L); H = 0.5; f = @(theta) (2 * H) / L ‐ (1 ‐ cos(theta))/theta; theta = fzero(f, [1e‐6, pi/2]); R = L / (2 * theta); C = 2 * R * sin (theta); fprintf ('The unloaded radius is %f\n', R); fprintf ('The ends of the spring are %f apart\n', C); end Script File: Executed by using green arrow (in editor window) or by entering the name of the file (in the command window). Variables are shared with the command window. Facilitates debugging (good) Executing a script can cause undesired side effects (bad) Any required functions that cannot be implemented as anonymous functions must be placed in separate files (can be a nuisance) Function with no Inputs and no Outputs: Executed in exactly the same way. Variables are not shared with the command window. Makes debugging harder (bad) Avoids undesired side effects (good) File may include subfunctions (very useful) 10 ...
View Full Document
This note was uploaded on 10/18/2010 for the course ECOR2606 2606 taught by Professor Aaaa during the Spring '10 term at Carleton CA.
- Spring '10