Using MATLAB to write a function that implements Newton's method in two dimensions

Question

I am trying to write a function that implements Newton's method in two dimensions and whilst I have done this, I have to now adjust my script so that the input parameters of my function must be f(x) in a column vector, the Jacobian matrix of f(x), the initial guess x0 and the tolerance where the function f(x) and its Jacobian matrix are in separate .m files.

As an example of a script I wrote that implements Newton's method, I have:

n=0;            %initialize iteration counter     
eps=1;          %initialize error     
x=[1;1];        %set starting value

%Computation loop     
while eps>1e-10&n<100 
    g=[x(1)^2+x(2)^3-1;x(1)^4-x(2)^4+x(1)*x(2)];         %g(x)      
    eps=abs(g(1))+abs(g(2));                             %error     
    Jg=[2*x(1),3*x(2)^2;4*x(1)^3+x(2),-4*x(2)^3+x(1)];   %Jacobian     
    y=x-Jg\g;                                            %iterate     
    x=y;                                                 %update x     
    n=n+1;                                               %counter+1     
end 

n,x,eps       %display end values

So with this script, I had implemented the function and the Jacobian matrix into the actual script and I am struggling to work out how I can actually create a script with the input parameters required.

Thanks!

Answer 1

If you don't mind, I'd like to restructure your code so that it is more dynamic and more user friendly to read.

Let's start with some preliminaries. If you want to make your script truly dynamic, then I would recommend that you use the Symbolic Math Toolbox. This way, you can use MATLAB to tackle derivatives of functions for you. You first need to use the syms command, followed by any variable you want. This tells MATLAB that you are now going to treat this variable as "symbolic" (i.e. not a constant). Let's start with some basics:

syms x;
y = 2*x^2 + 6*x + 3;
dy = diff(y); % Derivative with respect to x.  Should give 4*x + 6;
out = subs(y, 3); % The subs command will substitute all x's in y with the value 3
                  % This should give 2*(3^2) + 6*3 + 3 = 39

Because this is 2D, we're going to need 2D functions... so let's define x and y as variables. The way you call the subs command will be slightly different:

syms x, y; % Two variables now
z = 2*x*y^2 + 6*y + x;
dzx = diff(z, 'x'); % Differentiate with respect to x - Should give 2*y^2 + 1
dzy = diff(z, 'y'); % Differentiate with respect to y - Should give 4*x*y + 6
out = subs(z, {x, y}, [2, 3]); % For z, with variables x,y, substitute x = 2, y = 3
                               % Should give 56

One more thing... we can place equations into vectors or matrices and use subs to simultaneously substitute all values of x and y into each equation.

syms x, y;
z1 = 3*x + 6*y + 3;
z2 = 3*y + 4*y + 4;
f = [z1; z2];
out = subs(f, {x,y}, [2, 3]); % Produces a 2 x 1 vector with [27; 25]

We can do the same thing for matrices, but for brevity I won't show you how to do that. I will defer to the code and you can see it then.

Now that we have that established, let's tackle your code one piece at a time to truly make this dynamic. Your function requires the initial guess x0, the function f(x) as a column vector, the Jacobian matrix as a 2 x 2 matrix and the tolerance tol.

Before you run your script, you will need to generate your parameters:

syms x y; % Make x,y symbolic
f1 = x^2 + y^3 - 1; % Make your two equations (from your example)
f2 = x^4 - y^4 + x*y;
f = [f1; f2]; % f(x) vector

% Jacobian matrix
J = [diff(f1, 'x') diff(f1, 'y'); diff(f2, 'x') diff(f2, 'y')];

% Initial vector
x0 = [1; 1];

% Tolerance:
tol = 1e-10;

Now, make your script into a function:

% To run in MATLAB, do: 
% [n, xout, tol] = Jacobian2D(f, J, x0, tol);
% disp('n = '); disp(n); disp('x = '); disp(xout); disp('tol = '); disp(tol);

function [n, xout, tol] = Jacobian2D(f, J, x0, tol)

% Just to be sure...
syms x, y;

% Initialize error
ep = 1; % Note: eps is a reserved keyword in MATLAB

% Initialize counter
n = 0;

% For the beginning of the loop
% Must transpose into a row vector as this is required by subs
xout = x0';

% Computation loop     
while ep > tol && n < 100 
   g = subs(f, {x,y}, xout);   %g(x)
   ep = abs(g(1)) + abs(g(2)); %error     
   Jg = subs(J, {x,y}, xout);  %Jacobian  
   yout = xout - Jg\g;         %iterate     
   xout = yout;                %update x     
   n = n + 1;                  %counter+1     
end

% Transpose and convert back to number representation
xout = double(xout');

I should probably tell you that when you're doing computation using the Symbolic Math Toolbox, the data type of the numbers as you're calculating them are a sym object. You probably want to convert these back into real numbers and so you can use double to cast them back. However, if you leave them in the sym format, it displays your numbers as neat fractions if that's what you're looking for. Cast to double if you want the decimal point representation.

Now when you run this function, it should give you what you're looking for. I have not tested this code, but I'm pretty sure this will work.

Happy to answer any more questions you may have. Hope this helps.

Cheers!