|
0
|
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 |
