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