mmcdara

6710 Reputation

18 Badges

8 years, 153 days

MaplePrimes Activity


These are questions asked by mmcdara

I was developping an answer to @jennierubyjane  but the question disappeared meanwhile.

Here is the answer hopping for someone to restore the initial question (if it still matters?).
To  @jennierubyjane : comments and advices are in Maple 1D font within your original worksheet.

EDITED : Answer has been displaced here  233262-How-Do-I-Solve-Error-in-Dsolvenumericbvpconvertsys

I want to compute the series expansion of i3_r wrt (x, y, z) at point (x=y=z=0):

i2   := (x,y) -> -(1/2)*I*(exp(I*x)*(sin(x)/x)-exp(I*y)*(sin(y)/y))/(x-y):
i3_r := -(1/2)*I*(i2(y,z)-i2(y,x))/(z-x);

My first attempt was to compute this mulltiple series expansion this way:

ordre := 3:
sx := convert( series(i3_r, x, ordre), polynom);
sy := convert( series(sx  , y, ordre), polynom);
sz := convert( series(sy  , z, ordre), polynom);

But this gives me sy=sz=0 whatever the expansion order.

I then do this:

sx :=              convert(series(i3_r , x, ordre), polynom):
sy := add(map(u -> convert(series(u    , y, ordre), polynom), [op(expand(sx))])):
sz := add(map(u -> convert(series(u    , z, ordre), polynom), [op(expand(sy))]));

and obtained non zero results for both sy and sz (but are they are correct ?).

Could you explain me what happens and tell me how to find the series expansion of i3_r wrt (x, y, z) ?

TIA

In the code above expression #1 contains the term D(Phi)(X): how can I replace the "abstract" function Phi by a "specific" one  phi = k*X(t) ?

restart:
alias(X = X(t)):
alias(Phi=Phi(X(t))):
kin := 1/2*M*diff(X, t)^2 + 1/2*m*diff(Phi, t)^2;  #1

# attempt to replace Phi by phi=k*X(t))
alias(phi=k*X(t)):
eval(kin, Phi=phi);                                #2 ... I expected phi would replace Phi

The attached file  How_to_eval.mw contains a more detailed explanation and a more general situation.

Thanks in advance

 

Let  

f := beta/(2*sigma*GAMMA(1/beta))*exp(-(abs(x-mu)/sigma)^beta);

where mu::real, beta > 0, sigma > 0.
How can we obtain the expression of F?

F := int(f, x=-infinity..s)


Reverse problem (a priori simpler):
It's known that 

F := 1/2+signum(x-mu)/(2*GAMMA(1/beta))*(GAMMA(1/beta)-GAMMA(1/beta, abs((x-mu)/sigma)^beta))

How can we check that diff(F, x)=f?
Even with the assumptions on beta, mu and sigma and additional assumptions x>0 or x<0, I can't verify that diff(F, x)=f.

 

I want to estimate numerically the value of P (that is the probability that f exceeds the value 12 when x, y and z are uniformly distributed within the box [-Pi, Pi]3).

restart
f := sin(x)+7*sin(y)^2+0.1*z^4*sin(x);
Omega := [x, y, z] =~ [(-Pi..Pi)$3]:
                                2        4       
               sin(x) + 7 sin(y)  + 0.1 z  sin(x)


h := Heaviside(f-12):
P := Int(h, Omega);


Here is a simple Monte Carlo estimation of P.

f_MC  := x -> sin(x[1])+7*sin(x[2])^2+0.1*x[3]^4*sin(x[1]); 
h_MC  := x -> Heaviside(f_MC(x) - 12);
omega := -Pi, Pi;
                         
P_MC := proc(N)
  local Z := Statistics:-Sample(Uniform(omega), [N, 3]):
  local F := Vector(N, n -> h_MC(Z[n])):
  local K := add(F):
  local P := evalf(add(F)/N):
  local q := Statistics:-Quantile(Normal(0, 1), 0.005, numeric):
  local e := -q*sqrt(P*(1-P)/N):
  printf("A bilateral 99%% confidence interval that Prob(f > 12) is %1.3e +/- %1.2e\n", P, e);
end proc:

P_MC(10^6)
A bilateral 99% confidence interval that Prob(f > 12) is 1.643e-02 +/- 3.27e-04

 

I was hoping to get a value for P using  evalf/Int with some method.
But I didn't succeed with any of the methods for multiple integration:

  • The following command returns 0 for n=1 and 2 and Int(h, Omega) if n >=3
    evalf( Int(h, Omega, 'epsilon=1e-n', 'method=_MonteCarlo') );
    

     

  • And all my attempts with methods from the Cuba library have resulted in an unevaluated Int(h, Omega) integral. 
     

Could you help me to estimate P using  evalf/Int ?
TIA

First 12 13 14 15 16 17 18 Last Page 14 of 44