ndiff

(tmp_for_tar/misc/ndiff.m)

Function Synopsis

df = ndiff( func, ... )

Help text

       df = ndiff( func [, options], argument(s) )

 Numerical differentiation of 'func', by symmetric (default) or
 assymetric difference scheme.

 By default, differentiation is with respect to the first argument
 found after the eventual options.

 If this variable (call it 'x') is of size RxC, and the output 'y' of
 'func' is SxD, then 'df' will be (S*D)x(R*C).

                  df(i,j) = dy(i)/dx(j)

 For quicker differentiation, see cdiff and compare_speed_ncdiff. 

 OPTIONS :
 ---------
 "dx" , small_value : Value of the small perturbation.
                                                      Default = 1e-6 

 "wrt", variable    : With Respect To : Name or number, in argument
                      list, of the differentiation variable. "varnum"
                      can be used instead of "wrt". 
                      Octave>=2.1.XX is needed for specifying a
                      variable by its name. 
                                                          Default = 1 

 "rstack"  : allow stacking of arguments in rows. This means that

             [ func(x1) ; func(x2) ] == func( [ x1 ; x2 ] )

             (for a one-arg function). When this option is passed,
             'func' is called once only, so this may provide some
             speed-up. 

 "cstack"  : allow stacking of arguments in columns. This means that

             [ func(x1) , func(x2) ] == func( [ x1 , x2 ] )

             (for a one-arg function). When this option is passed,
             'func' is called once only. 

 "nostack" : is ignored.

 "assym"   : Perform calculation using an assymetric finite difference
             scheme rather than the symmetric scheme. With this
             option, 'func' is evaluated at R*C+1 points rather than
             the 2*R*C that are necessary for the symmetric scheme.
             Beware that error is then linear in the value of 'dx'
             rather than quadratic.

 "args"    : When 'ndiff' finds the "args" option, it stops reading
             options. All remaining arguments are passed to 'func'.
             Without "args", 'ndiff' stops after reading the first
             argument that it does not understand. This arg and the
             following are passed to 'func'.

 EXAMPLES :
 ----------
 function y = foo(x,y,z) ... end

 df = ndiff("foo","wrt","y0",x0,y0,z0)

 ## Returns an approximation of the differential of 'foo' wrt. 2nd
 ## variable, at (x0,y0,z0). 

 df = ndiff("foo",x0,y0,z0)

 ## Returns an approximation of the differential of 'foo' wrt. 'x', at
 ## (x0,y0,z0). 

 Last modified: April 2001

Listing of function file tmp_for_tar/misc/ndiff.m

##       df = ndiff( func [, options], argument(s) )
##
## Numerical differentiation of 'func', by symmetric (default) or
## assymetric difference scheme.
##
## By default, differentiation is with respect to the first argument
## found after the eventual options.
##
## If this variable (call it 'x') is of size RxC, and the output 'y' of
## 'func' is SxD, then 'df' will be (S*D)x(R*C).
##
##                  df(i,j) = dy(i)/dx(j)
##
## For quicker differentiation, see cdiff and compare_speed_ncdiff. 
## 
## OPTIONS :
## ---------
## "dx" , small_value : Value of the small perturbation.
##                                                      Default = 1e-6 
##
## "wrt", variable    : With Respect To : Name or number, in argument
##                      list, of the differentiation variable. "varnum"
##                      can be used instead of "wrt". 
##                      Octave>=2.1.XX is needed for specifying a
##                      variable by its name. 
##                                                          Default = 1 
##
## "rstack"  : allow stacking of arguments in rows. This means that
##
##             [ func(x1) ; func(x2) ] == func( [ x1 ; x2 ] )
##
##             (for a one-arg function). When this option is passed,
##             'func' is called once only, so this may provide some
##             speed-up. 
##
## "cstack"  : allow stacking of arguments in columns. This means that
##
##             [ func(x1) , func(x2) ] == func( [ x1 , x2 ] )
##
##             (for a one-arg function). When this option is passed,
##             'func' is called once only. 
##             
## "nostack" : is ignored.
## 
## "assym"   : Perform calculation using an assymetric finite difference
##             scheme rather than the symmetric scheme. With this
##             option, 'func' is evaluated at R*C+1 points rather than
##             the 2*R*C that are necessary for the symmetric scheme.
##             Beware that error is then linear in the value of 'dx'
##             rather than quadratic.
## 
## "args"    : When 'ndiff' finds the "args" option, it stops reading
##             options. All remaining arguments are passed to 'func'.
##             Without "args", 'ndiff' stops after reading the first
##             argument that it does not understand. This arg and the
##             following are passed to 'func'.
## 
## EXAMPLES :
## ----------
## function y = foo(x,y,z) ... end
## 
## df = ndiff("foo","wrt","y0",x0,y0,z0)
##
## ## Returns an approximation of the differential of 'foo' wrt. 2nd
## ## variable, at (x0,y0,z0). 
##
## df = ndiff("foo",x0,y0,z0)
##
## ## Returns an approximation of the differential of 'foo' wrt. 'x', at
## ## (x0,y0,z0). 
 

## Author:        Etienne Grossmann  <etienne@isr.ist.utl.pt>
## Last modified: April 2001

function df = ndiff( func, ... )


verbose = 0;
prudent = 1 ;
assym = 0;
rstack = cstack = 0 ;

dx = 1e-6 ;

filename = "ndiff" ;
opt0 = " verbose assym rstack nostack cstack " ;
opt1 = " dx wrt varnum " ;
nargin -= 1 ;
nargin_orig = nargin ;
opt_quiet = 1 ;

read_options ;

nvars = nargin + 1 ;		# Number of args to func

vskip = nargin_orig - nvars ;	# 1st arg to func on arg list

if isstr(last_option_name) && strcmp(last_option_name,"args"),
  nvars-- ;
  vskip++ ;
end
				# Default varnum
if exist("varnum","var")!=1, 
  if exist("wrt","var")==1, varnum = wrt ;
  else                varnum = 1 ;
  end
end

if isstr(varnum),		# find variable's number
  varname = varnum ;
  ## printf("Looking for argument called '%s'!\n",varnum) ;
  for cnt = 1:vskip+nvars+1,
    aname =  sprintf("%s",argn(cnt,:));
    ## printf("arg %d is called '%s'\n",cnt,aname);
    if strcmp(aname,varname), 
      ## printf("found!\n") ;
      if isstr(varnum),
	varnum = cnt-vskip-1 ;
      else
	printf(["ndiff wrt '%s' Whoa!! multiple args w/ that name\n",\
		"      will use arg number %d\n"], varname,varnum) ;
	break 
      end
    end
  end
end


sayif(verbose,strcat("ndiff : differentiating          func='%s'\n",\
		     "        number of variables is  nvars= %d\n",\
		     "        diff is wrt var number varnum= %d\n",\
		     "        assymetric difference?  assym= %d\n",\
		     "        step size                  dx= %g\n"),
      func,nvars,varnum,assym,dx);

ckcnt = 1;			# Check
for i= [ !isstr(varnum), varnum>0, !rstack || !cstack ],
  if !i,
    printf("ndiff : check %d fails\n",ckcnt) ;
    keyboard ;
  end
  ckcnt++ ;
end


				# Second loop thru args
va_start() ;			# skip options
cnt = 1 ;
while cnt++ <= vskip, va_arg(); end

				# read the other args
varstr = sprintf("var%d=va_arg();",1:nvars) ;
varstr = strrep(varstr,sprintf("var%d",varnum),"dvar") ;
sayif(verbose,"ndiff : reading args as ---\n%s\n--\n",varstr) ;
eval(varstr) ;
dvar0 = dvar ;

				# Function call string
## ",var1,var2,...,varN"
funcstr = sprintf(",var%d",1:nvars) ;
sayif(verbose,"ndiff : funcstr=%s\n",funcstr) ;

## ",var1,dvar,...,varN"
funcstr = strrep(funcstr,sprintf("var%d",varnum),"dvar") ;
sayif(verbose,"ndiff : funcstr=%s\n",funcstr) ;

## "func(var1,dvar,...,varN);"
funcstr = sprintf("%s(%s);",func,funcstr(2:length(funcstr)));
sayif(verbose,"ndiff : funcstr=%s\n",funcstr) ;

sz = size(dvar0) ;
nv = prod(sz) ;			# Number of scalars in dvar

				# Points of evaluation
## if rstack, dvar = dvar' ; end
  
if assym,  
  ## keyboard
  dvar = [zeros(nv,1), dx*eye(nv) ] + kron(dvar(:),ones(1,nv+1));
  np = nv+1 ;
else         
  dvar = [-dx*eye(nv), dx*eye(nv) ] + kron(dvar(:),ones(1,2*nv));
  np = 2*nv ;
end

				# Evaluate function

				# rsz will be size of one point's image
if cstack, 
  dvar = reshape(dvar,sz(1),sz(2)*np);
  result = eval(funcstr);	# rsz(1) x (rsz(2)*np)
  rsz = [size(result) ./ [1,np]] ;
  ## keyboard
  ixy = [ones(1,np)              ;\ # start row
	 rsz(1)*ones(1,np)       ;\ # end row
	 1:rsz(2):np*rsz(2)      ;\ # start col
	 rsz(2):rsz(2):np*rsz(2) ]; # end col

elseif rstack,
  scramble = reshape(1:nv,sz(1),sz(2))' ;
  dvar = reshape(dvar(scramble,:),sz(2),sz(1)*np)';
  result = eval(funcstr);	# (rsz(1)*np) x (rsz(2))
  rsz = [size(result) ./ [np,1]] ;
  ## keyboard

  ixy = [ 1:rsz(1):np*rsz(1)     ;\ # start row
	 rsz(1):rsz(1):np*rsz(1) ;\ # end row
	 ones(1,np)              ;\ # start col
	 rsz(2)*ones(1,np)       ]; # end col


else
				# Must do 'np' function calls. Massage
				# funcstr to do a single eval().

  ## "func(var1,dvar,...,varN),"
  funcstr = strrep(funcstr,";",",") ;
  ## sayif(verbose,"ndiff : funcstr='%s'\n",funcstr) ;

  ## "func(var1,reshape(dvar(%d:%d,%d),SZ1,SZ2),...,varN),"
  funcstr = strrep(funcstr,"dvar",                                    \
		   ["reshape(dvar(%d:%d,%d)",sprintf(",%d,%d)",sz)]);
  ## sayif(verbose,"ndiff : funcstr='%s'\n",funcstr) ;

  ## "func(var1,reshape(dvar(1:NV,1),SZ1,SZ2),...,varN)," , 
  ## "func(var1,reshape(dvar(1:NV,2),SZ1,SZ2),...,varN)," , ...
  ## "func(var1,reshape(dvar(1:NV,NP),SZ1,SZ2),...,varN)"
  funcstr = sprintf(funcstr,[ones(1,np);nv*ones(1,np);1:np]) ;
  ## sayif(verbose,"ndiff : funcstr='%s'\n",funcstr) ;
  funcstr = ["[",funcstr(1:length(funcstr)-1),"];"] ;
  sayif(verbose,"ndiff : funcstr='%s'\n",funcstr) ;
  ## keyboard
  result = eval(funcstr) ;
  rsz = [size(result) ./ [1,np]] ;
  ## keyboard
  ixy = [ones(1,np)              ;\ # start row
	 rsz(1)*ones(1,np)       ;\ # end row
	 1:rsz(2):np*rsz(2)      ;\ # start col
	 rsz(2):rsz(2):np*rsz(2) ]; # end col

end

				# Compute derivatives
if assym,
  ipt = [ones(1,nv) ; 2:nv+1    ] ;
  dx = 0.5*dx ;
else 
  ipt = [      1:nv ; nv+1:2*nv ] ;
end
## keyboard
df = zeros(prod(rsz),nv) ;
for i = 1:nv,
  df(:,i) = [result(ixy(1,ipt(2,i)):ixy(2,ipt(2,i)),   \
		    ixy(3,ipt(2,i)):ixy(4,ipt(2,i))) - \
	     result(ixy(1,ipt(1,i)):ixy(2,ipt(1,i)),   \
		    ixy(3,ipt(1,i)):ixy(4,ipt(1,i))) ](:) ;
end
df = df ./(2*dx) ;