acer

31297 Reputation

29 Badges

19 years, 110 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@salim-barzani You could add some additional constraints.

(I've also declared gamma as local at the top-level, which has been mentioned to you a few times before. And I've used PDEtools:-dcoeffs instead of your copy&paste of the collect result, in order to get the coefficient expressions that you want equated to zero. And I've changed some subs calls to eval. Those changes don't change the final result, but they make the code less dodgy.)

For example,

restart

_local(gamma)

PDEtools:-undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

Fode := (-delta*eta^2+alpha*eta)*(diff(diff(U(xi), xi), xi))-U(xi)*(2*eta*gamma*theta*(delta*eta-alpha)*U(xi)^2+eta^2*delta*k^2+(-alpha*k^2-2*delta*k)*eta+2*k*alpha+delta) = 0

F := sum(a[i]*G(xi)^i, i = 0 .. 1)

a[0]+a[1]*G(xi)

D1 := diff(F, xi)

S := (diff(G(xi), xi))^2 = G(xi)^4+A[2]*G(xi)^2+A[1]

S1 := diff(G(xi), xi) = sqrt(G(xi)^4+A[2]*G(xi)^2+A[1])

E1 := eval(D1, S1)

a[1]*(G(xi)^4+A[2]*G(xi)^2+A[1])^(1/2)

D2 := diff(E1, xi)

E2 := eval(D2, S1)

(1/2)*a[1]*(4*G(xi)^3*(G(xi)^4+A[2]*G(xi)^2+A[1])^(1/2)+2*A[2]*G(xi)*(G(xi)^4+A[2]*G(xi)^2+A[1])^(1/2))/(G(xi)^4+A[2]*G(xi)^2+A[1])^(1/2)

K := U(xi) = F

K1 := diff(U(xi), xi) = E1

K2 := diff(U(xi), xi, xi) = E2

L := eval(Fode, {K, K1, K2})

L1 := normal(L)

-2*G(xi)^3*delta*eta^2*gamma*theta*a[1]^3+2*G(xi)^3*alpha*eta*gamma*theta*a[1]^3-6*G(xi)^2*delta*eta^2*gamma*theta*a[0]*a[1]^2+6*G(xi)^2*alpha*eta*gamma*theta*a[0]*a[1]^2-6*G(xi)*delta*eta^2*gamma*theta*a[0]^2*a[1]+6*G(xi)*alpha*eta*gamma*theta*a[0]^2*a[1]-2*delta*eta^2*gamma*theta*a[0]^3-2*G(xi)^3*delta*eta^2*a[1]-G(xi)*delta*eta^2*k^2*a[1]+2*alpha*eta*gamma*theta*a[0]^3+2*G(xi)^3*alpha*eta*a[1]+G(xi)*alpha*eta*k^2*a[1]-G(xi)*delta*eta^2*A[2]*a[1]-delta*eta^2*k^2*a[0]+G(xi)*alpha*eta*A[2]*a[1]+2*G(xi)*delta*eta*k*a[1]+alpha*eta*k^2*a[0]-2*G(xi)*alpha*k*a[1]+2*delta*eta*k*a[0]-G(xi)*delta*a[1]-2*alpha*k*a[0]-delta*a[0] = 0

extra := {eta <> 0, a[1] <> 0}

{eta <> 0, a[1] <> 0}

eqs := {`~`[`=`](PDEtools:-dcoeffs((lhs-rhs)(L1), G(xi)), 0)}

{-6*delta*eta^2*gamma*theta*a[0]*a[1]^2+6*alpha*eta*gamma*theta*a[0]*a[1]^2 = 0, -2*delta*eta^2*gamma*theta*a[1]^3+2*alpha*eta*gamma*theta*a[1]^3-2*delta*eta^2*a[1]+2*alpha*eta*a[1] = 0, -2*delta*eta^2*gamma*theta*a[0]^3+2*alpha*eta*gamma*theta*a[0]^3-delta*eta^2*k^2*a[0]+alpha*eta*k^2*a[0]+2*delta*eta*k*a[0]-2*alpha*k*a[0]-delta*a[0] = 0, -6*delta*eta^2*gamma*theta*a[0]^2*a[1]+6*alpha*eta*gamma*theta*a[0]^2*a[1]-delta*eta^2*k^2*a[1]+alpha*eta*k^2*a[1]-delta*eta^2*A[2]*a[1]+alpha*eta*A[2]*a[1]+2*delta*eta*k*a[1]-2*alpha*k*a[1]-delta*a[1] = 0}

COEFFS := solve(`union`(eqs, extra), {alpha, eta, a[0], a[1]}, explicit)

COEFFS := map[2](remove, evalb, [COEFFS])

[{alpha = delta*(eta^2*k^2+eta^2*A[2]-2*eta*k+1)/(eta*k^2+eta*A[2]-2*k), a[0] = 0, a[1] = -1/(-gamma*theta)^(1/2)}, {alpha = delta*(eta^2*k^2+eta^2*A[2]-2*eta*k+1)/(eta*k^2+eta*A[2]-2*k), a[0] = 0, a[1] = 1/(-gamma*theta)^(1/2)}]

verifiedsols := table([seq(COEFFS[i]={},i=1..nops(COEFFS))]):
for thisCF in COEFFS do
  thisFode := simplify(eval(Fode, thisCF));
  if thisFode = (0=0) then
    verifiedsols[thisCF] := {U(xi)=U(xi)};
  else
    S2 := dsolve(S, G(xi));
    for thisS2 in S2 do
      thisCFS2 := simplify(eval(eval(K, thisCF), thisS2));
      if odetest(thisCFS2, thisFode)=0 then
        verifiedsols[thisCF] := verifiedsols[thisCF] union {thisCFS2};
      end if;
    end do;
  end if;
end do:

#map(u->print(u=verifiedsols[u]), [indices(verifiedsols,'nolist')]):

GG := [indices(verifiedsols,'nolist')];

[{alpha = delta*(eta^2*k^2+eta^2*A[2]-2*eta*k+1)/(eta*k^2+eta*A[2]-2*k), a[0] = 0, a[1] = -1/(-gamma*theta)^(1/2)}, {alpha = delta*(eta^2*k^2+eta^2*A[2]-2*eta*k+1)/(eta*k^2+eta*A[2]-2*k), a[0] = 0, a[1] = 1/(-gamma*theta)^(1/2)}]

GG[1], verifiedsols[GG[1]];
simplify(eval(Fode,GG[1]));
if % <> (0=0) then
  odetest(verifiedsols[GG[1]], eval(Fode,GG[1]));
end if;

{alpha = delta*(eta^2*k^2+eta^2*A[2]-2*eta*k+1)/(eta*k^2+eta*A[2]-2*k), a[0] = 0, a[1] = -1/(-gamma*theta)^(1/2)}, {U(xi) = -JacobiSN((1/2)*(2*(A[2]^2-4*A[1])^(1/2)-2*A[2])^(1/2)*xi+c__1, 2^(1/2)*(A[2]^2*A[1]-A[1]*A[2]*(A[2]^2-4*A[1])^(1/2)-2*A[1]^2)^(1/2)/(A[2]*(A[2]^2-4*A[1])^(1/2)-A[2]^2+2*A[1]))*A[1]*2^(1/2)/((-gamma*theta)^(1/2)*(A[1]*(-A[2]+(A[2]^2-4*A[1])^(1/2)))^(1/2))}

delta*eta*(2*U(xi)^3*gamma*theta-U(xi)*A[2]+diff(diff(U(xi), xi), xi))/((k^2+A[2])*eta-2*k) = 0

0

dsolve(simplify(eval(Fode,GG[1])), U(xi));
odetest(%, eval(Fode,GG[1]));

U(xi) = c__2*(A[2]/(gamma*theta*c__2^2-gamma*theta+A[2]))^(1/2)*JacobiSN(((gamma*theta-A[2])^(1/2)*xi+c__1)*(A[2]/(gamma*theta*c__2^2-gamma*theta+A[2]))^(1/2), c__2*(-(gamma*theta-A[2])*gamma*theta)^(1/2)/(gamma*theta-A[2]))

0

GG[2], verifiedsols[GG[2]];
simplify(eval(Fode,GG[2]));
if % <> (0=0) then
  odetest(verifiedsols[GG[2]], eval(Fode,GG[2]));
end if;

{alpha = delta*(eta^2*k^2+eta^2*A[2]-2*eta*k+1)/(eta*k^2+eta*A[2]-2*k), a[0] = 0, a[1] = 1/(-gamma*theta)^(1/2)}, {U(xi) = JacobiSN((1/2)*(2*(A[2]^2-4*A[1])^(1/2)-2*A[2])^(1/2)*xi+c__1, 2^(1/2)*(A[2]^2*A[1]-A[1]*A[2]*(A[2]^2-4*A[1])^(1/2)-2*A[1]^2)^(1/2)/(A[2]*(A[2]^2-4*A[1])^(1/2)-A[2]^2+2*A[1]))*A[1]*2^(1/2)/((-gamma*theta)^(1/2)*(A[1]*(-A[2]+(A[2]^2-4*A[1])^(1/2)))^(1/2))}

delta*eta*(2*U(xi)^3*gamma*theta-U(xi)*A[2]+diff(diff(U(xi), xi), xi))/((k^2+A[2])*eta-2*k) = 0

0

dsolve(simplify(eval(Fode,GG[2])), U(xi));
odetest(%, eval(Fode,GG[2]));

U(xi) = c__2*(A[2]/(gamma*theta*c__2^2-gamma*theta+A[2]))^(1/2)*JacobiSN(((gamma*theta-A[2])^(1/2)*xi+c__1)*(A[2]/(gamma*theta*c__2^2-gamma*theta+A[2]))^(1/2), c__2*(-(gamma*theta-A[2])*gamma*theta)^(1/2)/(gamma*theta-A[2]))

0

Download odetest_acc.mw
 

Your phrasing, "...in your case why calculate seperatly not all our parameter toghether map must be [4]?" doesn't make sense to me, sorry.

@Andiguys I don't know why you state that you included Kitonum's code, since your worksheet is completely missing the code for the procedure that does the contour labels.

Also, there are several other serious mistakes. You didn't actually store the results in a Matrix M, yet you then try to use the unassigned name M as if it were such a Matrix -- to access the unstored results. And you incorrectly attempt to call the labelling procedure with the result from contourplot instead of the actual expression as first the argument. And so on.

In the worksheet attached here I have put the code for the procedure that constructs this final plot (with labels on the curves) in the Startup Code of the worksheet. If you want to use it in another worksheet then you'll have to copy it over (or functional equivalent); merely calling it (as you tried before) cannot work, unless you also make the code accessible.

Not all your code highlighting will render well when the sheet is inlined on this site.

restart

`&pi;m` := proc (e, delta) options operator, arrow; b*y-(1-delta-e-tau01)*y*Cv-h*delta*y-Ce*k2*y-Clm*e^2-Cex*tau01*y-Am*e*y-Rm*e*y+R0m*e^2*y^2-piecewise(0 < s2*y*(delta0-e-delta), s2*y*(delta0-e-delta), 0)-piecewise(0 < s1*(tau0-e-tau01-k2)*y, s1*(tau0-e-tau01-k2)*y, 0)-g1*e^2*y^2-g2*delta^2*y^2 end proc

`&Pi;m` := eval(`&pi;m`(e, delta), [h = (1/2)*l+(1/2)*Cv+(1/2)*s2, k2 = -y*(Rer-Ce-l+Aer)/(2*(-R0er*y^2+Clr))])

DATA := [b = 150, Clm = 20000, Am = .5, Clr = 30000, Aer = 3.5, Rm = 48, Rer = 30, R0m = 5*10^(-5), R0er = 10^(-5), s1 = 4, s2 = 10, g1 = 0.5e-1, g2 = 0.9e-1, y = 2400, tau0 = .4, l = 40, delta0 = .25, tau01 = 0.5e-1, Cex = 5]

TRC := proc (e, delta) options operator, arrow; eval(`&Pi;m`, DATA) end proc

t := {80, 85, 95, 100}; ts := {1, 2, 3, 4, 5, 6, 7}

M := Matrix(nops(t)*nops(ts), 3); rr := 0; for Cv in t do for Ce in ts do s := Optimization:-Maximize(TRC(e, delta), e = 0 .. 1, delta = 0 .. 1, assume = nonnegative); stemp := s[1]; etemp := s[2][1]; deltatemp := s[2][2]; rr := rr+1; M[rr, 1 .. 3] := `<|>`(Cv, Ce, stemp) end do end do

R := Array(ArrayTools:-Reshape(M,[7,4,3]),datatype=float[8]):

func := Interpolation:-SplineInterpolation([[$1..7],[80,85,95,100]],R[..,..,3]):

conts := [seq(min(R[..,..,3])..max(R[..,..,3]),(max(R[..,..,3])-min(R[..,..,3]))/8)];

ContoursWithLabels(func(x, y), x = 1 .. 7, y = 80 .. 100, contours = conts, Coloring = [colorstyle = HUE, colorscheme = ["Blue", "Gold"], style = surface], TextOptions = [font = [HELVETICA, BOLD, 9], color = black], GraphicOptions = [thickness = 0], ImplicitplotOptions = [gridrefine = 2], size = [700, 600])

Download Question_plot_cv_ac.mw

@Andiguys Why don't you ask state all your requirements when you first post your Question?

You can use Kitonum's code for 2D contour plotting with labelled curves here. (You could also easily find that by doing a simple search on this site.) You can use that code it with the interpolating function in my earlier response above.

@Andiguys 

You have two choices, as far as using an external file goes and how you can reference it in the string passed to the ImportMatrix command. 1) use an absolute/full path, or 2) use a path relative to the currentdir(...) setting.

On Windows the choice 1) absolute path would look like, say,
  "C://Users/Andiguys/foo/data.txt"
or wherever you have actually put the file on your computer.

Suppose that the file is at "C://Users/Andiguys/foo/data.txt" and that currentdir(...) is set to, say,
   "C://Users/Andiguys"
then for the choice 2) the relative path would be just, "foo/data.txt" .

In my example I had currentdir(...) already set to the very folder/directory in which the file resided. So the relative path, for me, was just the filename by itself.

You don't have to use an external file. Your data could just start out life in a Matrix, say. I can't tell you how to adjust your code for that since you didn't give us your code.

@Andiguys I'm not entirely sure what you mean by, "...a 2D graph with x and y axes, where the z value is displayed along the curve at specific points."

I am supposing that you mean a 2D contour plot.

I'll do something a little different below, and construct an interpolating function from your data.

That can be used to obtain the same 3D surface as before. But it can also be used to produce a 2D contour plot.

In fact, it can be used like any other "function" of two variables, for plotting or polling pairs of values in your data's independent ranges. So even if my supposition about what you now want is mistaken then the interpolant can likely be used to attain the goal.

restart;

with(plots):

M := ImportMatrix("data.txt", source=csv)[..,2..4]:

R := Array(ArrayTools:-Reshape(M,[7,4,3]),datatype=float[8]):


Create an interpolating function from the data.

func := Interpolation:-SplineInterpolation([[$1..7],[80,85,95,100]],R[..,..,3]):


Show that the interpolant agrees with the surfdata surface.

display(
  plot3d(func(x,y),y=80..100,x=1..7,style=point,color=blue,
         grid=[5,5],symbolsize=15,symbol=solidcircle),
  surfdata(R, color=green)
);

conts := [seq(min(R[..,..,3])..max(R[..,..,3]),0.05)];

[HFloat(1.364583951), HFloat(1.414583951), HFloat(1.464583951), HFloat(1.5145839509999999), HFloat(1.564583951), HFloat(1.614583951), HFloat(1.664583951), HFloat(1.7145839509999998), HFloat(1.764583951)]

CP2 := contourplot(func(x,y), x=1..7, y=80..100,
                   contours=conts, coloring=["Blue","Gold"]):

CP2;


Show that the 2D contour plot agrees with the contours of the 3D surface.

display(
  plot3d(func(x,y),x=1..7,y=80..100,
         style=contour,contours=conts,
         colorscheme=["zgradient",["Blue","Gold"]]),
  plottools:-transform((x,y)->[x,y,1.3])(CP2),
  thickness=4
);


some_data_plot_3.mw

@C_R How is the second approach the best? The first is much faster and simpler, and completely avoids the float-roundoff (Pi) edge idiosyncrasy & need for higher grid values.

Cannot the individual portions all be exported to STL just as easily as the conjoined result from the second approach? I would expect that they can (even if done separately).

@C_R The first approach seems better all round. Why would you use the second approach?

The two kinds of gap can be reduced by increasing the corresponding grid option values.

You may be able to get by there with NULL instead of undefined (I had other variants where it had an effect). The key bit is adaptmesh=false.

But the first way is so much simpler, and faster and lightweight, that I don't understand any preference for using the second.

While it's not nearly as fast and simple as the above, the original methodology can be hammered into something that omits those conjoining artefacts.
1) Using adaptmesh=false to prevent plot3d from trying extra hard to stitch together pieces adjacent in the data structure.
2) split by undefined rather than NULL
3) increase the grid resolution.

restart;

rmod:= (t,n)->(round@((x->x)-floor))(n*t):

f:=(t,p)-> ifelse(is(rmod(t,4/(2*Pi))+rmod(p,4/(2*Pi))=1),1,NULL):

p1:=plot3d(f,0..2*Pi,0..Pi,coords=spherical,scaling=constrained,
           color="black",grid=[200,200],style=surface,adaptmesh=false):

p2:=plots:-sphereplot(0.999,theta=0..2*Pi,phi=0..Pi,color=white,style=surface):

plots:-display(p1,p2);


Download Chessboard_sphere_ac2.mw

Maple automatically transforms 1/2^(1/4) and 2^(-1/4) into 1/2*2^(3/4) .

That happens as a result of what's called automatic simplification, and cannot be avoid just by unevaluation. That is, it'd also happen for '1/2^(1/4)', ie. wrapped in single right-quotes (aka uneval quotes).

It should be possible programatically to turn it into an inert form and have that be displayed without the automatic simplification. Then you could get the output rendered as follows (or with 2^(-1/4) in the numerator):

restart;

with(InertForm):

'1/2^(1/4)';

(1/2)*2^(3/4)

'2^(-1/4)';

(1/2)*2^(3/4)

foo := something / 2 %^ (1/4):

value(foo);

(1/2)*something*2^(3/4)

Display(foo,'inert'=false);

something/Typesetting:-_Hold([`%^`(2, 1/4)])


Download inert_p_ex.mw

ps. Are you also wanting to simplify wrt that B^2 inside the square-root? Is B real?
Also, you have a term like e^(i*psi(x,t)), so is that deliberate that you don't use I and exp?

@salim-barzani Right the simplify command will not produce the results you want in Maple 2021.

That's why I showed you how it might be done using applyrule.

Is it the performance of the Maple GUI which is slow, or computation of results that is slow (even when run in both systems with all display of output results suppressed, etc)?

If it's the computation, then what computation?

@salim-barzani I've already showed you (above) that your given examples can be handled automatically by Maple 2024 by just the simplify command.

Such simplifications can also usually be restricted (eg. simplify(...,trig), etc) so that they don't alter other type of term.

You ask whether all your other Maple 2021 code (which you don't provide here) will run the ok in Maple 2024. It's likely that it will. But nobody can make such a prediction with 100% guarantee, without seeing the code.

Please put your followup examples on this here, instead of in wholly separate new Question threads,

@Ronan I added a response to your earlier Question thread on rotated text in a 2D plot.

@erik10 A terser variant on that is something like,

  Circle1 := ifelse(CircleOrNot="yes", circle(...), NULL);

or, when CircleOrNot is definitely assigned one of true or false,

  Circle1 := ifelse(CircleOrNot, circle(...), NULL);

First 7 8 9 10 11 12 13 Last Page 9 of 574