It is an equation in x^2, so I would look at 

eq:=x^4 -(3*m+2)*x^2 + 3*m+1; subs(x=sqrt(X), %); factor(%);

          (X-1)*(X-3*m-1)

and its real, positive (--> x real) solutions.

It is an equation in x^2, so I would look at 

eq:=x^4 -(3*m+2)*x^2 + 3*m+1; subs(x=sqrt(X), %); factor(%);

          (X-1)*(X-3*m-1)

and its real, positive (--> x real) solutions.

It would be better to provide your entire sheet

upload it using the 'green error' at the end of the menu bar of the editor of this forum which you see when you post or reply

Essentially: exp(-709) ~ 0 and the procedure tries to find the left and right most z.

First it starts in z = -1 and searches the negative axis to solve log(integrand) = -709
(and takes -1 as solution in case it does not suceed).

Then it searches from -1 to +infinity, avoiding a previous solution (else takes +1).

This is a bit of vague hope, just based on what you 'reported' and what I have seen
by starring at examples.

The rest is 'just coding', nothing special.

That solving to get the range is expensive (and I have no good idea to make it cheaper
like some approximative functions), but something like that is needed:

Your task is very peaked in some very small interval. And I guess NAG has to give up a
search (as it can not know where to continue - using digits=5 tells it almost for sure
that the fct is zero). While the 'sinc' method stretches the interval (IIRC) and Maple
is not fooled by exp(-709) = 0 like hardware methods.


May be you find a scaling where the very interval is of magnitude 1 and not of 1e-3,
and a 'guess', where the center of the peak(s) might be ?

And there exist other integrators as well, if that would not help.

Essentially: exp(-709) ~ 0 and the procedure tries to find the left and right most z.

First it starts in z = -1 and searches the negative axis to solve log(integrand) = -709
(and takes -1 as solution in case it does not suceed).

Then it searches from -1 to +infinity, avoiding a previous solution (else takes +1).

This is a bit of vague hope, just based on what you 'reported' and what I have seen
by starring at examples.

The rest is 'just coding', nothing special.

That solving to get the range is expensive (and I have no good idea to make it cheaper
like some approximative functions), but something like that is needed:

Your task is very peaked in some very small interval. And I guess NAG has to give up a
search (as it can not know where to continue - using digits=5 tells it almost for sure
that the fct is zero). While the 'sinc' method stretches the interval (IIRC) and Maple
is not fooled by exp(-709) = 0 like hardware methods.


May be you find a scaling where the very interval is of magnitude 1 and not of 1e-3,
and a 'guess', where the center of the peak(s) might be ?

And there exist other integrators as well, if that would not help.

And here is the next Java issue: http://seclists.org/fulldisclosure/2013/Jul/172

Thx, that sounds more interesting - so which of the variables/parameters do vary? In in what expected ranges? That would narrow the task. And how much do you expect in terms of 'exactness' - like the 5 digits you used or so?

Thx, that sounds more interesting - so which of the variables/parameters do vary? In in what expected ranges? That would narrow the task. And how much do you expect in terms of 'exactness' - like the 5 digits you used or so?

It is not, that it would be difficult to install Java for the browser(s) from Oracle.

The major point is: for reasons of security and privacy I will not do (in my standard,
which is FF on Win).

As it is (almost) neccessary to allow Javascript for many pages. But in FF that means
to allow Java as well (already selective choices for sited need s th like NoScript as
plugin). Which is a bad idea:

https://panopticlick.eff.org/ or similar

Without Javascript i am unique under ~ 3*10^3 user, with Javascript that would need at
least ~ 3*10^3 * 10^3 users. And having Java it would allow to even access my internal
IP if working within a company account. No.

I only allow that for 1 very special browser installation. And else ignore any site,
which provides solutions needing Java.

I have no really good idea.

It is not clear why you need thousands of those integrals nor why not just let it run
some hours or days - besides the fact of ideas of checking some results :-( 

Where - BTW - I would guess that double exponential would be better ... so what.

Genericall one problem is, that your 'mean's Zo=t, Xo and Yo may be arbitrary (while
they are not in your example). Thus finding peaks systematically is weird (not clear,
whether there is moer than one).

And all that becomes even more complicated, as you introduce additional variables to
modify the initial 'formal' equation.

Having all that (and having lost patients) I would try (ignoring run times in a first
turn, that would not help if results are wrong):

Write the integrand as exp( logIntegrand ), where

logIntegrand:=
  '-1/2*(Zo-2*z)^2 / Sigz^2
   -1/2*(Xpo*z+Xo)^2 / (Sigx^2+Sigxp^2*z^2)
   -1/2*(Yo+Ypo*z)^2 / (Sigy^2+Sigyp^2*z^2)'             # exp term 
  - '1/2*ln((Sigx^2+Sigxp^2*z^2)*(Sigy^2+Sigyp^2*z^2)*Sigz^2)'; # sqrt

In your provided example(s) you will find it to be roughly parabolic (and there would
be a way to estimate ranges in brute, but reasonable ways - in general it is not clear).

Now 'solve' the task, that abs(logIntegrand) <  709 for z, given your data.

Note that this might be empty (i.e. there are constellations, which are valid, but for
all values of z they are beyond) and for me it is not clear, that this is a connected
interval (different 'means' ?).

But for that very range: this are your ranges of integration, where NAG should work (and
for lower precision use something appropriate instead of exp(+-709) = extremal).

If you do not do it in advance then either the integrators will try that and may fail
you are sure they wil not fail (pray for correct results).


To illustrade the problem in your given *simple* case: take t=1e-4 (almost your tasks)
and your test data. Then the integrand is 

  128.036879932896*exp(-32786.8852459016*z^2)*exp(-1/2*(.83e-1*z+.1e-3)^2/(.2272e-9+
  .271750000000000e-6*z^2))/((.2272e-9+.271750000000000e-6*z^2)*(.61968e-14+
  .749333333333333e-7*z^2))^(1/2) 

Now use plot(%, z= -0.002 ..0) to show it (roughly almost 0, peaked in -0.001 having
a value of ~ 2*1e+13 (!)).

If you feed such for an integrator for +-infinity you need lots of trust, that it will
find that very location (searching over 'all' reals).

The method sinc or double exp will scale to a different intervall, but the problem is
still the 'same'.

I have no really good idea.

It is not clear why you need thousands of those integrals nor why not just let it run
some hours or days - besides the fact of ideas of checking some results :-( 

Where - BTW - I would guess that double exponential would be better ... so what.

Genericall one problem is, that your 'mean's Zo=t, Xo and Yo may be arbitrary (while
they are not in your example). Thus finding peaks systematically is weird (not clear,
whether there is moer than one).

And all that becomes even more complicated, as you introduce additional variables to
modify the initial 'formal' equation.

Having all that (and having lost patients) I would try (ignoring run times in a first
turn, that would not help if results are wrong):

Write the integrand as exp( logIntegrand ), where

logIntegrand:=
  '-1/2*(Zo-2*z)^2 / Sigz^2
   -1/2*(Xpo*z+Xo)^2 / (Sigx^2+Sigxp^2*z^2)
   -1/2*(Yo+Ypo*z)^2 / (Sigy^2+Sigyp^2*z^2)'             # exp term 
  - '1/2*ln((Sigx^2+Sigxp^2*z^2)*(Sigy^2+Sigyp^2*z^2)*Sigz^2)'; # sqrt

In your provided example(s) you will find it to be roughly parabolic (and there would
be a way to estimate ranges in brute, but reasonable ways - in general it is not clear).

Now 'solve' the task, that abs(logIntegrand) <  709 for z, given your data.

Note that this might be empty (i.e. there are constellations, which are valid, but for
all values of z they are beyond) and for me it is not clear, that this is a connected
interval (different 'means' ?).

But for that very range: this are your ranges of integration, where NAG should work (and
for lower precision use something appropriate instead of exp(+-709) = extremal).

If you do not do it in advance then either the integrators will try that and may fail
you are sure they wil not fail (pray for correct results).


To illustrade the problem in your given *simple* case: take t=1e-4 (almost your tasks)
and your test data. Then the integrand is 

  128.036879932896*exp(-32786.8852459016*z^2)*exp(-1/2*(.83e-1*z+.1e-3)^2/(.2272e-9+
  .271750000000000e-6*z^2))/((.2272e-9+.271750000000000e-6*z^2)*(.61968e-14+
  .749333333333333e-7*z^2))^(1/2) 

Now use plot(%, z= -0.002 ..0) to show it (roughly almost 0, peaked in -0.001 having
a value of ~ 2*1e+13 (!)).

If you feed such for an integrator for +-infinity you need lots of trust, that it will
find that very location (searching over 'all' reals).

The method sinc or double exp will scale to a different intervall, but the problem is
still the 'same'.

For peaked integrands it is not unlikely the routines may fail.

One should try to 'search' for the range or estimate it (and
sometimes scale.

Will ponder about it.

For peaked integrands it is not unlikely the routines may fail.

One should try to 'search' for the range or estimate it (and
sometimes scale.

Will ponder about it.

Just a minor question: does it really need a Java plugin for the browser?

At least it says so, if calling the pages ...

Up to the convention the result is 0.895019433040495

Find a sheet attached.

MP_spherical.mws

MP_spherical.pdf 


Edited for additional comment

It seems that this is related to you former question "Help with 2D Numerical Integration",
http://www.mapleprimes.com/questions/147530-Help-With-2D-Numerical--Integration

I would not say, that all the symmetries do carry over for that and it is not clear, what
your inputs for the spericals may be in general (though converting to hypergeometrics
and the simplifying should do a lot).

However I think that the main numerical problem (besides fast evaluation) is you first
term in your expression "Vsi", 1/sqrt((1+epsilon-sin(theta)*sin(phi+1/6*Pi))), and you
divide by almost zero for small epsilon (at the boundaries (?))

Have not checked the other 2 terms.