gradp.m (2783B)
1 % GRADP.ARC 2 % (C) Copyright 1988-1991 by Aptech Systems, Inc. 3 % All Rights Reserved. 4 % 5 % Purpose: Computes the gradient vector or matrix (Jacobian) of a 6 % vector-valued function that has been defined in a procedure. 7 % Single-sided (forward difference) gradients are computed. 8 % 9 % Format: g = gradp(&f,x0); 10 % 11 % Input: f -- scalar, procedure pointer to a vector-valued function: 12 % 13 % f:Kx1 -> Nx1 14 % 15 % It is acceptable for f(x) to have been defined in terms of 16 % global arguments in addition to x, and thus f can return 17 % an Nx1 vector: 18 % 19 % proc f(x); 20 % retp( exp(x*b) ); 21 % endp; 22 % 23 % x0 -- Kx1 vector of points at which to compute gradient. 24 % 25 % Output: g -- NxK matrix containing the gradients of f with respect 26 % to the variable x at x0. 27 % 28 % Remarks: GRADP will return a ROW for every row that is returned by f. 29 % For instance, if f returns a 1x1 result, then GRADP will 30 % return a 1xK row vector. This allows the same function to be used 31 % where N is the number of rows in the result returned by f. 32 % Thus, for instance, GRADP can be used to compute the 33 % Jacobian matrix of a set of equations. 34 % 35 % Example: proc myfunc(x); 36 % retp( x .* 2 .* exp( x .* x ./ 3 )); 37 % endp; 38 % 39 % x0 = { 2.5, 3.0, 3.5 }; 40 % y = gradp(&myfunc,x0); 41 % 42 % 82.98901842 0.00000000 0.00000000 43 % y = 0.00000000 281.19752975 0.00000000 44 % 0.00000000 0.00000000 1087.95414117 45 % 46 % It is a 3x3 matrix because we are passing it 3 arguments and 47 % myfunc returns 3 results when we do that. The off-diagonals 48 % are zeros because the cross-derivatives of 3 arguments are 0. 49 % 50 % Globals: None 51 52 function grdd=gradp(f,x0,preci,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10); 53 54 if isempty(preci); preci=1e-8; end; 55 evalstr = ['f0 =' f '(x0']; 56 pstr=[]; 57 for i=1:nargin - 3 58 pstr = [pstr,',P',int2str(i)]; 59 end 60 evalstr = [evalstr pstr ');']; 61 62 eval(evalstr); 63 n = rows(f0); 64 k = rows(x0); 65 grdd = zeros(n,k); 66 67 % Computation of stepsize (dh) for gradient 68 69 ax0 = abs(x0); 70 if all(x0~=0); 71 dax0 = div(x0,ax0); 72 else; 73 dax0 = ones(k,1); 74 end; 75 dh = preci*mult(max([ax0 (1e-2)*ones(k,1)]')',dax0); 76 xdh = x0+dh; 77 dh = xdh-x0; % This increases precision slightly 78 79 arg = kron(x0,ones(1,k)) + diag(dh); 80 for i=1:k; 81 eval(['grdd(:,i)=' f '(arg(:,i)' pstr ')-f0;']); 82 end; 83 grdd = div(grdd,dh'); 84