Hints for Assignment 4

• Problem 1 and 2: To find the number of function evaluations, use an m-file like this:

  function y = f1(x)
  global nf            % declare nf as global variable
  y = x.*cos(x).^3;    % must use ".*","./",".^" since x may be vector
  nf = nf + length(x); % update total number of function evaluations

  Then you can use quad and find the number of function evaluations as follows:

  global nf           % declare nf as global (IMPORTANT!!)
  nf = 0              % reset nf
  q = quad('f1',0,1,.001) % or:  q = adquad('f1',0,1,.001)
  nf                  % show number of function evaluations

  Note that you have to reset nf before each call of quad.

• Problem 2:In the m-file adquad.m you can evaluate func at a point x using y=feval(func,x). Your m-file adquad.m should look like this:

  function q = adquad(func,a,b,tol)
  fa = feval(func,a);
  ...
  if ...
    q = ...
  else
    q = adquad(func,...) + adquad(func,...)
  end

  Note this implementation is inefficient since it re-evaluates the function at the same points several times. Once you've got adquad working, rewrite it as follows: Use the following adquad0.m

  function q = adquad0(func,a,b,tol)
  fa = feval(func,a); fb = feval(func,b);
  q = adquad1(func,a,b,tol,fa,fb);

  and write a function adquad1 which only uses one new function evaluation fc=feval(func,c) at the midpoint c. If the estimated error in adquad1 is too large, it calls itself again for the two subintervals, passing on the already computed function values in the last two arguments.

• Problems 5 and 6: Stop the Newton iteration when ||x^(k+1)-x^(k)||₂ becomes smaller than some tolerance.
• Problem 5: In Matlab you can plot the contours where, e.g., cos(x²) - sin(xy) - .5 = 0 for x in [-5,5] and y in [-4,4] as follows:

  [X,Y]=meshgrid(-5:.1:5,-4:.1:4);
  contour(X,Y,cos(X.^2)-sin(X.*Y)-.5,[0 0])  % use .*,./,.^ !

  You can superimpose several contours by using hold on after the first contour command.

• Problem 6: You need to write m-files f.m and jac.m for the function and the Jacobian. The file f.m should look like this:

  function w=f(v)
  w = zeros(19,1);
  w(1) = ...
  for i=2:18
    w(i) = ...
  end
  w(19) = ...

  The file jac.m should look like this:

  function A=jac(v)
  A = sparse(19,19); % 19x19 zero matrix as sparse structure
  ...          % row 1 
  for j=2:18   % rows 2,...,18
    A(j,j-1) = ...
    A(j,j)   = ...
    A(j,j+1) = ...
  end
  ...          % row 19

  Use plot(0:.05:1,[0;v;0]) to plot the angle of the rod vs. the arclength.

  For finding 5 different solutions: Try values of C between -2 and 2. For small values of C the vector v will converge to the zero vector, and the iteration stops with a vector v with very small entries (less than 1e-5). Besides the zero vector, you should find four more solutions with larger entries.

If you have