Dkunb

35 Reputation

3 Badges

1 years, 34 days

MaplePrimes Activity


These are replies submitted by Dkunb

@mmcdara I tried to run the code based on your code with different parameter, it gave me overflow. What should I do to fix this issue?

And how can I try to find the correct method and options for numerical integration? ( for example? )

Thanks in advance

Negativity_(v13_beta_gamma_mu).mw

@mmcdara I stopped the process because it seemed to take a long time. Thank you so much!

@acer  Thank you!

@Carl Love  Thank you!

@Carl Love For the purpose of the last expression, your expression is ok. Thank you!

@Carl Love rho is a function of (p,q) and rho is zero on the surface.(boundary of the volume)

@acer What should I change in your code if I do not want to use cutoff? Without cutoff it gives me overflow.
 

restart;

F:=kappa->kappa:
f:=(alpha,delta)->exp(-abs(F(kappa))^2*(1+delta^2)/2-abs(F(kappa))*alpha)/abs(F(kappa)):

L:=(alpha,delta,Lambda) ->
  (lambda^2*exp(-alpha^2/2)/4)*(Int(f(alpha,delta),kappa= -infinity..-Lambda)+Int(f(alpha,delta),kappa= Lambda..infinity)):

CodeTools:-Usage(evalf(L(4,1,0.001)));

memory used=1.88MiB, alloc change=0 bytes, cpu time=157.00ms, real time=404.00ms, gc time=0ns

 

0.8209373770e-3*lambda^2

(1)

X:=Int(exp(-kappa^2*(delta^2 + 1)/2)*cos(kappa*beta)*cos(kappa*gg)/(kappa),
       kappa = Lambda .. infinity):

Y:='Int(-exp(-kappa^2*(delta^2 + 1)/2)*cos(kappa*beta)*Im(exp(kappa*gg*I)*erf( (gg+kappa*I)/sqrt(2)))/(kappa), kappa = Lambda .. infinity)':

theJ:='sqrt(X^2+Y^2)/2*exp(-alpha^2/2)*lambda^2':

J:=unapply(subs(infinity=6, theJ), [alpha,delta,Lambda,beta,gg]): # cut off

#J:=unapply(theJ, [alpha,delta,Lambda,beta,gg]):

thisJ:=unapply(subsindets(J(alpha,1,0.001,beta,3),
                          'specfunc(anything,Int)',
                          'q -> Int(op(q), digits=15, epsilon = 1e-4, method = _d01akc)'),[alpha,beta]):

#thisJ:=unapply(subsindets(J(alpha,1,0.001,beta,3),
                          'specfunc(anything,Int)',
                          'q -> Int(op(q), digits=15, epsilon = 1e-4, method = _d01amc)'),[alpha,beta]):

evalf(thisJ(4,8));

Error, (in evalf/int) numeric exception: overflow

 

N:=proc(beta,alpha) option numeric; local res;
     Digits:=15;

     if not [beta,alpha]::[numeric,numeric] then
     return 'procname'(args); end if;
     res:=thisJ(alpha,beta) - L(alpha,1,0.001);
     res:=evalf[15](res);
     res:=evalf[10](res);
     if type(res/lambda^2,numeric) then
     res/lambda^2; else undefined; end if;
end proc:

forget(evalf):
CodeTools:-Usage( evalf(N(4,8)) );

Error, (in evalf/int) numeric exception: overflow

 

#forget(evalf);
#CodeTools:-Usage( plots:-contourplot(N(beta,alpha), beta=0..10, alpha=0..2, grid=[7,7]) );

#forget(evalf);
#CodeTools:-Usage( plots:-contourplot(N(beta,alpha), beta=0..10, alpha=0..10,
#                                     contours=[-1e-8, 1e-8], grid=[17,17]) );

N:=subsop(4=NULL,eval(N)): # clear N's remember table
forget(evalf);
time[real](assign('PP',[seq([beta,fsolve(alpha->N(beta,alpha),1..10)],
                            beta=0..10,1.0)]))*`seconds real time`;

132.154*`seconds real time`

(2)

#eval(PP,1);

PPSPL:=Interpolation:-Interpolate(PP[..,1],PP[..,2]):

plot(PPSPL, 0..10, view=[0..10,0..10],
     filled=true, color=gray, labels=[beta,alpha], title=(gamma=3),
     transparency=0, background="aa0000");

 

N:=subsop(4=NULL,eval(N)): # clear N's remember table
forget(evalf);
time[real](assign('PNZ',
                  [seq([bval,
                        RootFinding:-NextZero(subs(beta=bval,
                                                   alpha->N(beta,alpha)),
                                              1.0)],bval=0..10,1.0)]) )*`seconds real time`;

123.808*`seconds real time`

(3)

#eval(PNZ,1);

PNZSPL:=Interpolation:-Interpolate(PNZ[..,1],PNZ[..,2]):

plot(PNZSPL, 0..10, view=[0..10,0..10],
     filled=true, color=gray, labels=[beta,alpha], title=(gamma=3),
     transparency=0, background="aa0000");

 

N:=subsop(4=NULL,eval(N)):
forget(evalf);
Grid:-Set(N);Grid:-Set(L);Grid:-Set(thisJ);Grid:-Set(J);Grid:-Set(f);Grid:-Set(F);
time[real](assign('PGR',
                  [Grid:-Seq['tasksize'=max(1,floor(10/kernelopts('numcpus')))]([beta,
                                                                                 fsolve(N(beta,alpha),
                                                                                        alpha=1..10)],
                                                       beta=0..10,1.0)]) )*`seconds real time`;

61.065*`seconds real time`

(4)

#eval(PGR,1);

PGRSPL:=Interpolation:-Interpolate(PGR[..,1],PGR[..,2]):

plot(PGRSPL,0..10, view=[0..10,0..10],
     filled=true, color=gray, labels=[beta,alpha], title=(gamma=3),
     transparency=0, background="aa0000");

 

 


 

Download Negativity_v1_acU0_dk1.mw

 

@acer  Thank you very much for your explanation and advice on posting question!!

I used ArrayInterpolation. I wondered what causes the feature around (6,6). I think I know why beacuase of your explanation.

I tried to understand your codes but it has been difficult for me to understand it becuase I am still new to Maple. 
 

restart;

with(CurveFitting)

[ArrayInterpolation, BSpline, BSplineCurve, Interactive, LeastSquares, Lowess, PolynomialInterpolation, RationalInterpolation, Spline, ThieleInterpolation]

(1)

with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

(2)

alpha := <seq(0..10,evalf(10/50))>:
beta := <seq(0..10,evalf(10/50))>:

excelfile:= FileTools:-JoinPath(["C:","Users","aimer","OneDrive","Desktop","Msc Thesis","Maple ref","N_data.xlsx"]);

"C:\Users\aimer\OneDrive\Desktop\Msc Thesis\Maple ref\N_data.xlsx"

(3)

NN:=ImportMatrix(excelfile,source=Excel):

#?ImportMatrix;

#NN:=ImportMatrix(matlabData, source=MATLAB);

#currentdir();

 

#contourplot(ArrayInterpolation([beta,alpha],NN,[x,y]),x=0..10,y=0..10,contours=[0]);

NN2 := proc(x,y):
  if (type(x,numeric) and type(y,numeric)) then:
     ArrayInterpolation([beta,alpha],NN,[[x],[y]])[1,1];
  else:
     'procname'(args):
  fi:
end proc:

#NN2:=(x,y)->ArrayInterpolation([beta,alpha],NN,[[x],[y]]);

NN2(1,1);

HFloat(-0.14351754)

(4)

 

contourplot(NN2(x,y),x=0..10,y=0..7.5,contours=[0],filledregion=true,coloring=["Gray","Red"], labels=[typeset('beta'),typeset('alpha')],title=[gamma=3],
     transparency=0, background="aa0000");

 

#listcontplot(ArrayInterpolation([beta,alpha],NN,[<seq(0.01*i,0..10)>,<seq(0.01*i,0..10)>]);

 

#?ArrayInterpolation

?contourplot

 


 

Download N_ArInterpol_contplot.mw

@mmcdara Thank you again!!

@mmcdara Thank you so much!

@Axel Vogt Thank you so much!

@acer Thank you with all my heart!!

@Axel Vogt I would like to know what ' ' in Y is doing and " %" in -Usage("%"=) is doing. 

How should I find N (beta, gamma) using your J?

Thank you in advance,

@tomleslie Thank you so much!!

1 2 Page 1 of 2