acer

26592 Reputation

29 Badges

17 years, 0 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@ijuptilk You have now stated that you expect,

  ((2*I)^2)/0.002;

to return,

   -2000.000000 * I

instead of,

   -2000.000000

Why do you expect that?

The original .mw worksheet attachment has been removed from this Question. Here is an earlier copy.

Case3_Ijuptilk_Contourplot_pmin_290822_orig.mw

It is very poor etiquette on the part of whoever removed the original.

@Dkunb You should explain -- properly, clearly, and adequately -- what the coloring and lines represent on those plots.

@Dkunb 

The internal plotting routines will attempt floating-point evaluation, using evalf (or evalhf) as appropriate.

That's how the plot can still be made, even if N is returning only inert Int calls from L.

You are free to call evalf(N(...)) to get any individual/scalar float result, if you prefer. Or you could add an evalf call in the L procedure itself. Of you could call evalf(L(...)) from within N. It's your choice.

@ijuptilk The code inside doCalc is not particulalrly good. It is not fast, for any single x-y input point. So the 5x5 grid is not very fast, and takes a few minutes on my machine.

And now you talk of an 81x81 grid. That is 260 times as many x-y points as a 5x5 grid. Naturally, using the same approach that would be much slower.

You had posted versions of some of it in several older postings, eg. here. I had looked at some of those and made a few suggestions which you appear since to have mostly not taken up. Eg. I suggested parallelization using the Grid package as an alternative to using plot3d to construct the 2-D data.

Parallelization is possible, but it won't give you an enormous speedup. If you raise the grid to 81x18 and make it 250 times slower then a 3-8 times speedup from parallelization won't help much. It's also a back-to-front approach to improving efficiency. You should focus first on the timing of a single call to doCalc, and only after that is good look at parallelization.

You also have a persistent habit of not explaining what exactly it is that you are trying to compute in adequate detail. Without adequate detail it is very difficult to make suggestions that pertain to anything but small aspects of the existing approach. I suspect that large portions of your approach could be replaced, but haven't received enough explanation to do so.

There is also some detail on the root-finding portion, here. And I also deleted my responses from at least one other thread (in which some details appeared), in frustration over the lack of complete detail, the way you kept changing and adding details, your omission in several prior Questions (until after people had put effort into your problem) of the important detail that you are using a very old Maple version, and your reluctance to utilize some suggested improvements.

I likely won't be able to spend more time on improving it.

@ijuptilk I ran the following in Maple 18.02.

I cast the data to a Matrix, before passing to listcontplot.

I used method=QR in the call to LinearSolve below. I don't know enough about what you're computing to pass judgement on whether that a valid approach. What do you think ought to happen in the case that your linear system is, say, underdetermined?

You still haven't bothered to upload an attachment with a worksheet that contains any corrections to the x-axis (range) that you mentioned in a Reply comment on your own Question.

restart:

kernelopts(version);

`Maple 18.02, X86 64 LINUX, Oct 20 2014, Build ID 991181`

Procedure

 

``

doCalc:= proc( xi , u)  #u is the \bar(H): normalize magnetic field magnitude,
                          # where H = bar(H)*H__a
                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.1e-3,
                       eta__1:= 0.240e-1, chi:= 1.219*10^(-6),
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       theta_init:= theta0*sin(Pi*z/d),
                       H__a:= Pi*sqrt(k__1/chi)/d,
                       H := u*H__a,     
                       c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2)),
                       w := chi*H^2*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),
                       n:= 50,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,                        nstar,soln3, Imagroot1, Imagroot2, AreTherePurelyImaginaryRoots;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

  g:= q-(1-alpha)*tan(q)+  c*tan(q) + w*q:
  f:= subs(q = x+I*y, g):
  b1:= evalc(Re(f)) = 0:
  b2:= y-(1-alpha)*tanh(y) -(alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 0;
  qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
  OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

  ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
  qstarTemporary:= min(ModifiedOddAsym);
  indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
  qstar2:= OddAsymptotes(indexOfqstar2);

# Compute complex roots

  AreThereComplexRoots:= type(true, 'truefalse');
  try
   soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
   soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
   qcomplex1:= subs(soln1, x+I*y);
   qcomplex2:= subs(soln2, x+I*y);
   catch:
   AreThereComplexRoots:= type(FAIL, 'truefalse');
  end try;

AreTherePurelyImaginaryRoots:= type(true, 'truefalse');
  try
# Compute the rest of the Roots
  soln3:= simplify~(fsolve({ b2}, { y = 0 .. infinity}),zero);  
  Imagroot1:=subs(soln3, I*y);
  Imagroot2:= -Imagroot1;
  catch:
   AreTherePurelyImaginaryRoots:= type(FAIL, 'truefalse');
  end try;

## odd and all asymptotes
  OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
       
 if (xi > 0) then
  if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2,op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2, Imagroot2, op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op      (Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif not AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [Imagroot2, op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q =               AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif not AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] ..      AllAsymptotes[i+1])), i = 1 .. n)];
  end if:
else
if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2,op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2, Imagroot2, op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op      (Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif not AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q =                          AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif not AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] ..      AllAsymptotes[i+1])), i = 1 .. n)];
  end if:
end if:

# Remove the repeated roots if any & Redefine n

  qq:= MakeUnique(gg):
  m:= numelems(qq):

## Compute the `τ_n`time constants

   for i to m do
   p[i] := simplify(evalf(gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2)),':-zero'):
   end do:

## Compute θ_n from initial conditions

  for i to m do
  Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
  end do:

## Compute integral coefficients for off-diagonal elements θ_n matrix

   printlevel := 2:
   for i to m do
       for j from i+1 to m do
           b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
           b[j, i] := b[i, j]:
               aa[i, j] := b[i, j]:
               aa[j, i] := b[j, i]:
        end do :
   end do:

## Calculate integral coefficients for diagonal elements in theta_n matrix

  for i to m do
  aa[i, i] := int(Efun[i]^2, z = 0 .. d):
  end do:

## Calculate integrals for RHS vector

  for i to m do
  F[i] := int(theta_init*Efun[i], z = 0 .. d):
  end do:

## Create matrix A and RHS vector B

  A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
  B := Vector([seq(F[i], i = 1 .. m)]):

## Calculate inverse of A and solve A*x=B

  #Ainv := 1/A:
  #r := MatrixVectorMultiply(Ainv, B):
  r := LinearSolve(A,B);

## Define Theta(z,t)
  theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):

## Compute v_n for n times constant

  for i to m do
  v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):
  end do:

## Compute v(z,t) from initial conditions
  for i to m do
  Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
  end do:

## Define v(z,t)
   v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):

## compute minimum value of tau
  minp:=min(seq( Re(p[i]), i = 1 .. m) );
  member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):

## Return all the plots
     return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, H, H__a,qq;
     end proc:

Run the calculation for supplied value of 'xi

 

NULL

``

NULL

NULL

NULL

NULL

NULL

NULL

NULL

ans:=[doCalc(0.1, 2)]:
evalf(ans[9]);
evalf(ans[10]);

69.69853431

34.84926715

 

Contour plot for pmin (minimum  of tau) for different xi and H values

 

NULL

NULL

with(plots):
d:= 0.2e-3:

testproc := proc(j, u, k) option remember;
   local calcvals;
   if not [j,u,k]::[numeric,numeric,posint] then
      return 'procname'(args);
   end if;
   calcvals:=doCalc(j,u);
   evalf(calcvals[k]);
end proc:

NN := 8;

8

(gridji,rngj,rngi):=[5,5],0..3,-2.0..-0.0;
PP[gridji,4]:=plot3d(max(-1e2,testproc(i,j,4)), j=rngj, i=rngi, grid=gridji):
PP[gridji,7]:= plot3d(Im(testproc(i,j,7)), j=rngj, i=rngi, grid=gridji):

[5, 5], 0 .. 3, -2.0 .. -0.

## Surface plot using 3d
PP[gridji,4]:

testproc(-1.0,1.5,4);

-119.9451877

## This is the pmin values
data[gridji,4]:=op([1,3],PP[gridji,4]);

data[[5, 5], 4] := Matrix(5, 5, {(1, 1) = -4.56109372700000, (1, 2) = -16.5517306200000, (1, 3) = -8.36945680000000, (1, 4) = -9.14002617600000, (1, 5) = 0.723410141700000e-2, (2, 1) = -4.40785476000000, (2, 2) = -14.6977719300000, (2, 3) = -7.86757454700000, (2, 4) = -100., (2, 5) = 0.723450034100000e-2, (3, 1) = -45.1423258900000, (3, 2) = -11.0010802700000, (3, 3) = -100., (3, 4) = -35.2069491100000, (3, 5) = -59.0582965200000, (4, 1) = -16.5967943100000, (4, 2) = -7.75166363700000, (4, 3) = -21.5341635200000, (4, 4) = -15.0367168400000, (4, 5) = -18.1717834500000, (5, 1) = -8.80334277800000, (5, 2) = -39.1054788600000, (5, 3) = -10.0221677400000, (5, 4) = -8.34415057400000, (5, 5) = -9.22785869200000})

## Convert Array to Matrix to allow the use listcontplot
ConMatrix := Matrix(data[gridji,4]);

ConMatrix := Matrix(5, 5, {(1, 1) = -4.56109372700000, (1, 2) = -16.5517306200000, (1, 3) = -8.36945680000000, (1, 4) = -9.14002617600000, (1, 5) = 0.723410141700000e-2, (2, 1) = -4.40785476000000, (2, 2) = -14.6977719300000, (2, 3) = -7.86757454700000, (2, 4) = -100., (2, 5) = 0.723450034100000e-2, (3, 1) = -45.1423258900000, (3, 2) = -11.0010802700000, (3, 3) = -100., (3, 4) = -35.2069491100000, (3, 5) = -59.0582965200000, (4, 1) = -16.5967943100000, (4, 2) = -7.75166363700000, (4, 3) = -21.5341635200000, (4, 4) = -15.0367168400000, (4, 5) = -18.1717834500000, (5, 1) = -8.80334277800000, (5, 2) = -39.1054788600000, (5, 3) = -10.0221677400000, (5, 4) = -8.34415057400000, (5, 5) = -9.22785869200000})

## Assign the min and max value of pmin to minv and maxv
(minv,maxv):=[min[defined],max[defined]](data[gridji,4])[];

-100., 0.723450034100000042e-2

## GRB color and the contours
(color1,color2):="Red","Yellow":
conts:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv]/2.0;

 

[-50.00000000, -44.44404253, -38.88808506, -33.33212758, -27.77617011, -22.22021264, -16.66425516, -11.10829770, -5.552340220, 0.3617250170e-2]

colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),
                                               i=1..N)])([ColorTools:-Color(color2)[]],
                                                         [ColorTools:-Color(color1)[]],
                                                         nops(conts))):

##

## Use the Conmatrix and the conts to plot contours and transform the axes
subsindets(plots:-listcontplot(ConMatrix, contours=conts,filledregions,thickness=0),
           specfunc(anything,THICKNESS),u->THICKNESS(0.2)):
LCP[gridji,4]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),
                                            lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%):
display(LCP[gridji,4],
        seq(plot(-2.0,1..1,legend=evalf[4](conts[i]),thickness=15,
                 color=colorlist[nops(conts)-i+1]),
            i=1..nops(conts)),
        size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

## This deals with the Imaginary part

## Surface plot using plot3d for Imginary part
PP[gridji,7]:

##

data[gridji,7]:=op([1,3],PP[gridji,7]):

## Convert Array to matrix
AMatrix := Matrix(data[gridji,7]);

AMatrix := Matrix(5, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = 0., (4, 2) = 0., (4, 3) = 0., (4, 4) = 0., (4, 5) = 0., (5, 1) = 0., (5, 2) = 0., (5, 3) = 0., (5, 4) = 0., (5, 5) = 0.})

(minv,maxv):=[min,max](data[gridji,7])[];

0., 0.

(color1,color2):="Red","Yellow":
conts1:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv]/~1.5;

[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]

colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),
                                               i=1..N)])([ColorTools:-Color(color2)[]],
                                                         [ColorTools:-Color(color1)[]],
                                                         nops(conts1))):

subsindets(plots:-listcontplot(AMatrix,
                               contours=conts1-~1e-9,filledregions,thickness=0),
           specfunc(anything,THICKNESS),u->THICKNESS(0.2)):
LCP[gridji,7]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),
                                            lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%):
display(LCP[gridji,7],
        seq(plot(-2.0,1..1,legend=evalf[4](conts1[i]),thickness=15,
                 color=colorlist[nops(conts)-i+1]),
            i=1..nops(conts1)),
        size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

 

``

Download Case3_Ijuptilk_Contourplot_pmin_290822_acc.mw

@ijuptilk You should provide a link to your revised worksheet with the changes to the x-axis range that your Question mentioned in a later revision.

Otherwise people might edit it in different places than you intend.

@Ronan If I make the target a valid .xlsx (Excel) file then on my Maple 2021.2 for Linux it opens the file in Libre Office when I click the link. And for a .png image file it opens in my machine's image viewer.

I don't know how your machine handles (or allows configuring of) file associations.

@Dkunb I'll have to wait until tomorrow to get a moment to answer your query about the details.

How long did your original grid=[25,25] contourplot take? How long does the above approach take? Is it adequate?

@Ronan Within a string a single backslash denotes that the next character will be escaped.

So use a double backslash (on MS-Windows) to get an effective backslash. (Or four, to get two, etc.)

That's not a Maple specific thing. (See here, here, etc, etc.)

I was just able to submit a new Question, and edit its text body.

(I've since deleted that test submission.)

@Carl Love Thanks.

Perhaps that means that somebody has edited the Question's body and defaced it by removing the original example.

I hope it wasn't the original poster. Deleting one's own question -- after anyone's spent time looking into it -- is very inconsiderate behaviour.

@Carl Love Where is the question?

@Rouben Rostamian  

That obtained solution procedure U only allows requests for arguments within the (earlier) specified range. Ie,

  U(-1e-5, 0.5);
    Error, (in solnproc) requested y value must be in the range HFloat(0.0)..HFloat(10.0)

And the call -D[1](U)(0,0.5) invokes fdiff to approximate the derivative numerically. It uses a simple symmetric differencing scheme with h=Float(1,-Digits/2), and thereupon attempts to call,
  U(-0.5e-5, 0.5)
which throws that same error. That error is caught and absorbed, in this situation.

With default working precision you could get this approximation,

  -D[1](U)(0.5e-5,0.5);
                    0.5443737835

Of course one is free to approximate the derivative manually, in some other fashion.

@MalakMMK You should show us all that you've done.

You can upload your worksheet using the green up-arrow in the Mapleprimes editor.

2 3 4 5 6 7 8 Last Page 4 of 494