updated lecture 07

lecture_07/bisect.m

@@ -0,0 +1,37 @@
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xrold = xr;
+  xr = (xl + xu)/2;
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});

lecture_07/car_payments.m

@@ -0,0 +1,17 @@
+function amount_left = car_payments(monthly_payment,price,apr,no_of_years,plot_bool)
+  interest_per_month = apr/12;
+  number_of_months = no_of_years*12;
+  principle=price;
+  P_vector=zeros(1,number_of_months);
+  for i = 1:number_of_months
+    principle=principle-monthly_payment;
+    principle=(1+interest_per_month)*principle;
+    P_vector(i)=principle;
+  end
+  amount_left=principle;
+  if plot_bool
+    plot([1:number_of_months]/12, P_vector)
+    xlabel('time (years)')
+    ylabel('principle amount left ($)')
+  end
+end

lecture_07/falsepos.m

@@ -0,0 +1,39 @@
+function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xrold = xr;
+  xr = (xl + xu)/2;
+  % xr = (xl + xu)/2; % bisect method
+  xr=xu - (func(xu)*(xl-xu))/(func(xl)-func(xu)); % false position method
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});

lecture_07/fzerosimp.m

@@ -0,0 +1,41 @@
+function b = fzerosimp(xl,xu)
+a = xl; b = xu; fa = f(a); fb = f(b);
+c = a; fc = fa; d = b - c; e = d;
+while (1)
+  if fb == 0, break, end
+  if sign(fa) == sign(fb) %If needed, rearrange points
+    a = c; fa = fc; d = b - c; e = d;
+  end
+  if abs(fa) < abs(fb)
+    c = b; b = a; a = c;
+    fc = fb; fb = fa; fa = fc;
+  end
+  m = 0.5*(a - b); %Termination test and possible exit
+  tol = 2 * eps * max(abs(b), 1);
+  if abs(m) <= tol | fb == 0.
+    break
+  end
+  %Choose open methods or bisection
+  if abs(e) >= tol & abs(fc) > abs(fb)
+    s = fb/fc;
+    if a == c %Secant method
+      p = 2*m*s;
+      q = 1 - s;
+    else %Inverse quadratic interpolation
+      q = fc/fa; r = fb/fa;
+      p = s * (2*m*q * (q - r) - (b - c)*(r - 1));
+      q = (q - 1)*(r - 1)*(s - 1);
+    end
+    if p > 0, q = -q; else p = -p; end;
+    if 2*p < 3*m*q - abs(tol*q) & p < abs(0.5*e*q)
+      e = d; d = p/q;
+    else
+      d = m; e = m;
+    end
+  else %Bisection
+    d = m; e = m;
+  end
+  c = b; fc = fb;
+  if abs(d) > tol, b=b+d; else b=b-sign(b-a)*tol; end
+  fb = f(b);
+end

lecture_07/incsearch.m

@@ -0,0 +1,37 @@
+function xb = incsearch(func,xmin,xmax,ns)
+% incsearch: incremental search root locator
+%   xb = incsearch(func,xmin,xmax,ns):
+%     finds brackets of x that contain sign changes
+%     of a function on an interval
+% input:
+%   func = name of function
+%   xmin, xmax = endpoints of interval
+%   ns = number of subintervals (default = 50)
+% output:
+%   xb(k,1) is the lower bound of the kth sign change
+%   xb(k,2) is the upper bound of the kth sign change
+%   If no brackets found, xb = [].
+if nargin < 3, error('at least 3 arguments required'), end
+if nargin < 4, ns = 50; end %if ns blank set to 50
+% Incremental search
+x = linspace(xmin,xmax,ns);
+f = func(x);
+nb = 0; xb = []; %xb is null unless sign change detected
+%for k = 1:length(x)-1
+%  if sign(f(k)) ~= sign(f(k+1)) %check for sign change
+%    nb = nb + 1;
+%    xb(nb,1) = x(k);
+%    xb(nb,2) = x(k+1);
+%  end
+%end
+sign_change = diff(sign(f));
+[~,i_change] = find(sign_change~=0);
+nb=length(i_change);
+xb=[x(i_change)',x(i_change+1)'];
+
+if isempty(xb) %display that no brackets were found
+  fprintf('no brackets found\n')
+  fprintf('check interval or increase ns\n')
+else
+  fprintf('number of brackets: %i\n',nb) %display number of brackets
+end