Advanced Search

mmcdara

6730 Reputation

18 Badges

8 years, 186 days
Contact mmcdara

MaplePrimes Activity


These are questions asked by mmcdara

Numerical integration using Monte-Carlo ...

mmcdara 6730 Maple
integration numeric monte-carlo
November 23 2024
2 5

I have a lot of questions about using evalf/Int with Monte-Carlo like methods.
They are red written in the attached file and concern several different points.
I would like you to answer each of them and not to focus on a specific one.
Thanks in advance.

method=_MonteCarlo uses the Fortran procedure d01gbc (maybe Maple rewritten?): the original NAG procedure is described here d01gbc

> 

restart
> 

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

 (1)
> 

domain := 0.8..3.;
f := x -> 1/(1 + sinh(2*x)*log(x)^2);

.8 .. 3.

  

proc (x) options operator, arrow; 1/(1+sinh(2*x)*log(x)^2) end proc

 (2)
> 

infolevel[`evalf/int`] := 4:

# I understand _CubaSuave used a random sample of size 134000 of the (x, y) integration domain.
# Am I right?
 
evalf(Int(f(x), [x=domain, y=0..1], method=_CubaSuave, epsilon=1e-4));
printf("\n%s\n", cat("-"$100));

# Did _CubaDivonne used a random sample of size 3256 ?
 
evalf(Int(f(x), [x=domain, y=0..1], 'method=_CubaDivonne', epsilon=1e-4));
printf("\n%s\n", cat("-"$100));

# What is the true number of points _MonteCarlo used?
 
evalf(Int(f(x), [x=domain, y=0..1], method=_MonteCarlo, epsilon=1e-4));

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: transformed original integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

  

2.2/(1+sinh(4.4*x+1.6)*ln(2.2*x+.8)^2)

  

cuba: integration completed successfully
cuba: # of integrand evaluations: 134000
cuba: estimated (absolute) error: 6.73167e-05
cuba: chi-square probability that the error is not reliable: 1
cuba: number of regions that the domain was divided into: 134

  

HFloat(0.6770776956970137)

  


----------------------------------------------------------------------------------------------------
Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: transformed original integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

  

2.2/(1+sinh(4.4*x+1.6)*ln(2.2*x+.8)^2)

  

cuba: integration completed successfully
cuba: # of integrand evaluations: 3256
cuba: estimated (absolute) error: 6.37672e-05
cuba: chi-square probability that the error is not reliable: 0
cuba: number of regions that the domain was divided into: 22

  

HFloat(0.6768372810943345)

  


----------------------------------------------------------------------------------------------------
Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

trying d01gbc (nag_multi_quad_monte_carlo)
d01gbc: epsrel=.1e-3; minpts=0; maxpts=500000000; method=2; cont=0
d01gbc: procedure for evaluation is:
proc (X) 1/(1.+sinh(2.*X[1])*ln(X[1])^2) end proc

d01gbc: result=.676840185506968228
d01gbc: relerr=.103631215397982221e-4; usedpts=1044
result=.676840185506968228

  

.6768401855

 (3)
> 

# By the way, why is infolevel[`evalf/int`] output discarded?

evalf(Int(f(x), [x=domain, y=0..1], method=_CubaSuave, epsilon=1e-4));
printf("\n%s\n", cat("-"$100));
evalf(Int(f(x), [x=domain, y=0..1], method=_MonteCarlo, epsilon=1e-4));

HFloat(0.6770776956970137)

  


----------------------------------------------------------------------------------------------------

  

.6768401855

 (4)
> 

# Let us suppose I don(t care of the accuracy of the result as I am capable to assess it
# by some other way:
#      How can I use a Monte-Carlo integration method with a given number of points?
#      Why there is no estimation of the integral returned?

infolevel[`evalf/int`] := 0:
infolevel[`evalf/int`] := 4:
evalf(
  Int(
    f(x), [x=domain, y=0..1]
    , method=_CubaVegas
    , methodoptions=[minimalpoints = 10^3, maximalpoints = 10^3]
  )
);

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: transformed original integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

  

2.2/(1+sinh(4.4*x+1.6)*ln(2.2*x+.8)^2)

  

cuba: could not attain requested accuracy
cuba: # of integrand evaluations: 1000
cuba: estimated (absolute) error: 0.025423
cuba: chi-square probability that the error is not reliable: -999
evalf/int: error from Control_multi was:
"could not attain requested accuracy; try increasing epsilon or absepsilon or maximalpoints"

  

Int(Int(1/(1.+sinh(2.*x)*ln(x)^2), x = .8 .. 3.), y = 0. .. 1.)

 (5)
> 

infolevel[`evalf/int`] := 0:
infolevel[`evalf/int`] := 4:
evalf(
  Int(
    f(x), [x=domain, y=0..1]
    , method=_CubaSuave
    , methodoptions=[minimalpoints = 10^3, maximalpoints = 10^3]
  )
);

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: transformed original integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

  

2.2/(1+sinh(4.4*x+1.6)*ln(2.2*x+.8)^2)

  

cuba: could not attain requested accuracy
cuba: # of integrand evaluations: 1000
cuba: estimated (absolute) error: 0.025423
cuba: chi-square probability that the error is not reliable: -999
cuba: number of regions that the domain was divided into: 1
evalf/int: error from Control_multi was:
"could not attain requested accuracy; try increasing epsilon or absepsilon or maximalpoints"

  

Int(Int(1/(1.+sinh(2.*x)*ln(x)^2), x = .8 .. 3.), y = 0. .. 1.)

 (6)
> 

infolevel[`evalf/int`] := 0:
infolevel[`evalf/int`] := 4:
evalf(
  Int(
    f(x), [x=domain, y=0..1]
    , method=_MonteCarlo
    , methodoptions=[minimalpoints = 10^3, maximalpoints = 10^3]
  )
);

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

  

1/(1+sinh(2*x)*ln(x)^2)

  

evalf/int: error from Control_multi was:
"NAG d01gbc expects epsilon >= 0.5e-4, but received %1", .5000000000e-9

  

Int(Int(1/(1.+sinh(2.*x)*ln(x)^2), x = .8 .. 3.), y = 0. .. 1.)

 (7)
> 

# Example 8 of the reference given in the question text

evalf(Int(4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2, [seq(x||i=0..1, i=1..4)], method=_MonteCarlo, epsilon=1e-2))

.5753724460

 (8)
> 

# Why is not the seeed updated ?

for j from 1 to 3 do
  Threads:-Sleep(10):
  evalf(Int(4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2, [seq(x||i=0..1, i=1..4)], method=_MonteCarlo, epsilon=1e-2))
end do;

.5753724460

  

.5753724460

  

.5753724460

 (9)
> 

# Here a new seed is forced at each iteration but the estimations of the integral are always the same.
# This is impossible as these estimationsare the realizations of a random variable: so they cannot be identical.
#
# So, does evalf/Int always use the same internal seed?
# More importantly: does it really use a random generator?
for j from 1 to 3 do
  seed := randomize(rand()());
  J := evalf(Int(4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2, [seq(x||i=0..1, i=1..4)], method=_MonteCarlo, epsilon=1e-2)):
  print('seed' = seed,  'J' = J);
end do:

Control_multi: integrating on [0, 0, 0, 0] .. [1, 1, 1, 1] the integrand

  

4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2

  

trying d01gbc (nag_multi_quad_monte_carlo)
d01gbc: epsrel=.1e-1; minpts=0; maxpts=500000000; method=2; cont=0
d01gbc: procedure for evaluation is:
proc (X) 4.*X[1]*X[3]^2*exp(2.*X[1]*X[3])/(1.+X[2]+X[4])^2 end proc
d01gbc: result=.574090186993690188
d01gbc: relerr=.675204841605886539e-2; usedpts=8064
result=.574090186993690188

  

seed = 233366062458, J = .5740901870

  

seed = 887991988815, J = .5740901870

  

seed = 416683078956, J = .5740901870

 (10)
> 

# This is what one expects from Monte-Carlo integration

use Statistics in
  for i from 1 to 10 do
    seed := randomize(rand()());
    S := Sample(Uniform(op(domain)), 100):
    print('seed' = seed,  'J' = Mean(f~(S)) * (- `-`(op(domain))) );
  end do:
end use:

seed = 36995932795, J = HFloat(0.6280945376385327)

  

seed = 976943479321, J = HFloat(0.7320842356827296)

  

seed = 542869221880, J = HFloat(0.6519958887669823)

  

seed = 303692769322, J = HFloat(0.6134976803250747)

  

seed = 443233702046, J = HFloat(0.6304929479488379)

  

seed = 881136125112, J = HFloat(0.6528328283998008)

  

seed = 708639694457, J = HFloat(0.6601745158455014)

  

seed = 675811574035, J = HFloat(0.7640187396123023)

  

seed = 164464632871, J = HFloat(0.7949627835070758)

  

seed = 604897894346, J = HFloat(0.5922357495886873)

 (11)
 

 

Download Numerical_integration_Using_MC_methods.mw

A few problems with limit: how to fix th...

mmcdara 6730 Maple 2015
integration eval limit
October 02 2024
1 1

It seems that Maple needs more help than necessary:

> 

restart:
> 

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

 (1)
> 

expr:= A+B*limit(f(x), x=+infinity);
eval(expr, limit(f(x), x=+infinity)=1)

A+B*(limit(f(x), x = infinity))

  

A+B

 (2)
> 

expr:= A+B*limit(2*f(x), x=+infinity);

eval(expr, limit(f(x), x=+infinity)=1);     # Shouldn't this return A+2*B
eval(expr, limit(2*f(x), x=+infinity)=2);   # Can I avoid doing this?

A+B*(limit(2*f(x), x = infinity))

  

A+B*(limit(2*f(x), x = infinity))

  

A+2*B

 (3)
> 

expr:= A+B*limit(f(x)^2, x=+infinity);

eval(expr, limit(f(x), x=+infinity)=1);     # Shouldn't this return A+B
eval(expr, limit(f(x)^2, x=+infinity)=1);   # Can I avoid doing this?

A+B*(limit(f(x)^2, x = infinity))

  

A+B*(limit(f(x)^2, x = infinity))

  

A+B

 (4)
> 

expr:= A+B*limit(2*f(x)^2, x=+infinity);

eval(expr, limit(f(x)^2, x=+infinity)=1);    # Shouldn't this return A+2*B
eval(expr, limit(2*f(x)^2, x=+infinity)=2);  # Can I avoid doing this?

A+B*(limit(2*f(x)^2, x = infinity))

  

A+B*(limit(2*f(x)^2, x = infinity))

  

A+2*B

 (5)
 

 

Download limits.mw

Why don't the commands labelled "Shouldn't this return.." do the job?

TIA

Why do I get this unevaluated expression...

mmcdara 6730 Maple 2015
integral substitution indefinite-integration
September 21 2024
1 12

I have an expression equal to the sum of N terms of the form Int(fn=1..N(x), x) and I want to replace each fn(x) by its Taylor (or series) expansion.

When the integrals are definite, like J1 below, I can easily obtain a new expression (K1) where the integrand has been replaced by some expansion.
But when the integral is indefinite, like J2, I get an evaluated expression for K2.

It seems I have to do some gymnastic (J3 --> K3) to get what I want

> 

restart
> 

J1 := Int(sin(p*x), x=0..1);
K1 := eval(J1, Int = ((a, b) -> Int(mtaylor(a, x=0, 5), b)));

Int(sin(p*x), x = 0 .. 1)

  

Int(p*x-(1/6)*p^3*x^3, x = 0 .. 1)

 (1)
> 

# undefined integration

J2 := Int(sin(p*x), x);

`Expected result` = Int(p*x-(1/6)*p^3*x^3, x);

K2 := eval(J2, Int = ((a, b) -> Int(mtaylor(a, x=0, 5), b)));

Int(sin(p*x), x)

  

`Expected result` = Int(p*x-(1/6)*p^3*x^3, x)

  

eval(Int(sin(p*x), x), {Int = (proc (a, b) options operator, arrow; Int(mtaylor(a, x = 0, 5), b) end proc)})

 (2)
> 

# undefined integration using Intat

J3 := Intat(op(1, J2), op(2, J2)=y);
eval(%, Intat = ((a, b) -> Intat(mtaylor(a, x=0, 5), b))):

K3 := IntegrationTools:-Change(convert(%, Int), y=x, x)

Intat(sin(p*x), x = y)

  

Int(p*x-(1/6)*p^3*x^3, x)

 (3)
 

 

Download Integration.mw

Why do I get this unevaluatedform for K2?
Do I have to use Intat to get K3?

Thanks in advance

Why is orthopoly not seen as a package w...

mmcdara 6730 Maple 2015
package procedure orthopoly
August 26 2024
1 6


For years I observe that package orthopoly is not considered as a package within a procedure.
Finally I have decided to ask for clarifications: Can someone explain me why procedure f generates an error?
 

> 

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

 (1)


Without any procedure

> 

restart
> 

H(2, x)

H(2, x)

 (2)
> 

with(orthopoly):
> 

H(2, x)

4*x^2-2

 (3)


Within a procedure

> 

restart
> 

type(orthopoly, package)

true

 (4)
> 

f := proc(m)
  uses orthopoly:
  H(m, x)
end proc:

Error, `orthopoly` is not a module or member

  
> 

g := proc(m)
  orthopoly:-H(m, x)
end proc:

g(2)

4*x^2-2

 (5)
> 

 


Download meaning.mw

TIA

 

How to upload a large zip file?...

mmcdara 6730 MaplePrimes
zip attachment
August 20 2024
1 3


A few times ago I asked a question about the best way to upload a bunch of files and was told that I had to create a zip file and upload it.

The post I want to create is about the packing of "balls" in "unit boxes" (in arbitrary dimension). I have about 10 Maple worksheets and more than one thousand result files (likely about a hundred of hours of computation).
The size of the zip file is about 260 Mb and the link cannot be inserted (it does not even appear in the list of the uploadable files).

What do you recommend?

TIA

1 2 3 4 5 6 7 Last Page 1 of 44