Keywords - Function groups - @ A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

Library: nummath
See also: nmjacobian

Quantlet: nmnewton
Description: Newton-Raphson method for solving system func(x)=0

Usage: z = nmnewton(func, fder, x0{, epsf, epsx, maxiter, checkcond})
Input:
func name of function(s) (string or vector of strings); strings describe functions func_i, their common root is to be found; format of func should be: 1. func is one string describing a function whose output is n x 1 vector (in this case, derivative of func, fder, has to be given) 2. func is a vector of strings, each of them describes one of functions in the system (i.e. its output is a real number) The function should have just one parameter x, which is: a. vector n x 1, if the fder is given (see below) b. matrix n x k (for any k >= 1) containing x = (x1,x2,...,xk), if the jacobian is to be computed automatically; xi (n x 1 vector; i=1..k) represent points in which the function should be evaluated.
fder derivatives of func; there are several possible formats: 1. name of function(s) (string or vector of strings); strings describe gradients of func_i; fder is either one string describing the function computing jacobian of func and its output is n x n vector, or fder is a vector of strings, each of them describing the gradient of one of functions in the system (i.e. its output is n x 1 vector) 2. empty string; in this case the jacobian will be computed automatically using the quantlet nmjacobian with a default step h 3. a number or n x 1 vector h; the jacobian will be computed automatically using the quantlet nmjacobian with the given step h
x0 n x 1 vector, the initial estimate of the solution
epsf optional number; iterative process ends if sum of absolute values of func will be less or equal to epsf; default value is 1e-7
epsx optional number; iterative process ends if sum of absolute values of corrections to x will be less or equal to epsx; default value is 1e-7
maxiter optional number; maximal number of iterations; default value is 100
checkcond optional number; if given, the jacobian will be tested in each step; a warning will be written if its condition number is greater then checkcond
Output:
z.x n x 1 vector, estimated solution
z.fval n x 1 vector, values of func in the point x
z.fderval n x n matrix, values of gradients of func in the point x

Example:
library("nummath")
;
; definition of vector function f(x1,x2) =(x1+x2,x1-x2)
; and its jacobian
;
proc(func)=f(x)
  func=#(x[1,]+x[2,],x[1,]-x[2,])
endp
proc(fder)=d(x)
  fder=#(1,1)~#(1,-1)
endp
;
nmnewton("f","d",#(1,1))

Result:
Contents of _tmp.x
[1,]        0
[2,]        0

Contents of _tmp.fval
[1,]        0
[2,]        0

Contents of _tmp.fderval
[1,]        1        1
[2,]        1       -1
Example:
library("nummath")
;
; definition of function f(x1,x2) =(sin(x1+x2),cos(x1-x2))
; and of gradients of its components
;
proc(func)=f(x)
  func=#(sin(x[1,]+x[2,]),cos(x[1,]-x[2,]))
endp
proc(fder1)=d1(x)
  fder1=#(cos(x[1,]+x[2,]),cos(x[1,]+x[2,]))
endp
proc(fder2)=d2(x)
  fder2=#(-sin(x[1,]-x[2,]),sin(x[1,]-x[2,]))
endp
;
nmnewton("f",#("d1","d2"),#(1,-1))
; the nearest solution is(pi/4,-pi/4)

Result:
Contents of _tmp.x
[1,]   0.7854
[2,]  -0.7854

Contents of _tmp.fval
[1,]        0
[2,]  6.1257e-17

Contents of _tmp.fderval
[1,]        1        1
[2,]       -1        1
Example:
library("nummath")
;
; definition of functions f1(x1,x2) = x1^2 - 2*x2   - 2
;                         f2(x1,x2) = x1   -   x2^3 - 1
; and of their gradients
;
proc(func1)=f1(x)
  func1=x[1,]^2-2*x[2,]-2
endp
proc(func2)=f2(x)
  func2=x[1,]-x[2,]^3-1
endp
proc(fder1)=d1(x)
  fder1=#(2*x[1,],-2)
endp
proc(fder2)=d2(x)
  fder2=#(1,-3*x[2,]^2)
endp
sol = nmnewton(#("f1","f2"),#("d1","d2"),#(5,5),1e-15,1e-15,200)
sol

Result:
Contents of sol.x
[1,]        2
[2,]        1

Contents of sol.fval
[1,]        0
[2,]        0

Contents of sol.fderval
[1,]        4       -2
[2,]        1       -3
Example:
library("nummath")
;
; definition of functions f1(x1,x2) = x1^2 - 2*x2   - 2
;                         f2(x1,x2) = x1   -   x2^2 - 1
; jacobian will be computed automatically
;
proc(func1)=f1(x)
  func1=x[1,]^2-2*x[2,]-2
endp
proc(func2)=f2(x)
  func2=x[1,]-x[2,]^3-1
endp
sol = nmnewton(#("f1","f2"),"",#(0,0))
sol

Result:
Contents of sol.x
[1,] -1.7807e-15
[2,]       -1

Contents of sol.fval
[1,]  1.5987e-14
[2,]  2.2204e-14

Contents of sol.fderval
[1,]        0       -2
[2,]        1       -3



Author: L. Cizkova, 20010506 license MD*Tech
(C) MD*TECH Method and Data Technologies, 05.02.2006