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.082,
S2 = -0.959,
S3 = 0.131,
and
A=domain([-30,-30,-30,-30, -30, -30.0, -30.0, -30, -30],
[ 30, 30, 30, 30, 30, -0.01, 30, -0.001, 32] )
=[-30,10] x [-30,10] x ... x [-30,-0.001] x [-30,32]
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");
test9_64_demo();
The starting point of the holonomic gradient descent method is taken at
Pos = [6.695,0.048,3.502,
6.459,22.902,
-13.693,
1.712,-28.99,29.992]
This starting point is found by wood's method and apply the holonomic gradient
descent for the output (test8_64()).
During the first execution, the big data file "tmp-sei13c-r1.ab" is
automatically downloaded by "wget". It will take a while.
The file "tmp-sei13c-r1.ab" is the Pfaffian system for the Fisher-Bigham
integral on S^2.
The output of this function will be
[873,2,[7.065,-0.032,3.422,5.339,24.922,-13.693,1.642,-31.99,31.992],[0.4373096253840751950,17.29022464267128727,1.826327122835943145,13.41174509764664094,0.05896860012506528815,-0.2036217013092417161]]
point = [[x11,7.065],[x12,-0.032],[x13,3.422],[x22,5.339],[x23,24.922],[x33,-13.693],[y1,1.642],[y2,-31.99],[y3,31.992]]
int value=4126880502191184.670
This means that the target function takes the minimal
0.4373096253840751950
at the point
point = [[x11,7.065],[x12,-0.032],[x13,3.422],[x22,5.339],[x23,24.922],[x33,-13.693],[y1,1.642],[y2,-31.99],[y3,31.992]]
-31.99 and 31.992 are on the border of the variables y2 and y3.
demo10_64() shows that the target value changes only 10^(-5) when we enlarge
the domain.
-----------------------------------------------