Simple Scilab routine for Newton's method

Please note: the functions below ignore a number of subtleties at the aim of being simple and (hopefully) understandable.

Newton's method is used to solve a system of equations f(x)=0 (here x is a d-dimensional vector and f is a d-dimensional vector of functions). There are 3 functions below:

  1. The function f(x), which should be written by the user (a simple function of 2 variables is used in the code below).
  2. The function newtstep(x,h). This takes a single Newton step starting from x. Derivatives of f are computed using the simplest difference formula with "stepsize" h (in the 1 dimensional case this means approximating f'(x) by (f(x+h)-f(x))/h ).
  3. The function newt(x,h,epsilon). This runs Newton's method, starting from x, using h as stepsize in the derivatives and stopping when the norm of f is less than epsilon. DANGER: Newton's method might not converge, and no stopping criterion has been put in! 

function a=f(x)
a(1)=x(1)*(x(1)*x(1)-3*x(2)*x(2))-1;
a(2)=x(2)*(3*x(1)*x(1)-x(2)*x(2));
endfunction

function a=newtstep(x,h)
n=size(x,1);              // get the dimension of the COLUMN vector submitted
F=f(x);                   // compute f(x)
M=zeros(n,n);             // prepare to make the matrix of derivatives
for i=1:n 
  v=zeros(n,1);
  v(i)=h;
  M(:,i)=(f(x+v)-F)/h;
end
a=x-M\F;                  // do the Newton step
endfunction

function y=newt(x,h,epsilon)
y=x;
z=norm(f(y));             // stop when this is below epsilon
while z>epsilon
  y=newtstep(y,h);            
  z=norm(f(y));
end
endfunction

Example of use

I put the three functions above into a file called "newt". Here is the transcript of a scilab session in which I execute the file "newt" and then use the commands to find the three roots of f=0 (for the f above): 

                           ===========
                           S c i l a b
                           ===========
 
 
                          scilab-2.6
                  Copyright (C) 1989-2001 INRIA 
 
 
Startup execution:
  loading initial environment
 
-->exec newt;
 
-->newt([1.5;0.5],1e-6,1e-6) 
 ans  =
 
!   1.        !
! - 1.203E-10 !
 
-->newt([-1.0;0.8],1e-6,1e-6)
 ans  =
 
! - 0.5       !
!   0.8660254 !
 
-->newt([-1.0;-0.8],1e-6,1e-6)
 ans  =
 
! - 0.5       !
! - 0.8660254 !
 
-->


Back to main numerical methods course page
Back to my main teaching page
Back to my main page