128 lines
3.0 KiB
Scilab
128 lines
3.0 KiB
Scilab
function [J]=diffext(f,x,options,varargin)
|
|
[nargout,nargin] = argn(0)
|
|
//---------------------------------------------------------------------------
|
|
//DIFFEXT Numerical approximation for hessian.
|
|
// The method is Richardson`s extrapolation.
|
|
// Sample call
|
|
// [D,err,relerr,n] = diffext('f',x,delta,toler)
|
|
// Inputs
|
|
// f name of the function
|
|
// x differentiation point
|
|
// options matrix of algorithm parameters
|
|
// delta error goal (1e-12) (suppressed MJ 02/27/02)
|
|
// toler relative error goal (1e-12)
|
|
// Return
|
|
// J Jacobian
|
|
//
|
|
// NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1995
|
|
//
|
|
// Modified F. Collard, August 2001
|
|
//---------------------------------------------------------------------------
|
|
if nargin > 2 then
|
|
if ~(options==[]) then
|
|
// %delta = options(1);
|
|
toler = options(2);
|
|
maxit = options(3);
|
|
else
|
|
// %delta = 1e-12;
|
|
toler = 1e-12;
|
|
maxit = 20;
|
|
end
|
|
else
|
|
%delta = 1e-12;
|
|
toler = 1e-12;
|
|
maxit = 20;
|
|
end
|
|
con = 1.4;
|
|
con2 = con*con;
|
|
big = 1e30;
|
|
safe = 2;
|
|
if nargin > 3
|
|
ff = evstr(f+'(x,varargin)');
|
|
else
|
|
ff = evstr(f+'(x)');
|
|
end
|
|
nx = size(x,1);
|
|
nf = size(ff,1);
|
|
J = zeros(nf,nx);
|
|
for xi = 1:nx
|
|
err = big*ones(nf,1);;
|
|
// relerr = big*ones(nf,1);
|
|
h = max(abs(x(xi))/10, 10*gstep_)*100*gstep_;
|
|
dx = zeros(nx,1);
|
|
dx(xi,1) = h;
|
|
if nargin > 3
|
|
fs = evstr(f+'(x+dx,varargin)');
|
|
fm = evstr(f+'(x-dx,varargin)');
|
|
else
|
|
fs = evstr(f+'(x+dx)');
|
|
fm = evstr(f+'(x-dx)');
|
|
end
|
|
D1 = (fs-fm)/(2*h);
|
|
mask = ones(nf,1);
|
|
j = 2;
|
|
while (1)
|
|
D2 = zeros(nf,j);
|
|
h = h/con;
|
|
dx(xi,1) = h;
|
|
if nargin > 3
|
|
fs = evstr(f+'(x+dx,varargin)');
|
|
fm = evstr(f+'(x-dx,varargin)');
|
|
else
|
|
fs = evstr(f+'(x+dx)');
|
|
fm = evstr(f+'(x-dx)');
|
|
end
|
|
for fi = 1:nf
|
|
if mask(fi)
|
|
err2 = big;
|
|
D2(fi,1) = (fs(fi)-fm(fi))/(2*h);
|
|
fac = con2;
|
|
for k = 2:j
|
|
D2(fi,k) = D2(fi,k-1)+(D2(fi,k-1)-D1(fi,k-1))/(fac-1);
|
|
fac = con2*fac;
|
|
errt = max(abs(D2(fi,k)-D2(fi,k-1)),abs(D2(fi,k)-D1(fi,k-1)));
|
|
if errt <= err2
|
|
err2 = errt;
|
|
deriv = D2(fi,k);
|
|
end
|
|
end
|
|
err(fi) = abs(D2(fi,j)-D1(fi,j-1));
|
|
// relerr(fi) = 2*err(fi)/(abs(D2(fi,j))+abs(D1(fi,j-1))+%eps);
|
|
// if (err(fi) < toler & relerr(fi) < %delta)| err(fi) > safe*err2
|
|
if err(fi) < toler | err(fi) > safe*err2
|
|
J(fi,xi) = deriv;
|
|
mask(fi) = 0;
|
|
end
|
|
end
|
|
end
|
|
if (mask == 0) then
|
|
break
|
|
end
|
|
j = j+1;
|
|
if j == maxit
|
|
error('DIFFEXT didn''t converge. Try to increase gstep_ (default 0.01)')
|
|
end
|
|
D1 = D2;
|
|
end
|
|
[m_err,i] = max(err);
|
|
if m_err > 1e-12
|
|
error('DIFFEXT obtains an accuracy > 1e-12. Try to increase gstep_ (default 0.01)')
|
|
end
|
|
// [m_err,i] = max(relerr);
|
|
// if m_err > 1e-12
|
|
// dyn_disp(err)
|
|
// dyn_disp(relerr)
|
|
// dyn_disp(D2)
|
|
// error('DIFFEXT obtains an relative accuracy > 1e-12. Try to increase gstep_ (default 0.01)')
|
|
// end
|
|
end
|
|
|
|
// 10/12/2001 MJ modified initial h
|
|
// 02/25/2002 MJ put equation look inside
|
|
|
|
|
|
|
|
|
|
|
|
|