acer

31431 Reputation

29 Badges

19 years, 134 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@C_R If you would really like to utilize an expression form calling sequence (without the need for uneval quotes each time you reference it) such as,
   plot( fr(t), t=0..1)
rather than an operator form calling sequence like,
   plot( fr, 0..1)
then you could instead construct fr in one of these ways:

   R := unapply('eval(r,dsol1(t))',t,numeric):

   R := proc(t)
       if not t::numeric then return 'procname'(t); end if;
       eval(r,dsol1(t));
   end proc:

The effect there is that both of those return fr(t) as the unevaluated function call itself when argument t is not of type numeric (eg. here, a symbol).

I suppose you already know that the numeric solving of this system faces some quirks & difficulties. Perhaps you are trying to investigate the numeric DE solver itself. If not, then without details it's hard to know whether you'd run into situations in which you might 1) want to specify stepsize or another option to dsolve, 2) augment the system if looking for integrals|derivatives versus post-calculation via evalf@D or evalf@Int, etc, 3) other...

@bathudaide Feel free to provide a URL for the author's code site.

(You asked before, "Can I upload it for discussion?", where "it" really does seem to mean "source code" from the previous sentence. So, no, you can't upload a copy of his source code here without permission. But you can provide a link to it.)

@bathudaide You asked, "I have found the source code for SDPTools.  Can I upload it for discussion?"

Only with the explicit permission of its author (or copyright holder).

I find the business of having deleted replies stick around distracting and unhelpful.

I also notice that they remain even if the response was deleted as spam (using the special item for that). So some Posts that attract multiple spammers get quite cluttered.

I think that those changes might also have coincided with the regression/bug in which responses without votes are ordered wrong by date. Such used to be ordered downwards with oldest higher, and now they are (poorly) ordered in the opposite direction.

@bathudaide Yes, that's what I mean about those two files.

And I get an equivalent error message, on Linux.

I added some update notes to my Comment above.

I notice that for Linux the name of the shared-binary in that site's FGb download is named "libfgbuni.so". And that webpage is 9 years old. And modern Maple for Linux itself already has its own "libfgbuni.so" shared-binary (from the same research people, IIANM).

So it's possible that for MS-Windows you might be already have a functioning .dll that exports the FGb external functions, in your installed Maple, and be able to simply unpack the Linux FGb download and install just its Maple Library portions. (It's so old that it's bundled as the pair FGblib.ind and FGblib.lib) And then augment libname as it instructs, by prepending it with the location you put those Maple Library files.

You'll know if it's gotten off the ground because this example would work:
   with(FGb):
   fgb_gbasis([x^2+y,y^2+x],0,[x],[y]);

[update] I now suspect that the above might not work at all, sorry!

It's not clear that all the functions' calling sequences have stayed the same between that old FGb package of 2015 on their site and the current sharEd-binary in modern Maple. When I try to utilize the .ind|.lib files from their site with even the Linux libfgbuni.so shared-binary then with(FGb) throws an error (path issues, perhaps in its ModuleLoad).

Also, the calling sequences of the exports in Maple 2024's Library's bundled fgbrs package seem different from what that old 2015 FGb package has. Eg.

restart;
FGb := fgbrs:
with(FGb):
fgb_gbasis([x^2+y,y^2+x],0,[x],[y]);
Error, invalid input: fgbrs:-fgb_gbasis expects its 2nd argument, vars1, to be of type list(name), but received 0
restart;
fgbrs:-fgb_gbasis([x^2+y,y^2+x],0,[x],[y]);
Error, invalid input: fgbrs:-fgb_gbasis expects its 2nd argument, vars1, to be of type list(name), but received 0

I'm out of ideas, sorry.

49 - 25 + 9 = 33

@Scot Gould I don't know whether the underlying problem would just get fixed, without a new Options/interface setting.

If you want both modes (now and then) then it might be easier to utilize an alternate OS user account for just the PDF export, rather than keep shunting versions of that preferences file.

@Andiguys You asked, "Can we change x-axis label to Ce from x and y label  to Cv?"

Yes, you could simply add the option,

   labels = ["Ce", "Cv"]

to the call to the procedure.

Question_plot_cv_ac2.mw

@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).

First 9 10 11 12 13 14 15 Last Page 11 of 577