Advanced Search

mmcdara

6740 Reputation

18 Badges

8 years, 191 days
Contact mmcdara

MaplePrimes Activity


These are questions asked by mmcdara

A strange result which requires explanat...

mmcdara 6740 Maple 2015
denom floating-point numer
December 06 2024
1 3

Can anyone explain me the reason of the last result?
Thanks in advance

> 

restart
> 

kernelopts(version)

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

 (1)
> 

a/n^b;
den := denom(%);
print(cat(`_`$50));

3/n^2;
den := denom(%);
print(cat(`_`$50));

1.23/n^1.65;
den := denom(%);
num := numer(%%);

a/n^b

  

n^b

  

__________________________________________________

  

3/n^2

  

n^2

  

__________________________________________________

  

1.23/n^1.65

  

1

  

1.23/n^1.65

 (2)
 

 

Download What-does-happen-here.mw

Numerical integration using Monte-Carlo ...

mmcdara 6740 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 6740 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 6740 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 6740 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

 

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