Mercurial > forge

comparison main/optim/dfdp.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 2ac2777b30bc
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 % Copyright (C) 1992-1994 Richard Shrager, Arthur Jutan and Ray Muzic
2 %
3 % This program is free software; you can redistribute it and/or modify
4 % it under the terms of the GNU General Public License as published by
5 % the Free Software Foundation; either version 2 of the License, or
6 % (at your option) any later version.
7 %
8 % This program is distributed in the hope that it will be useful,
9 % but WITHOUT ANY WARRANTY; without even the implied warranty of
10 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 % GNU General Public License for more details.
12 %
13 % You should have received a copy of the GNU General Public License
14 % along with this program; if not, write to the Free Software
15 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 function prt=dfdp(x,f,p,dp,func)
18 % numerical partial derivatives (Jacobian) df/dp for use with leasqr
19 % --------INPUT VARIABLES---------
20 % x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
21 % f=func(x,p) vector initialsed by user before each call to dfdp
22 % p= vec of current parameter values
23 % dp= fractional increment of p for numerical derivatives
24 % dp(j)>0 central differences calculated
25 % dp(j)<0 one sided differences calculated
26 % dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
27 % func=string naming the function (.m) file
28 % e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum')
29 %----------OUTPUT VARIABLES-------
30 % prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
31 %================================
32
33 m=size(x,1); if (m==1), m=size(x,2); end %# PAK: in case #cols > #rows
34 n=length(p); %dimensions
35 ps=p; prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero
36 for j=1:n
37 del(j)=dp(j) .*p(j); %cal delx=fract(dp)*param value(p)
38 if p(j)==0
39 del(j)=dp(j); %if param=0 delx=fraction
40 end
41 p(j)=ps(j) + del(j);
42 if del(j)~=0, f1=feval(func,x,p);
43 if dp(j) < 0, prt(:,j)=(f1-f)./del(j);
44 else
45 p(j)=ps(j)- del(j);
46 prt(:,j)=(f1-feval(func,x,p))./(2 .*del(j));
47 end
48 end
49 p(j)=ps(j); %restore p(j)
50 end
51 return