Write a function in MATLAB that returns the value of a function f (where f needs to be written by the user)

Question

Newton's method for 2x2 case:

$$J(x_n)s_n= f(x_n)$$ then update $$x_{n+1} = x_n - s_n$$ where $J$ is the Jacobian matrix formed from $f$. And $x_n$ is a 2-D vector.

Write a MATLAB function to implement Newton's method for a system of 2 nonlinear equations.

It should have the calling sequence [roots, count, resids, history] = newton(func, x0, tolerance)

The first input argument is a function handle to a function that returns the value of the function $f$ and the Jacobian matrix $J$ at the position $x$

e.g. function [f, J] = arbitraryfunc(x)

...

end

The function arbitraryfunc needs to be written by the user, including the expression for f and J. (i.e. the program should allow the user to input any functions). There is no need for J to be found symbolically from f.

My question are: 1. I am not sure how to fill in the blank ... ? 2. How to return the Jacobian of the function (inputted by the user)?

A detailed explanation is very much appreciated. Thanks.

Answer 1

So you want to write a subroutine that returns the jacobian matrix at point $x\in \mathbb{R}^2$ of a function handle @f in matlab. Also you wanna do this numerically, not symbolically.

Here is a simple finite difference subroutine that does the job for f being a function handle:

function J=numjacobian(f,x)

n=length(x); % n=2 in your case, but can be higher 
fx=feval(f,x); % evaluate the function at point x
step=1e-6; % difference step you choose, can be 1e-10 if you like

for i=1:n
   xstep = x;
   % compute the i-th component partial derivative 
   % numerically using first order forward difference approximation
   xstep(i)=x(i)+step;
   J(:,i)=(feval(f,xstep)-fx)/step;
end;
end

The input f is a function handle, x is any point. The output J is the jacobian matrix at x. Example: $f = (x^2+y^2, x^2-y^2)$:

  >>f = @(x) [x(:,1).^2 + x(:,2).^2, x(:,1).^2 - x(:,2).^2]; 
  >>numjacobian(f,[1 2])
    ans =
      2.00000100036846          4.00000100064801
      2.00000099992437         -4.00000100064801

---

For the Symbolic Toolbox approach, please refer to my answer here in this question: Solving a system with Newton's method in matlab? . A more object-oriented approach would be like you mentioned, allowing user to input any type of function: inline function, symbolic, function handle, etc. Converting symbolic function to function handle can be done by matlabFunction in MATLAB. And inline functions are string, hence easy to be converted to a symbolic function.