Example 1 (Magnetism data)
We apply the holonomic gradient descent method to the problem
 minimize f(x,y,r=1) * exp( - S11*x11 - S12*x12 - S13*x13 - S22*x22 -
          S23*x23 - S33*x33 - S1*y1 - S2*y2 - S3*y3 )
 subject to (x11,x12,x13,x22,x23,x33,y1,y2,y3) in A (in R^9)
where 
  S11 = 0.045,
  S12 = -0.075,
  S13 = 0.014,
  S22 = 0.921,
  S23 = -0.122,
  S33 = 0.034,
  S1 = 0.083,
  S2 = -0.959,
  S3 = -0.131,
and
  A=domain([-30,-30,-30,-30, -30,  -1.0, -0.1,  -30,  -30],
           [ 10, 10, 10, 10,  20,  -0.1,  1.0, -0.1, 10])
   =[-30,10] x [-30,10] x ... x [-30,-0.1] x [-30,10]
The function f(x,y,r) is the Fisher-Bingham integral on the sphere S^2.


---------------------------
We provide some programs to solve this problem on the computer algebra
system Risa/Asir and Maple 
at http://www.math.kobe-u.ac.jp/OpenXM/Math/Fisher-Bingham
(Programs are provided as they are in our study and keep our tries and have
 no cosmetics. 
A user friendly package of the holonomic gradient descent method will be provided later.)

Instruction to execute this program (on Debian GNU linux)  
1. Download all programs with the extension .rr
 (sei15c.rr, sei15.rr set15a.rr, sei13.rr)
2. If you do not have Risa/Asir (version >= 20100206), install it:
  2.a. Add the line 
        http://fe.math.kobe-u.ac.jp/KnoppixMath ./
       to /etc/apt/sources.list
  2.b. sudo su
  2.c  apt-get update
  2.d  dpkg --purge openxm
  2.e. apt-get install openxm
 You do not need Maple to execute the function test5b_demo().
3. Start the Risa/Asir, load sei15c.rr, and  execute test5b_demo():
   asir
   load("sei15c.rr");
   test5b_demo();

The starting point of the holonomic gradient descent method is taken at
 Pos=[0, 0, 0, 5.5, -20.0, -0.2, 0.5, -15.5, -17.5, 1].
This starting point is found by a numerical integration of the target function
on a coarse grid of the size 0.5 (see mhlml() in sei14c.ml).

During the first execution, the big data file "tmp-sei14c-r1.ab" is 
automatically downloaded by "wget".  It will take a while.
The file "tmp-sei14c-r1.ab" is the Pfaffian system for the Fisher-Bigham
integral on S^2.

The output of this function will be

[1194,0,[7.05,-2.05,9.95,6.45,-29.95,-0.95,0.95,-29.95,-29.95],[0.0003236947435145229824,0.01365914116495716502,0.0003193051763246178940,0.009267060184740950526,0.00002355671589603190298,-0.00001046867630351169487]]
point = [[x11,7.05],[x12,-2.05],[x13,9.95],[x22,6.45],[x23,-29.95],[x33,-0.95],[y1,0.95],[y2,-29.95],[y3,-29.95]]
int value=1380180060345998.244
I=0
base_replace_n
Runge-kutta
  snip 
I=8
base_replace_n
Runge-kutta
[[ 7.05 -2.05 9.95 6.45 -29.95 -0.95 0.95 -29.95 -29.95 ],[0.0003236947435145229824,0.01365914116495716502,0.0003193051763246178940,0.009267060184740950526,0.00002355671589603190298,-0.00001046867630351169487]]

This means that the target function takes the minimal 
0.0003236947435145229824
at the point
 [[x11,7.05],[x12,-2.05],[x13,9.95],[x22,6.45],[x23,-29.95],[x33,-0.95],
  [y1,0.95],[y2,-29.95],[y3,-29.95]]

-29.95 is on the border and the point is not in the interior.
If we enlarge the search domain beyond -30, the target value becomes
too small and we have a numerical instability.
-----------------------------------------------