I think that this does not make sense: you want to feed it to a C++, which is
to work in double precision, but your formulae are already too lengthy to feed
it (if they ever can be achieved).

For example: det_2 should not be computed by inverting J_2 and taking its det,
but taking det(J_2) and inverting. Which is the same, but much less expensive.

Already det(J_2) would fill ~ 5.000 pages and it is almost for sure, that this
neither can be feed to C++. Nor (after splitting) it will result in any value
which is reliable.

If you really want to go that way: to not use simplify/size (it will not give
you any advantage here, I tested it).

Perhaps some small advantage by codegen[makeproc] for the matrix and then
codegen[optimize] on that.

Now set up a proc, which computes the det of a 3x3 matrix. And feed it, having
numerical values. May be, some C++ libs already have that.

In this case you want to check the input(s) for fsolve.

For your Vector I get MM[1] = 660049668.737575*a[1](.1e-3) + ...
and you want to solve for a[1](.1e-3), ..., a[18](.1e-3)
What should that mean? You need variables.The last line has false curly brackets, the syntax is {seq(expr,k=1..18)}

It is easy to see, that the function has period = Pi and furthermore has a
point symmetry in Pi/2. So you can restrict to the interval 0 ... Pi/2 and
from that it is clear that there are at least 4 global extrema.

I have no explicit solution, but my guessing is: it is the first local extremum.

Here is brute (!) numerical approach.

f:=(cos(x)-cos((2*n+1)*x))/2/sin(x);
g:=eval(%, x= Pi/2*t); # so we search in t = 0 .. 1

nTst:=2000;
G:=eval(g, n=nTst);
2=numer(G);            # it is 2 at most 
fsolve(%, t=0, complex);
t0:=abs(Re(%));

diff(G,t)=0;
fsolve(%,t=t0);

#plot(G, t=0..nTst*t0);
plot(G, t=0..4*t0);

Edited 01.09.2013 to append a worksheet
http://www.mapleprimes.com/view.aspx?sf=151187/468374/MP_maximum_rational_.mws

Maple is a CAS - I do not understand why to expect a publication feature in top. It just makes it more complicated to be handle: some star instead a bar ... anyway: how about ` *` (holding a blank)?

Your sheet "crashes" for me: you only provide output ... but starring at it
I would aks: "solving for what variables?" Even if they would be linear the
result is parametric. And thus may be complicated.

Write the results into a variable (and only display parts which are needed) - nobody would read 1 Mio characters

?Int

Int( exp(x), x); value(%);

Find appended some coded suggestion. 

Using double precision I get ~ 40 msec for each of those integrals using NAG routines,
working on an usual office machine.

Not unusual - but one of the main problems is to get a reasonable (?!) result at all.


Decreasing required precision (either through epsilon = ... or digits = ...) in the
NAG integration routines does not change run times significantly. However that seems
to lead to errors only.

Shrinking the integration ranges to single precision also gives no advantages.

Expanding the range does seems not to avoid errors which may occure.

I played a bit in trying to make the code a bit more efficient (avoiding the symbolics
I used within it), but that would give only about 10 % in speed or so (but makes it
hard to change).

So only brushing it up for Compiler:-Compile and using that might have advantages.


Looking at some examples (in the enclosed sheet) it seems, that the problem is quite
fragile for inputs. Thus it may be, that a proper scaling or reformulation should be
done (for which one usually needs to know a 'bit' of the application theory and the
possible ranges of all the parameters and variables).

NumIntExample_work.mws
NumIntExample_work.pdf 

If you plot one of the integrals then you see, that it very peaked near zero (try some
interval +- 1e-3 and thus one can guess, that there is no need to integrate over the
full axis.

Indeed the integrand is of form exp(rational of degree (4,2))/sqrt(polynom of degree 4),
where the rational is ~ 2^15*z^2 + ...

Plotting that A:=rational around zero you see, that it like a quadric, but quickly exceeds
values, for which exp makes sense. With 15 Digits that is ~ +- 700.

Now fsolve(A = - 700, z) gives ~ 5*1e-3, beyond exp(A) will be zero and in theory will
dominate the denominator (one can confirm by checking your case of parameters, I think)

I have not look deeper, but that should be the way to look at the numerical problem
and then there should be no problem to use NAG in full precision, there is only a very
small range near zero where the task lives.

(and I have not look for the problem of increasing time)

Though I would not use isolate here:

(r^2*cos(theta)^2+r^2*sin(theta)^2)^2 = a^2*(r^2*cos(theta)^2-r^2*sin(theta)^2);
simplify(%);

                     4    2  2              2
                    r  = a  r  (2 cos(theta)  - 1)

isolate(%, r^2); 
                                     4
                      2             r
                     r  = ----------------------
                           2              2
                          a  (2 cos(theta)  - 1)

combine(%);
                                     4
                          2         r
                         r  = ---------------
                               2
                              a  cos(2 theta)

I was pondering (and that formulation of 'brute integration' certainly is the wrong way).

By symmetry one can reduce to the square over 0 .. Pi (and takes 4 times of that). I get

Int(2/(-x^2+1)^(1/2)*EllipticK(2*I*x*2^(1/2))-
  6*EllipticE(2*I*x*2^(1/2))/(-x^2+1)^(1/2)/(8*x^2+1),x = 0 .. 1);
         1
        /                     1/2                       1/2
       |   2 EllipticK(2 I x 2   )   6 EllipticE(2 I x 2   )
       |   ----------------------- - ----------------------- dx
       |           2     1/2            2     1/2     2
      /         (-x  + 1)            (-x  + 1)    (8 x  + 1)
        0

and after switching to hypergeomtrics and "lower b" for the 2nd term one has to show

1 = Int(1/(-x^2+1)^(1/2)*hypergeom([1/2, 1/2],[1],-8*x^2)-
        3/(-x^2+1)^(1/2)*hypergeom([1/2, 3/2],[1],-8*x^2),x = 0 .. 1);

         1
        /                                 2
       |   hypergeom([1/2, 1/2], [1], -8 x )
  1 =  |   ---------------------------------
       |                2     1/2
      /              (-x  + 1)
        0

                                            2
           3 hypergeom([1/2, 3/2], [1], -8 x )
         - ----------------------------------- dx
                         2     1/2
                      (-x  + 1)

Maple solves that integral - but in useless terms of MeijerG (only numerical).
BTW: at least the first term can be 'solved' using Prudnikov Vol 3 p. 153, 2.4.16.1

I have not done such for longer time. And all that stuff was / is (?) not documented
or explained nicely, it needs some trail-and-error.

Anyway: already in actual programs one would do some checks on data input before
processing - in your case you would check for numerical values.

Then one needs something like MapleToFloat - your first example is just an integer.

And there are commands like executing procedure and not a statement. May be a way
is to put all the critical things into some proc and just call that.


Overall it is somewhat painful work and a start to learn it are the examples which
come with Maple (have not looked, whether they are still delivered).


NB: I see that you have 64 bit ... just be aware that all fits together.

If you want that eq and want to edit it: why do you not look for some Math editor / add in for MS Word ?

I changed your inital data to 'exact' numbers (and l1=r, l2=s to keep it readable
and put it into a reasonable sheet, not the strange standard *.wm).

After looking closer one can reduce the problem in 2 steps: first solve for x and y
with Eq 3 and Eq4.

Use that to feed Eq5 to express s=l1 by r=l2 and lambda.

Now you are at Eq1 and and Eq2, but those are 'only' in r and lambda.

These are 'algebraic' equations in r and lambda, r enters linear and 'solve' can
provide a solution - there are 6 of them (some are conjugated I think).

Find a sketchy sheet appended, hope you can finish your task from that.

MP_nonlinear_4variab.mws

solve for x or r?

degree in x it is = 16 and in r it is = 15

PS: whatever ... I doubt it can be done. Assuming it is x then divide out
((r*x^2+x-1)*(r^2*x^2-r*x-r+1) ) giving degree = 12. Now setting r = 1
Maple tells it is irreducible.
Otherwise set x=1 giving degree = 12 in r and seeting x=1 it turns out
to be a square of degree 6 and that is irreducible as well