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