488f4cd8203dc64c9f8485c9c897a9d6d31865c6
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2008-2009 - INRIA - Michael Baudin
4 // Copyright (C) 2011 - DIGITEO - Michael Baudin
5 //
6 //  This file is distributed under the same license as the Scilab package.
7 // =============================================================================
8 // <-- CLI SHELL MODE -->
9 // <-- ENGLISH IMPOSED -->
10 // 1. Test with a scalar argument
11 function y = myfunction (x)
12     y = x*x;
13 endfunction
14 x = 1.0;
15 expected = 2.0;
16 // 1.1 With default parameters
17 computed = numderivative(myfunction, x);
18 assert_checkalmostequal ( computed , expected , 1.e-11 );
19 // 1.2 Test order 1
20 computed = numderivative(myfunction, x, [], 1);
21 assert_checkalmostequal ( computed , expected , 1.e-8 );
22 // 1.3 Test order 2
23 computed = numderivative(myfunction, x, [], 2);
24 assert_checkalmostequal ( computed , expected , 1.e-11 );
25 // 1.4 Test order 4
26 computed = numderivative(myfunction, x, [], 4);
27 assert_checkalmostequal ( computed , expected , 1.e-13 );
28 // 1.5 Compute second numderivative at the same time
29 Jexpected = 2.0;
30 Hexpected = 2.0;
31 [Jcomputed, Hcomputed] = numderivative(myfunction, x);
32 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-11 );
33 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
34 // 1.6 Test order 1
35 [Jcomputed, Hcomputed] = numderivative(myfunction, x, [], 1);
36 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-8 );
37 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-6 );
38 // 1.7 Test order 2
39 [Jcomputed, Hcomputed] = numderivative(myfunction, x, [], 2);
40 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-11 );
41 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
42 // 1.8 Test order 4
43 [Jcomputed, Hcomputed] = numderivative(myfunction, x, [], 4);
44 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-13 );
45 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-11 );
46 // 1.9 Configure the step
47 computed = numderivative(myfunction, x, 1.e-1);
48 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-13 );
49 // 1.10 Configure the step
50 [Jcomputed, Hcomputed] = numderivative(myfunction,x,1.e-1);
51 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-13 );
52 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-11 );
53 // 2. Test with a vector argument
54 function y = myfunction2 (x)
55     y = x(1)*x(1) + x(2) + x(1)*x(2);
56 endfunction
57 x = [1.0; 2.0];
58 Jexpected = [4. 2.];
59 Hexpected = [2. 1. 1. 0.];
60 // 2.1 With default parameters
61 computed = numderivative(myfunction2, x);
62 assert_checkalmostequal ( computed , Jexpected , 1.e-10 );
63 // 2.2 Test order 1
64 computed = numderivative(myfunction2, x, [], 1);
65 assert_checkalmostequal ( computed , Jexpected , 1.e-8 );
66 // 2.3 Test order 2
67 computed = numderivative(myfunction2, x, [], 2);
68 assert_checkalmostequal ( computed , Jexpected , 1.e-10 );
69 // 2.4 Test order 4
70 computed = numderivative(myfunction2, x, [], 4);
71 assert_checkalmostequal ( computed , Jexpected , 1.e-13 );
72 // 2.5 Compute second numderivative at the same time
73 [Jcomputed, Hcomputed] = numderivative(myfunction2, x);
74 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
75 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
76 // 2.6 Test order 1
77 [Jcomputed , Hcomputed] = numderivative(myfunction2, x, [], 1);
78 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-8 );
79 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-5 );
80 // 2.7 Test order 2
81 [Jcomputed, Hcomputed] = numderivative(myfunction2, x, [], 2);
82 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
83 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
84 // 2.8 Test order 4
85 [Jcomputed, Hcomputed] = numderivative(myfunction2, x, [], 4);
86 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-13 );
87 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-10, 1.e-10 );
88 //
89 // 2.9 Configure the step - Expansion of scalar h to the same size as x
90 [Jcomputed, Hcomputed] = numderivative(myfunction2, x, 1.e-1);
91 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-13 );
92 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-10, 1.e-10 );
93 // 2.10 Configure the step - Expansion of scalar h to the same size as x
94 h = %eps^(1/3)*abs(x);
95 [Jcomputed, Hcomputed] = numderivative(myfunction2, x, h);
96 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-8 );
97 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-5 , 1.e-5);
98 // 3. Test H_form
99 // 3.1 Test H_form = "default"
100 Jexpected = [4.0 2.0];
101 Hexpected = [2.0 1.0 1.0 0.0];
102 [Jcomputed, Hcomputed] = numderivative(myfunction2, x, [], [], "default");
103 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
104 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
105 // 3.2 Test H_form = "hypermat"
106 Jexpected = [4.0 2.0];
107 Hexpected = [2.0 1.0
108 1.0 0.0];
109 [Jcomputed, Hcomputed] = numderivative(myfunction2, x , [], [], "hypermat");
110 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
111 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
112 // 3.3 Test H_form = "blockmat"
113 Jexpected = [4.0 2.0];
114 Hexpected = [2.0 1.0
115 1.0 0.0];
116 [Jcomputed, Hcomputed] = numderivative(myfunction2, x, [], [], "blockmat");
117 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
118 assert_checkalmostequal ( Hcomputed , Hexpected , %eps );
119 // 5. Test h parameter
120 // Test a case where the default step h is very small ~ 1.e-9,
121 // but, because the function is very flat in the neighbourhood of the
122 // point, a larger step ~ 1.e-4 reduces the error.
123 // This means that this test cannot pass if the right step is
124 // not taken into account, therefore testing the feature "h is used correctly".
125 myn = 1.e5;
126 function y = myfunction3 (x)
127     y = x^(2/myn);
128 endfunction
129 x = 1.0;
130 h = 6.055454e-006;
131 Jexpected = (2/myn) * x^(2/myn-1);
132 Hexpected = (2/myn) * (2/myn-1) * x^(2/myn-2);
133 [Jcomputed, Hcomputed] = numderivative(myfunction3, x, 1.e-4, 1);
134 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-4 );
135 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-3 );
136 // 6. Test Q parameter
137 function y = myfunction4 (x)
138     y = x(1)*x(1) + x(2)+ x(1)*x(2);
139 endfunction
140 x = [1.; 2.];
141 Jexpected = [4. 2.];
142 Hexpected = [2. 1. 1. 0.];
143 //
144 rand("seed", 0);
145 Q = qr(rand(2, 2));
146 [Jcomputed, Hcomputed] = numderivative(myfunction4, x, [], [], [], Q);
147 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
148 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-8, 1.e-7 );
149 //
150 // 7. Test vector output y
151 function y = myexample(x)
152     f1 = sin(x(1)*x(2)) + exp(x(2)*x(3)+x(1));
153     f2 = sum(x.^3);
154     y = [f1; f2];
155 endfunction
157 function [g1, g2] = exactg(x)
158     g1(1) = cos(x(1)*x(2))*x(2) + exp(x(2)*x(3)+x(1));
159     g1(2) = cos(x(1)*x(2))*x(1) + exp(x(2)*x(3)+x(1))*x(3);
160     g1(3) = exp(x(2)*x(3)+x(1))*x(2);
161     g2(1) = 3*x(1)^2;
162     g2(2) = 3*x(2)^2;
163     g2(3) = 3*x(3)^2;
164 endfunction
165 // The exact Hessian
166 function [H1, H2] = exactH(x)
167     H1(1, 1) = -sin(x(1)*x(2))*x(2)^2 + exp(x(2)*x(3)+x(1));
168     H1(1, 2) = cos(x(1)*x(2)) - sin(x(1)*x(2))*x(2)*x(1) + exp(x(2)*x(3)+x(1))*x(3);
169     H1(1, 3) = exp(x(2)*x(3)+x(1))*x(2);
170     H1(2, 1) = H1(1, 2);
171     H1(2, 2) = -sin(x(1)*x(2))*x(1)^2 + exp(x(2)*x(3)+x(1))*x(3)^2;
172     H1(2, 3) = exp(x(2)*x(3)+x(1)) + exp(x(2)*x(3)+x(1))*x(3)*x(2);
173     H1(3, 1) = H1(1, 3);
174     H1(3, 2) = H1(2, 3);
175     H1(3, 3) = exp(x(2)*x(3)+x(1))*x(2)^2;
176     //
177     H2(1, 1) = 6*x(1);
178     H2(1, 2) = 0;
179     H2(1, 3) = 0;
180     H2(2, 1) = H2(1, 2);
181     H2(2, 2) = 6*x(2);
182     H2(2, 3) = 0;
183     H2(3, 1) = H2(1, 3);
184     H2(3, 2) = H2(2, 3);
185     H2(3, 3) = 6*x(3);
186 endfunction
187 x=[1; 2; 3];
188 [g1, g2] = exactg(x);
189 [H1, H2] = exactH(x);
190 Jexpected = [g1'; g2'];
191 Hexpected = [H1(:)'; H2(:)'];
192 // 7.1.1 Check Jacobian with default options
193 Jcomputed = numderivative(myexample, x);
194 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 );
195 // 7.1.2 Check Jacobian with order = 1
196 Jcomputed = numderivative(myexample, x, [], 1);
197 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-7 );
198 // 7.1.3 Check Jacobian with order = 2
199 Jcomputed = numderivative(myexample, x, [], 2);
200 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 );
201 // 7.1.4 Check Jacobian with order = 4
202 Jcomputed = numderivative(myexample, x, [], 4);
203 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
204 // 7.2.1 Check Jacobian and Hessian with default options
205 [Jcomputed, Hcomputed] = numderivative(myexample, x);
206 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 );
207 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-6 );
208 // 7.2.2 Check Jacobian and Hessian with order = 1
209 [Jcomputed, Hcomputed] = numderivative(myexample, x, [], 1);
210 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-7 );
211 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-4 );
212 // 7.2.3 Check Jacobian and Hessian with order = 2
213 [Jcomputed, Hcomputed] = numderivative(myexample, x, [], 2);
214 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 );
215 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-6 );
216 // 7.2.3 Check Jacobian and Hessian with order = 4
217 [Jcomputed, Hcomputed] = numderivative(myexample, x, [], 4);
218 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-10 );
219 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-4 , 1e-9);
220 //
221 // 7.3 Test with "blockmat"
222 Jexpected = [g1'; g2'];
223 Hexpected = [H1; H2];
224 [Jcomputed, Hcomputed] = numderivative(myexample, x, [], [], "blockmat");
225 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 );
226 assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-6 );
227 //
228 // 7.4 Test with "hypermat"
229 Jexpected = [g1'; g2'];
230 Hexpected = [];
231 Hexpected(:, :, 1) = H1;
232 Hexpected(:, :, 2) = H2;
233 [Jcomputed, Hcomputed] = numderivative(myexample, x, [], [], "hypermat");
234 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 );
235 assert_checkequal ( size(Hcomputed) , [3 3 2] );
236 // This is a limitation of assert (http://bugzilla.scilab.org/show_bug.cgi?id=9461)
237 // assert_checkalmostequal ( Hcomputed , Hexpected , 1.e-6 );
238 assert_checkalmostequal ( Hexpected(:, :, 1) , Hexpected(:, :, 1) , 1.e-6);
239 assert_checkalmostequal ( Hexpected(:, :, 2) , Hexpected(:, :, 2) , 1.e-6);
240 //
241 // 8. Check the number of function evaluations
242 function y = myFevalFun(x)
243     global FEVAL
244     FEVAL = FEVAL + 1;
245     y = sum(x.^3);
246 endfunction
247 n = 3;
248 x = ones(n, 1);
249 //
250 // 8.1 Jacobian with various orders
251 global FEVAL;
252 FEVAL = 0;
253 g = numderivative(myFevalFun, x, [], 1);
254 assert_checkequal ( FEVAL, n+1 );
255 //
256 FEVAL = 0;
257 g = numderivative(myFevalFun, x, [], 2);
258 assert_checkequal ( FEVAL, 2*n );
259 //
260 FEVAL = 0;
261 g = numderivative(myFevalFun, x, [], 4);
262 assert_checkequal ( FEVAL, 4*n );
263 //
264 // 8.2 Hessian with various orders
265 FEVAL = 0;
266 [g, H] = numderivative(myFevalFun, x, [], 1);
267 assert_checkequal ( FEVAL, (n+1)^2+n+1 );
268 //
269 FEVAL = 0;
270 [g, H] = numderivative(myFevalFun, x, [], 2);
271 assert_checkequal ( FEVAL, 4*n^2+2*n );
272 //
273 FEVAL = 0;
274 [g, H] = numderivative(myFevalFun, x, [], 4);
275 assert_checkequal ( FEVAL, 16*n^2+4*n );
276 //
277 // 9. Check error messages.
278 //
279 // 9.1 Cannot evaluate the function - Function case
280 //
281 x = [1.; 2.];
282 Q = qr(rand(2, 2));
283 instr = "J = numderivative(myexample, x, [], [], [], Q)";
284 lclmsg = "%s: Error while evaluating the function: ""%s""\n";
285 assert_checkerror (instr, lclmsg, [], "numderivative", msprintf(_("Invalid index.\n")));
286 // 9.2 Cannot evaluate the function - List case
287 function y = myfunction6 (x, p1, p2)
288     y = p1*x*x + p2;
289 endfunction
290 x = [1.; 2.];
291 Q = qr(rand(2, 2));
292 funf = list(myfunction6, 7., 8.);
293 instr = "J = numderivative(funf, x, [], [], [], Q)";
294 lclmsg = "%s: Error while evaluating the function: ""%s""\n";
295 assert_checkerror (instr, lclmsg, [], "numderivative", msprintf(_("Inconsistent row/column dimensions.\n")));
296 // 9.3 Various error cases
297 x = 2;
298 // Correct syntax: [J, H] = numderivative(myfunction, x)
299 // Number of input arguments
300 instr = "J = numderivative()";
301 lclmsg = "%s: Wrong number of input arguments: %d to %d expected.\n";
302 assert_checkerror (instr, lclmsg, [], "numderivative", 2, 6 );
303 // Wrong type of f
304 myfunction7 = "myfunction";
305 instr = "[J, H] = numderivative(myfunction7, x)";
306 lclmsg = "%s: Wrong type for argument #%d: Function or list expected.\n";
307 assert_checkerror (instr, lclmsg, [], "numderivative", 1 );
308 // Wrong type of x
309 xx = "";
310 instr = "[J, H] = numderivative(myfunction, xx)";
311 lclmsg = "%s: Wrong type for argument #%d: Matrix expected.\n";
312 assert_checkerror (instr, lclmsg, [], "numderivative", 2 );
313 // Wrong type of h
314 hh = "";
315 instr = "[J, H] = numderivative(myfunction, x, hh)";
316 lclmsg = "%s: Wrong type for argument #%d: Matrix expected.\n";
317 assert_checkerror (instr, lclmsg, [], "numderivative", 3 );
318 // Wrong type of order
319 oo = "";
320 instr = "[J, H] = numderivative(myfunction, x, [], oo)";
321 lclmsg = "%s: Wrong type for argument #%d: Matrix expected.\n";
322 assert_checkerror (instr, lclmsg, [], "numderivative", 4 );
323 // Wrong type of H_form
324 HHform = 12;
325 instr = "[J, H] = numderivative(myfunction, x, [], [], HHform)";
326 lclmsg = "%s: Wrong type for input argument #%d: String array expected.\n";
327 assert_checkerror (instr, lclmsg, [], "numderivative", 5 );
328 // Wrong type of Q
329 Q = "";
330 instr = "[J, H] = numderivative(myfunction, x, [], [], [], Q)";
331 lclmsg = "%s: Wrong type for argument #%d: Matrix expected.\n";
332 assert_checkerror (instr, lclmsg, [], "numderivative", 6 );
333 // Wrong size of f (list case)
334 myfunl = list(myfunction);
335 instr = "[J, H] = numderivative(myfunl, x)";
336 lclmsg = "%s: Wrong number of elements in input argument #%d: At least %d elements expected, but current number is %d.\n";
337 assert_checkerror (instr, lclmsg, [], "numderivative", 1, 2, 1 );
338 // Wrong size of x
339 x = [1. 2.; 3. 4.];
340 xx = x';
341 instr = "[J, H] = numderivative(myfunction2, xx)";
342 lclmsg = "%s: Wrong size for input argument #%d: Vector expected.\n";
343 assert_checkerror (instr, lclmsg, [], "numderivative", 2 );
344 // Wrong size of h
345 x = [1.0; 2.0];
346 h = [1; 2; 3];
347 instr = "[J, H] = numderivative(myfunction2, x, h)";
348 lclmsg = "%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n";
349 assert_checkerror (instr, lclmsg, [], "numderivative", 3, 1);
350 // Wrong size of order
351 x = [1.0; 2.0];
352 order = [1; 2; 4];
353 instr = "[J, H] = numderivative(myfunction2, x, [], order)";
354 lclmsg = "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n";
355 assert_checkerror (instr, lclmsg, [], "numderivative", 4, 1, 1);
356 // Wrong size of H_form
357 x = [1.0; 2.0];
358 H_form = ["blockmat" "hypermat"];
359 instr = "[J, H] = numderivative(myfunction2, x, [], [], H_form)";
360 lclmsg = "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n";
361 assert_checkerror (instr, lclmsg, [], "numderivative", 5, 1, 1);
362 // Wrong size of Q
363 x = [1.0; 2.0];
364 Q = ones(2, 3);
365 instr = "[J, H] = numderivative(myfunction2, x, [], [], [], Q)";
366 lclmsg = "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n";
367 assert_checkerror (instr, lclmsg, [], "numderivative", 6, 2, 2);
368 // 10. Check that a nonzero step is used for components of x which are zero.
369 // Check also that a scaled step is used.
370 x=[0; 0; 1.e-50];
371 [g1, g2] = exactg(x);
372 [H1, H2] = exactH(x);
373 Jexpected = [g1'; g2'];
374 Hexpected = [H1(:)'; H2(:)'];
375 // 10.1 Check Jacobian and Hessian with default options
376 Jcomputed = numderivative(myexample, x);
377 assert_checkalmostequal ( Jcomputed , Jexpected , 1.e-9 , 1.e-10);
378 // For the Jacobian, the step used is h = [%eps^(1/3); %eps^(1/3); %eps^(1/3)*1.e50]
379 // 11. Prove the numerical superiority of numderivative.
380 // Although the step provided by numderivative is not always
381 // optimal, it is often sufficiently accurate.
382 // Print the number of significant digits.
383 // x = 1.00D-100, numdiff = 0, derivative = 0,  numderivative = 11
384 // x = 1.000D-32, numdiff = 0, derivative = 0,  numderivative = 10
385 // x = 1.000D+32, numdiff = 5, derivative = 0,  numderivative = 11
386 // x = 1.00D+100, numdiff = 5, derivative = 0,  numderivative = 10
387 // x = 1,         numdiff = 7, derivative = 11, numderivative = 11
388 function y = myfunction10 (x)
389     y = x^3;
390 endfunction
391 function y = g10 (x)
392     y = 3*x^2;
393 endfunction
394 for x = [10^-100 10^-32 10^32 10^100 1]
395     exact = g10 (x);
396     g3 = numderivative(myfunction10, x);
397     d3 = assert_computedigits ( g3, exact );
398     assert_checktrue ( d3 > 9 );
399 end
400 // 12. Check that numderivative also accepts row vector x
401 function f = myfunction11(x)
402     f = x(1)*x(1) + x(1)*x(2);
403 endfunction
404 // 12.1 Check with row x
405 x = [5 8];
406 g = numderivative(myfunction11, x);
407 exact = [18 5];
408 assert_checkalmostequal ( g , exact , 1.e-9 );
409 // 12.2 Check with h
410 h = sqrt(%eps)*(1+1d-3*abs(x));
411 g = numderivative(myfunction11, x, h);
412 assert_checkalmostequal ( g , exact , 1.e-8 );
413 // 13. Check that we can derivate a compiled function
414 x = 1;
415 g = numderivative (sqrt, x);
416 assert_checkalmostequal ( g , 0.5 , 1.e-8 );
417 // 14.1 Check that numderivative works when f takes extra arguments
418 function y = f(x, A, p, w)
419     y = x'*A*x + p'*x + w;
420 endfunction
421 // with Jacobian and Hessian given
422 // by J(x) = x'*(A+A')+p' and H(x) = A+A'.
423 A = rand(3, 3);
424 p = rand(3, 1);
425 w = 1;
426 x = rand(3, 1);
427 h = 1;
428 [J, H] = numderivative(list(f, A, p, w), x, h, [], "blockmat");
429 assert_checkalmostequal( J , x'*(A+A')+p' );
430 assert_checkalmostequal( H , A+A' );
431 // 14.2 Same test with a different function
432 function y = G(x, p)
433     f1 = sin(x(1)*x(2)*p) + exp(x(2)*x(3)+x(1));
434     f2 = sum(x.^3);
435     y = [f1; f2];
436 endfunction
437 x = rand(3, 1);
438 p = 1;
439 df1_dx1 = cos(x(1)*x(2)*p)*x(2)*p+exp(x(2)*x(3)+x(1));
440 df1_dx2 = cos(x(1)*x(2)*p)*x(1)*p+exp(x(2)*x(3)+x(1))*x(3);
441 df1_dx3 = exp(x(2)*x(3)+x(1))*x(2);
442 df2_dx1 = 3*x(1)^2;
443 df2_dx2 = 3*x(2)^2;
444 df2_dx3 = 3*x(3)^2;
445 expectedJ = [df1_dx1   df1_dx2   df1_dx3;
446 df2_dx1   df2_dx2  df2_dx3 ];
447 df1_dx11 = -sin(x(1)*x(2)*p)*x(2)^2*p^2+exp(x(2)*x(3)+x(1));
448 df1_dx12 = -sin(x(1)*x(2)*p)*x(1)*x(2)*p^2+cos(x(1)*x(2)*p)*p+exp(x(2)*x(3)+x(1))*x(3);
449 df1_dx13 = exp(x(2)*x(3)+x(1))*x(2);
450 df1_dx21 = df1_dx12;
451 df1_dx22 = -sin(x(1)*x(2)*p)*x(1)^2*p^2+exp(x(2)*x(3)+x(1))*x(3)^2;
452 df1_dx23 = exp(x(2)*x(3)+x(1))*(1+x(3)*x(2));
453 df1_dx31 = df1_dx13;
454 df1_dx32 = df1_dx23;
455 df1_dx33 = exp(x(2)*x(3)+x(1))*x(2)^2;
456 df2_dx11 = 6*x(1);
457 df2_dx12 = 0;
458 df2_dx13 = 0;
459 df2_dx21 = 0;
460 df2_dx22 = 6*x(2);
461 df2_dx23 = 0;
462 df2_dx31 = 0;
463 df2_dx32 = 0;
464 df2_dx33 = 6*x(3);
465 expectedH = [df1_dx11   df1_dx12   df1_dx13   df1_dx21   df1_dx22   df1_dx23   df1_dx31   df1_dx32   df1_dx33;
466 df2_dx11   df2_dx12   df2_dx13   df2_dx21   df2_dx22   df2_dx23   df2_dx31   df2_dx32   df2_dx33 ];
467 [J, H] = numderivative(list(G, p), x);
468 assert_checkalmostequal( J , expectedJ );
469 assert_checkalmostequal( H , expectedH , [], 1e-7);