Trying to solve Simultaneous equations in matlab, cannot work out how to format the functions

Question

I was given a piece of Matlab code by a lecturer recently for a way to solve simultaneous equations using the Newton-Raphson method with a jacobian matrix (I've also left in his comments). However, although he's provided me with the basic code I cannot seem to get it working no matter how hard I try. I've spent many hours trying to introduce the 'func' function but to no avail, frequently getting the message that there aren't enough inputs. Any help would be greatly appreciated, especially with how to write the 'func' function.

  function root = newtonRaphson2(func,x,tol)

% Newton-Raphson method of finding a root of simultaneous    
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.    
% USAGE: root = newtonRaphson2(func,x,tol)    
% INPUT:    
% func = handle of function that returns[f1,f2,...,fn].
% x = starting solution vector [x1,x2,...,xn].
% tol = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root = solution vector.

if size(x,1) == 1; x = x'; end % x must be column vector    
for i = 1:30    
   [jac,f0] = jacobian(func,x);    
   if sqrt(dot(f0,f0)/length(x)) < tol   
      root = x; return
   end   

   dx = jac\(-f0);   
   x = x + dx;    
   if sqrt(dot(dx,dx)/length(x)) < tol    
      root = x; return    
   end
end

error('Too many iterations')

function [jac,f0] = jacobian(func,x)

% Returns the Jacobian matrix and f(x).

h = 1.0e-4;    
n = length(x);    
jac = zeros(n);    
f0 = feval(func,x);    
for i =1:n    
   temp = x(i);    
   x(i) = temp + h;    
   f1 = feval(func,x);    
   x(i) = temp;    
   jac(:,i) = (f1 - f0)/h;   
end

The simultaneous equations to be solved are:

sin(x)+y^2+ln(z)-7=0
3x+2^y-z^3+1=0
x+y+Z-=0

with the starting point (1,1,1).

However, these are arbitrary and can be replaced with anything, I mainly just need to know the general format. Many thanks, I know this may be a very simple task but I've only recently started teaching myself Matlab.

score 1 · Answer 1 · answered Feb 06 '15 at 12:10

You need to create a new file called myfunc.m (or whatever name you like) which takes a single input parameter - a column vector x - and returns a single output vector - a column vector y such that y = f(x).

For example,

function y = myfunc(x)

    y = zeros(3, 1);
    y(1) = sin(x(1)) + x(2)^2 + log(x(3)) - 7;
    y(2) = 3*x(1) + 2^x(2) - x(3)^3 + 1;
    y(3) = x(1) + x(2) + x(3);

end

You can then refer to this function as @myfunc as in

>> newtonRaphson2(@myfunc, [1;1;1], 1e-6);

The reason for the special notation is that Matlab allows you to call a function with no parameters by omitting the parens () that follow it. So for example, Matlab interprets myfunc as you calling the function with no arguments (so it tries to replace it with its result) whereas @myfunc refers to the function itself, rather than its result.

Alternatively you can write a function directly using the @ notation, as in

>> newtonRaphson2(@(x) exp(x) - 3*x, 2, 1e-2)
ans =
     1.5315
>> newtonRaphson2(@(x) exp(x) - 3*x, 1, 1e-2)
ans =
     0.6190

which are the two roots of the equation exp(x) - 3 * x = 0.

Edit - as an aside, your professor has terrible coding style (if the code in your question is a direct copy-paste of what he gave you, and you haven't mangled it along the way). It would be better to write the code like this, with indentation making it clear what the structure of the code is.

function root = newtonRaphson2(func, x, tol)
% Newton-Raphson method of finding a root of simultaneous
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.
%
% USAGE: root = newtonRaphson2(func,x,tol)
%
% INPUT:
%   func = handle of function that returns[f1,f2,...,fn].
%   x    = starting solution vector [x1,x2,...,xn].
%   tol  = error tolerance (default is 1.0e4*eps).
%
% OUTPUT:
%   root = solution vector.

    if size(x, 1) == 1; % x must be column vector
        x = x';
    end 

    for i = 1:30
        [jac, f0] = jacobian(func, x);

        if sqrt(dot(f0, f0) / length(x)) < tol
            root = x; return
        end

        dx = jac \ (-f0);
        x  = x + dx;

        if sqrt(dot(dx, dx) / length(x)) < tol
            root = x; return
        end
    end

    error('Too many iterations')

end

function [jac, f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).

    h = 1.0e-4;
    n = length(x);
    jac = zeros(n);
    f0 = feval(func,x);

    for i = 1:n
        temp = x(i);
        x(i) = temp + h;
        f1 = feval(func,x);
        x(i) = temp;
        jac(:,i) = (f1 - f0)/h;
    end

end

Trying to solve Simultaneous equations in matlab, cannot work out how to format the functions

1 Answers1