% discussed the use of FZERO and FMIN, particularly the how to specify the % function whose zero is being sought. % The help for INLINE says that INLINE(expr) will provide a function whose % value at a point is obtained by evaluating the expr (which must be a string). % If that string contains many variables, how does matlab know which one is % the independent variable? The help gives some rule >> a = 3.14; >> fname = inline('exp(a)-b') fname = Inline function: fname(a,b) = exp(a)-b % How will FZERO know what fname is a function of??? >> fzero(fname,1) ??? Error using ==> fzero FZERO cannot continue because user supplied inline object ==> fname failed with the error below. This may be due to changed syntax of FZERO. Please check HELP FZERO. Error using ==> inline/feval Not enough inputs to inline function. % We know that FZERO merely calls on FEVAL, so let's try it directly there: >> feval(fname,1) ??? Error using ==> inline/feval Not enough inputs to inline function. %% aha, the problem seems to be that our inline function mentions two variables %% a and b, and so it is not clear what variable to use. We can specify which %% variable we mean by saying so in an optional argument to INLINE: >> fname = inline('exp(a)-b','b') fname = Inline function: fname(b) = exp(a)-b >> fzero(fname,1) ??? Error using ==> fzero FZERO cannot continue because user supplied inline object ==> fname failed with the error below. This may be due to changed syntax of FZERO. Please check HELP FZERO. Error using ==> inline/feval Error in inline expression ==> exp(a)-b ??? Undefined function or variable 'a'. %% .., so the fact that we gave a a particular value earlier has no effect %% here. By the time FEVAL is given fname, nobody knows about the a . %% what if we use explicitly x for the independent variable? >> fname = inline('exp(x)-b') fname = Inline function: fname(b,x) = exp(x)-b fzero(fname,1) ??? Error using ==> fzero FZERO cannot continue because user supplied inline object ==> fname failed with the error below. This may be due to changed syntax of FZERO. Please check HELP FZERO. Error using ==> inline/feval Not enough inputs to inline function. %% Perhaps, INLINE can be given more than one statement??? fname = inline('b = 3.14; exp(x)-b') fname = Inline function: fname(b,x) = b = 3.14; exp(x)-b fzero(fname,1) ??? Error using ==> fzero FZERO cannot continue because user supplied inline object ==> fname failed with the error below. This may be due to changed syntax of FZERO. Please check HELP FZERO. Error using ==> inline/feval Not enough inputs to inline function. % So, that is not a solution, either. We conclude that the expression in the % inline statement better have only one `obvious' variable, and everything % else needs to be known. % This raises a related question: What if the function whose zero we are % seeking itself depends on some parameters? % For this, FZERO provides additional arguments. % E.g., fzero(fname,xo,[],p1,p2) % will work if fname is a function of three arguments, the first of which % is the independent variable, and the other two are parameters. % To be specific, let's construct a spline and then look for a zero of it: >> axis([0 10 -2 5]); grid >> xy = ginput(10) xy = 0.0461 1.4795 0.0230 1.9708 1.4286 -0.0556 2.4424 -1.4678 3.4562 -0.9766 4.7465 -0.4854 6.6129 0.3333 7.3272 2.1345 8.5023 3.7105 9.3088 4.1608 >> sp = spline(xy(:,1), xy(:,2)) sp = form: 'pp' breaks: [1x10 double] coefs: [9x4 double] pieces: 9 order: 4 dim: 1 >> fnplt(sp), grid %% This spline has a zero in the interval [6 .. 7], but the only way we know %% of evaluating it is at x is via ppval(sp,x). So, we cannot use >> fzero('ppval',6,[],sp) ??? Error using ==> fzero FZERO cannot continue because user supplied function ==> ppval failed with the error below. This may be due to changed syntax of FZERO. Please check HELP FZERO. Error using ==> diff Function 'diff' not defined for variables of class 'struct'.. %% .. but we can use INLINE as follows: >> fname = inline('ppval(sp,x)','x','sp') fname = Inline function: fname(x,sp) = ppval(sp,x) >> fzero(myfunc,6,[],sp) Warning: Cannot determine ... etc Zero found in the interval: [5.52, 6.48]. ans = 6.4321 %% ... or, we could make up an m-file, like the following >> type myppval function values = myppval(x,pp) values = ppval(pp,x) %% ... which, in effect, ensures that the independent variable is the first %% argument, and then supply it to FZERO: >> fzero(myppval,3,[],sp) ??? Input argument 'pp' is undefined. Error in ==> c:\courses\412\00\myppval.m On line 2 ==> values = ppval(pp,x); %% WHAT?? WHAT IS WRONG NOW??? %% have a look at myppval: >> edit myppval %% .. all looks well; so does a look at the help for PPVAL %% AHA, I didn't give FZERO the STRING 'myppval', but the VARIABLE myppval %% and that generated real difficulty for FEVAL. %% So, here is an important difference! The name we set INLINE equal to can %% be used directly in FZERO. The name of a function m-file must be given %% to FZERO in quotes. >> fzero('myppval',3,[],sp) Warning: Cannot determine ... Zero found in the interval: [5.52, 6.48]. ans = 6.4321 %% Now let's look at the history again: >> fzero('myppval',3,1,sp) ??? Error using ==> > Function '>' not defined for variables of class 'struct'. Error in ==> C:\MATLABR11\toolbox\matlab\funfun\fzero.m On line 127 ==> if trace > 1 %% On checking HELPDESK for FZERO, we find that the third argument better %% be a structure. We set its one and only field, Display: >> opt.Display = 'iter'; >> fzero('myppval',3,opt,sp) Func-count x f(x) Procedure 1 6 -0.476016 initial 2 5.83029 -0.563306 search 3 6.16971 -0.337088 search 4 5.76 -0.586421 search 5 6.24 -0.262335 search 6 5.66059 -0.607793 search 7 6.33941 -0.137686 search 8 5.52 -0.61854 search 9 6.48 0.079614 search Looking for a zero in the interval [5.52, 6.48] 10 6.37053 -0.0938679 interpolation 11 6.42976 -0.00381583 interpolation 12 6.43206 -0.000143295 interpolation 13 6.43215 1.73833e-08 interpolation 14 6.43215 -1.25061e-12 interpolation 15 6.43215 -2.22045e-16 interpolation 16 6.43215 3.88578e-15 interpolation Zero found in the interval: [5.52, 6.48]. ans = 6.4321 %%%%%%%% then we discussed finding minima. % Bisection now becomes Golden Search, as described in the book. % Methods corresponding to Newton's method use a parabola as a local model % whose unique minimum then gives the next approximation. Again, matlab % uses a hybrid method, in FMIN: % This time, interpolate to sin: >> x = linspace(0,2*pi,11); >> sp = spline(x,sin(x)); >> fnplt(sp) %% from the graph, we see a minimum in the interval [0 .. 7], so ... >> fmin(myppval,0,7,opt,sp) ??? Input argument 'pp' is undefined. Error in ==> c:\courses\412\00\myppval.m On line 2 ==> values = ppval(pp,x); %% Well, I forgot again to put the name of the m-file myppval into quotes :-) >> fmin('myppval',0,7,opt,sp) ??? Conversion to double from struct is not possible. Error in ==> C:\MATLABR11\toolbox\matlab\funfun\foptions.m On line 51 ==> OPTIONS(1:sizep)=parain(1:sizep); Error in ==> C:\MATLABR11\toolbox\matlab\funfun\fmin.m On line 40 ==> options = foptions(options); %% I used helpdesk to check the help for FMIN and was told there to use %% FMINBND instead (it turns out that FMIN expects an array for that fourth %% argument, opt, while FMINBND expects the same kind of structure that FZERO %% expects). >> fminbnd('myppval',0,7,opt,sp) Func-count x f(x) Procedure 1 2.67376 0.450925 initial 2 4.32624 -0.926291 golden 3 5.34752 -0.803922 golden 4 4.66884 -0.99879 parabolic 5 4.7142 -0.999739 parabolic 6 4.71281 -0.99974 parabolic 7 4.71274 -0.99974 parabolic 8 4.71277 -0.99974 parabolic 9 4.7127 -0.99974 parabolic Optimization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 ans = 4.7127 %% We notice from the history that FMINBND first looks for a suitable %% interval, using Golden Search, and then uses parabolic interpolation to zero %% in on the minimum. % We know that this should be (3/2)pi, so check it: >> ans - pi*(3/2) ans = 3.4789e-004 %... that's close enough, particularly since we are finding only a minimum of % the spline interpolant to sin . % Why not just finding a zero of f' instead of a minimum of f??? % Would need a program to compute values of f' (well, you know how to % differentiate a program :-). Also, a zero of f' is NOT guaranteed to be a % minimum; it could also be a maximum, or neither. So, after finding a zero of % f', would have to do further checking. %% finally, briefly discussed using elimination to produce the factors L and %% U of an appropriately row-permuted matrix A, using either partial pivoting %% or scaled partial pivoting. Here is a 2-by-2 example >> A = [2 5; >> 1 1]; % put A into a working array: >> W = A; n = 2; % initialize the vector for recording the interchanges: >> piv = 1:n; % start elimination: >> k = 1; % pick pivot row for xk by partial pivoting: since 2>1, it's the first row, % so, no interchange of rows. % Multiplier(s): >> W(k+1:n,1) = W(k+1:n,1)/W(k,k); % Elimination: >> W(k+1:n,k+1:n) = W(k+1:n,k+1:n) - W(k+1:n,1)*W(k,k+1:n); % Since there are only two unknowns, we are done. >> U = triu(W); L = eye(n) + tril(W,-1); >> P = eye(n); P = P(piv,:); % Check >> P*A - L*U ans = 0 0 0 0 % For scaled partial pivoting, would first compute the sizes: >> s = max(abs(A),[],2) s = 5 1 % put A into a working array: >> W = A; n = 2; % initialize the vector for recording the interchanges: >> piv = 1:n; % start elimination: >> k = 1; % pick pivot row for xk by scaled partial pivoting: % since 2/5 < 1/1, it's the second row, % so, we must interchange >> W([2,1],:) = W([1 2],:); piv([2 1]) = piv([1 2]); % Multiplier(s): >> W(k+1:n,1) = W(k+1:n,1)/W(k,k); % Elimination: >> W(k+1:n,k+1:n) = W(k+1:n,k+1:n) - W(k+1:n,1)*W(k,k+1:n); % Since there are only two unknowns, we are done. >> U = triu(W); L = eye(n) + tril(W,-1); >> P = eye(n); P = P(piv,:); % Check >> P*A - L*U ans = 0 0 0 0