Unformatted text preview: First some leftover from Monday. We ran into difficulty using the forward function. function der=forward(f,x,h,varargin) %Calculates the derivative of the function f with respect to x using a %forward difference approximation with step size h. varargin includes %variables used in the call to f der=(f(x+h,varargin{:})f(x,varargin{:}))/h; The function includes the varargin structure to allow additional parameters to be passed to the function. However, apparently the cost is that if there are no parameters at all, it does not work. forward(asin,0.5,0.1) ??? Error using ==> asin Not enough input arguments. Error in ==> asin at 15 [varargout{1:nargout}] = builtin('asin', varargin{:}); This can be remedied by using a function without the varargin or by inventing a dummy parameter >> [email protected](x,a) asin(x) asinn = @(x,a) asin(x) >> forward(asinn,0.5,0.1,0.1) ans = 1.1990 Then we can see the benefits of central differences by averaging the forward and backwards differences >> der=(forward(asinn,0.5,0.1)+forward(asinn,0.5,0.1))/2 der = 1.1599 >> der=(forward(asinn,0.5,0.01)+forward(asinn,0.5,0.01))/2 der = 1.1548 So that we are much closer to the exact derivative of 1.1547. Vectorization: There was a question after class about one of the homework problems involving vectorization. The sign function is of great help in promoting vectorization. For example, say we want a function basket(x) that will be zero if x is smaller in absolute magnitude than 0.3, and is equal to x otherwise. Say you want to plot basket(sin(x)) for 1<x<1. >> x=1:0.1:1 x = Columns 1 through 11 1.0000 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1000 0 Columns 12 through 21 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 >> y=sin(x) y = Columns 1 through 11 0.8415 0.7833 0.7174 0.6442 0.5646 0.4794 0.3894 0.2955 0.1987 0.0998 0 Columns 12 through 21 0.0998 0.1987 0.2955 0.3894 0.4794 0.5646 0.6442 0.7174 0.7833 0.8415 >> sig=sign(abs(x)0.3) sig = Columns 1 through 19 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Columns 20 through 21 1 1 >> basket=abs(y).*(sig+1)/2 basket = Columns 1 through 11 0.8415 0.7833 0.7174 0.6442 0.5646 0.4794 0.3894 0 0 0 0 Columns 12 through 21 0 0 0 0.3894 0.4794 0.5646 0.6442 0.7174 0.7833 0.8415 >> plot(x,basket) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 Chapter 5: Bracketing methods. We have already played with the bisection method. function [xm,ym]= cosbisect_while(xl,xu,accuracy) %performs bisection towards finding the root of cos(x) in the %interval (xl,xu), until the interval is smaller than accuracy and outputs the value of the %midpoint and the function while (xuxl)>accuracy xm=(xl+xu)/2; ym=cos(xm); if sign(ym)<0 xu=xm; else xl=xm; end end >> cosbisect_while(0,2,0.01) ans = 1.5703 >> The textbook also discusses the False Position method, which is a type of secant method. Matlab has a builtin function fzero, which combines bisection, False Position, and another method called inverse quadratic interpolation to home in on a root. Here are the first few line that we get from help command on fzero help fzero FZERO Scalar nonlinear zero finding. X = FZERO(FUN,X0) tries to find a zero of the function FUN near X0, if X0 is a scalar. It first finds an interval containing X0 where the function values of the interval endpoints differ in sign, then searches that interval for a zero. FUN accepts real scalar input X and returns a real scalar function value F, evaluated at X. The value X returned by FZERO is near a point where FUN changes sign (if FUN is continuous), or NaN if the search fails. X = FZERO(FUN,X0), where X0 is a vector of length 2, assumes X0 is an interval where the sign of FUN(X0(1)) differs from the sign of FUN(X0(2)). An error occurs if this is not true. Calling FZERO with an interval guarantees FZERO will return a value near a point where FUN changes sign. >> fzero(@cos,1,2) ans = 1.5708 If we want to follow the progress of fzero, we can create a function that calculates cos(x) but displays its value. function [cs]=pcos(x) cs=cos(x); x,cs We now call fzero and get a log of how it found the root. Note that because it uses the secant method and quadratic interpolation it gets to the right answer to five digits in about four iterations. Try it with only one value instead of an interval, and you will see that it has much harder time. >> fzero(@pcos,[0.5,2]) x = 0.5000 cs = 0.8776 x = 2 cs = 0.4161 x = 1.5175 cs = 0.0533 x = 1.5723 cs = 0.0015 x = 1.5708 cs = 6.7207e007 x = 1.5708 cs = 2.3841e013 x = 1.5708 cs = 6.1232e017 x = 1.5708 cs = 6.0490e016 x = 1.5708 cs = 6.1232e017 ans = 1.5708 ...
View
Full Document
 Spring '09
 RAPHAELHAFTKA
 fzero, Rootfinding algorithm

Click to edit the document details