1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num. differentiation
17 ## = "function df = dfunc (var1,..,dvar,..,varN) .. endfunction
19 ## Returns a string of octave code that defines a function 'dfunc' that
20 ## returns the derivative of 'func' with respect to it's 'wrt'th
23 ## The derivatives are obtained by symmetric finite difference.
25 ## dfunc()'s return value is in the same format as that of ndiff()
27 ## func : string : name of the function to differentiate
29 ## wrt : int : position, in argument list, of the differentiation
30 ## variable. Default:1
32 ## N : int : total number of arguments taken by 'func'.
33 ## If N=inf, dfunc will take variable argument list.
36 ## dfunc : string : Name of the octave function that returns the
37 ## derivatives. Default:['d',func]
39 ## stack : string : Indicates whether 'func' accepts vertically
40 ## (stack="rstack") or horizontally (stack="cstack")
41 ## arguments. Any other string indicates that 'func'
42 ## does not allow stacking. Default:''
44 ## dx : real : Step used in the symmetric difference scheme.
45 ## Default:10*sqrt(eps)
47 ## See also : ndiff, eval, todisk
49 function c = cdiff (func,wrt,nargs,dfunc,stack,dx)
57 if nargin<4 || strcmp(dfunc,""),
60 printf(["cdiff : Warning : name of derivative not specified\n",\
61 " and canonic name '%s' is already taken\n"],\
66 if nargin<5, stack = "" ; end
67 if nargin<6, dx = 10*sqrt(eps) ; end
70 ## build argstr = "var1,..,dvar,...var_nargs"
72 argstr = sprintf("var%i,",1:nargs);
74 argstr = [sprintf("var%i,",1:wrt),"...,"];
77 argstr = strrep(argstr,sprintf("var%i",wrt),"dvar") ;
78 argstr = argstr(1:length(argstr)-1) ;
80 if strcmp("cstack",stack) , # Horizontal stacking ################
82 calstr = "reshape (kron(ones(1,2*ps), dvar(:))+[-dx*eye(ps),dx*eye(ps)], sz.*[1,2*ps])";
83 calstr = strrep(argstr,"dvar",calstr) ;
84 calstr = sprintf("%s(%s)",func,calstr) ;
86 calstr = sprintf(strcat(" res = %s;\n",
87 " pr = prod (size(res)) / (2*ps);\n",
88 " res = reshape (res,pr,2*ps);\n",
89 " df = (res(:,ps+1:2*ps)-res(:,1:ps)) / (2*dx);\n"),
93 elseif strcmp("rstack",stack), # Vertical stacking ##################
95 calstr = "kron(ones(2*ps,1),dvar)+dx*[-dv;dv]" ;
96 calstr = strrep(argstr,"dvar",calstr) ;
97 calstr = sprintf("%s(%s)",func,calstr) ;
99 calstr = sprintf(strcat(" dv = kron (eye(sz(2)), eye(sz(1))(:));\n",\
101 " sr = size(res)./[2*ps,1];\n",\
102 " pr = prod (sr);\n",\
103 " df = (res(sr(1)*ps+1:2*sr(1)*ps,:)-res(1:sr(1)*ps,:))/(2*dx);\n",\
104 " scramble = reshape (1:pr,sr(2),sr(1))';\n",\
105 " df = reshape (df',pr,ps)(scramble(:),:);\n"),\
107 ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
108 else # No stacking ########################
109 calstr = sprintf("%s (%s)",func,argstr) ;
110 ## "func(var1,dvar%sdv(:,%d:%d),...,varN),"
111 ## calstr = strrep(calstr,"dvar","dvar%sdv(:,(i-1)*sz(2)+1:i*sz(2))")(:)';
113 calstr = strrep(calstr,"dvar","dvar%sdv")(:)';
115 ## func(..,dvar+dv(:,1:sz(2)),..) - func(..)
116 calstr = strcat(calstr,"-",calstr) ; ## strcat(calstr,"-",calstr) ;
117 calstr = sprintf(calstr,"+","-") ;
119 ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
120 calstr = sprintf(strcat(" dv = zeros (sz); dv(1) = dx;\n",\
122 " sr = size (df0);\n",\
123 " df = zeros(prod (sr),ps); df(:,1) = df0(:);\n",\
125 " dv(i) = dx; dv(i-1) = 0;\n",\
126 " df(:,i) = (%s)(:);\n",\
133 ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
135 ## "func(var1,reshape(dvar(1:NV,1),SZ1,SZ2),...,varN)," ,
136 ## "func(var1,reshape(dvar(1:NV,2),SZ1,SZ2),...,varN)," , ...
137 ## "func(var1,reshape(dvar(1:NV,NP),SZ1,SZ2),...,varN)"
138 ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
140 argstr = strrep (argstr, "...", "varargin");
141 calstr = strrep (calstr, "...", "varargin{:}");
143 c = sprintf(strcat("function df = %s (%s)\n",\
144 " ## Numerical differentiation of '%s' wrt to it's %d'th argument\n",\
145 " ## This function has been written by 'cdiff()'\n",\
147 " sz = size (dvar);\n",\
148 " ps = prod (sz);\n",\