Advanced Search

Dkunb

65 Reputation

4 Badges

2 years, 253 days
Contact Dkunb

MaplePrimes Activity


These are replies submitted by Dkunb

more questions...

Dkunb 65
April 10 2023
0 0

@C_R  How to know the assumption 1/5<x<1 does not fit to the 3rd and 4th root which seem to be complex other than plotting??

Thanks  a lot !

Thank you!...

Dkunb 65
April 09 2023
0 0

@mmcdara Thank you!

Stupid me! Thanks!!...

Dkunb 65
April 07 2023
0 0

@dharr 

Stupid me! Thanks!!

Thanks!!...

Dkunb 65
March 17 2023
0 0

@Carl Love 

Thanks a lot...

Dkunb 65
March 07 2023
0 0

@Hullzie16 

Thank you!...

Dkunb 65
December 30 2022
0 0

@Axel Vogt 

Overflow issue...

Dkunb 65
December 29 2022
0 0

@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

Thank you!...

Dkunb 65
December 22 2022
0 0

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

@acer  Thank you!...

Dkunb 65
November 23 2022
0 0

@acer  Thank you!

@Carl Love  Thank you!...

Dkunb 65
November 23 2022
0 0

@Carl Love  Thank you!

Yes, it is ok...

Dkunb 65
September 19 2022
0 0

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

rho is a function of (p,q)...

Dkunb 65
September 18 2022
0 0

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

Without cutoff...

Dkunb 65
September 14 2022
0 0

@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

 

Thank you again!...

Dkunb 65
September 13 2022
0 0

@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

Thank you again!!...

Dkunb 65
September 13 2022
0 0

@mmcdara Thank you again!!

1 2 3 4 Page 2 of 4