acer

32395 Reputation

29 Badges

19 years, 344 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@salim-barzani I find your description of what you want to be confusing and incomplete. Please try and explain exactly what you want, with more explcit nouns and fewer vague pronouns.

You can get a 3D plot from your second example,

plot3d(G(l, a, 1, 1, 1, 1, 5), l=-2..2, a=0..1.5,
       colorscheme="turbo",
       style=surface, lightmodel=Light1,
       orientation=[-110,75,0]);

Or you could Explore the 2D plot, wrt parameter `a`.

Or you could get 2D contourplot or densityplot. And so on.

Bgraph_acc.mw

In your middle example (that shows error messages) your procedure G is malformed.

The mistake is that you have subscripted names entered using two different 2D Input syntax mechanisms. The procedure's parameters are defined as names with double-underscore (eg. alpha__3) while the body of the procedure has indexed names (eg. alpha[3]).

If you convert the 2D Input to 1D Input then the synax muddle becomes more visually apparent. Yes, there is a line of code to do assignments of one form to another, but that's fragile in the face of the way you do operator construction.

Bgraph1_ac.mw

You had kind of problem a few times before. I think that you'd be better off using just one of those two forms, consistently, throughout your worksheet, than mixing them. Especially if you use 2D Input syntax for operator definitions, etc.

To use Explore you just have to combine steps 4 and 5. The separation into steps 1-5 is not all necessary.

I upvote.

Are you considering putting this version onto the Application Center and/or the Maple Cloud?

@dharr 

[edit] I have made a correction in my code below, and have edited this Reply.

Using Maple 2024.2, the univariate numeric integration method=_d01amc appears to be able to handle directly the semi-finite integral (with your evalc'd form) faster than does method=_d01akc for the oscillatory integrand transformed to a finite range.

The plot3d with grid=[40,40] can then be done in  less than 3.3sec on my machine, in contrast to 9.6sec with that _d01akc.

I have not looked to see whether there was anything faster still.

restart;

rsol:=1/2/Pi^(1/2)*(1/4*x*Int(sin(t-zeta)/zeta^(7/2)*exp(-1/4*x^2/zeta)*(x^2
-6*zeta),zeta = 0 .. t, method=_Dexp)+Pi^(1/2)*(erf(1/2*(2*t+x)/t^(1/2))*exp(x+t)-erf(1/2*(2*
t-x)/t^(1/2))*exp(-x+t)-exp(x+t)+exp(-x+t))):

params := {x = 0.5, t = 0.1}:

ans_exact := evalf(eval(rsol, params));

.6581132195

V := -1/2*I*exp(k*x*I - k^2*t)*(1/(k - I) + 1/(k + I) - k*(1/(k^2 + I) + 1/(k^2 - I)))/Pi:
k1 := r*exp(1/8*I*Pi): k2 := r*exp(7/8*I*Pi): dk1 := exp(1/8*I*Pi): dk2 := exp(7/8*I*Pi):
integrand1 := eval(V, k = k1)*dk1 - eval(V, k = k2)*dk2:

integrand2 := simplify(evalc(integrand1)):

forget(`evalf/int`); forget(evalf);
Fokas_int2 := CodeTools:-Usage(evalf(Int(eval(integrand2, params), r = 0 .. infinity, method=_d01amc)));
Fokas2 := evalf(eval(exp(-x/sqrt(2))*cos(t - x/sqrt(2)), params)) + Fokas_int2;

memory used=139.68KiB, alloc change=0 bytes, cpu time=3.00ms, real time=2.00ms, gc time=0ns

-0.2162433743e-1

.6581132200

approx_u2 := proc(x,t) local temp;
  temp := Int(eval(integrand2, [:-x=x,:-t=t]), r = 0 .. infinity, method=_d01amc);
  evalf(eval(exp(-x/sqrt(2))*cos(t - x/sqrt(2)), [:-x=x,:-t=t]) + temp);
end proc:

approx_u2(0.5,0.1);

.6581132200

forget(`evalf/int`); forget(evalf);
CodeTools:-Usage(plot3d(approx_u2, 0 .. 3, 0 .. 2*Pi, grid = [40, 40], axes = boxed, labels = ["x", "t", "u(x,t)"], title = "Fokas Method of solution of Half-Line Heat Equation", shading = zhue, size=[300,300]));

memory used=289.58MiB, alloc change=36.00MiB, cpu time=3.52s, real time=3.31s, gc time=365.05ms

asympt(integrand2, r):

Suggesting the following change of variables to u ranging from 0 to 1

itr := u = 1 - exp(-r^2*t/sqrt(2)):
eval(itr,r=0), eval(itr,r=infinity) assuming t>0:
tr := r = solve(itr, r)[1]:

int3 := PDEtools:-dchange(tr,Int(integrand2, r = 0..infinity), [u], simplify) assuming t>0:

integrand3 := IntegrationTools:-GetIntegrand(int3):
#Fokas_int3 := t/Pi*Int(integrand3, u = 0..1, method = _d01akc):

forget(`evalf/int`); forget(evalf);
Fokas_int3 := CodeTools:-Usage(evalf(eval(t/Pi*Int(integrand3, u = 0..1, method = _d01akc), params)));
evalf(eval(exp(-x/sqrt(2))*cos(t - x/sqrt(2)), params)) + Fokas_int3;

memory used=152.34KiB, alloc change=0 bytes, cpu time=6.00ms, real time=6.00ms, gc time=0ns

-0.2162433742e-1

.6581132200

approx_u := unapply(exp(-x/sqrt(2))*cos(t - x/sqrt(2)) + t/Pi*Int(integrand3, u = 0..1, method = _d01akc), x,t, numeric):

evalf(approx_u(0.5,0.1));

.6581132200

forget(`evalf/int`); forget(evalf);
CodeTools:-Usage(plot3d(approx_u, 0 .. 3, 0 .. 2*Pi, grid = [40, 40], axes = boxed, labels = ["x", "t", "u(x,t)"], title = "Fokas Method of solution of Half-Line Heat Equation", shading = zhue, size=[300,300]));

memory used=268.52MiB, alloc change=0 bytes, cpu time=9.82s, real time=9.59s, gc time=377.55ms

forget(`evalf/int`); forget(evalf); CodeTools:-Usage(plot3d(proc (a, b) options operator, arrow; approx_u(a, b)-approx_u2(a, b) end proc, 0 .. 3, 0.1e-1 .. 2*Pi, grid = [40, 40], axes = boxed, labels = ["x", "t", ""], size = [300, 300]))

memory used=407.28MiB, alloc change=0 bytes, cpu time=12.17s, real time=11.71s, gc time=730.92ms

Download Fokas_acc.mw

@sand15 Here's another variant, for Maple 2015.

I showed two uses with value, just to illustrate that while only one of them might display the nicer way for output, both could be used to compute with in some ways.

Yes, I realize that some variants produce that particular error (evaluation of an rtable indexed by a name rather than integer value). That's indeed why I looked again for Maple2015-specific workarounds...

restart;

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Dec 20 2015, Build ID 1097895`

with(Statistics):

 

Dist := proc(N::{posint, name}, alpha::uneval)
  local K;
  uses Statistics;

  K := numelems(eval(alpha));
  
  Distribution( ':-PDF' = (z -> Product((''alpha''[k]+z[k])^N, k=1..K)) )
end proc:

 

L := [3,17,22];
X2 := RandomVariable(Dist(N, L)):

[3, 17, 22]

res2 := PDF(X2, y);

Product((L[k]+y[k])^N, k = 1 .. 3)

value(res2);

(3+y[1])^N*(17+y[2])^N*(22+y[3])^N

value(eval(res2,1));

(3+y[1])^N*(17+y[2])^N*(22+y[3])^N

eval(res2,1);

Product((L[k]+y[k])^N, k = 1 .. 3)

value(eval(res2, L=[1/2,1/3,1/4]));

(1/2+y[1])^N*(1/3+y[2])^N*(1/4+y[3])^N


Download Product_error_ac_2015_2.mw

Do you mean like with a legend, or caption, or inlined textplot, or...?

How do you want to "link" them?

Did it work for you?

@Muhammad Usman 

restart;

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Dec 20 2015, Build ID 1097895`

with(plots): with(plottools,polygon): with(ColorTools,Color):

aa := 4: bb := 2:  

f := -(x^2/(aa^2)+y^2/(bb^(2))-1) * aa^2*bb^2/(aa^2+bb^2);

-(1/5)*x^2-(4/5)*y^2+16/5

N := 250:
PD := densityplot(piecewise(x<sqrt(-4*y^2+4),f,undefined),
            x = 0 .. aa/2, y = 0 .. bb/2, style=surface,
            colorscheme=["zgradient",["Blue","Red"],'colorspace'="HSV"],
            axes = boxed, title = "Contour plot over quarter ellipse",
            grid = [N,N]):

PC := contourplot(f, x = 0 .. sqrt(-4*y^2+4), y = 0 .. bb/(2),
            contours = 15, color=black, thickness=2,
            axes = boxed, grid = [100,100]):

BND := plot([sqrt(-4*y^2+4),y,y=0..bb/2],thickness=3,color=black):

Q := plot3d(f, x = 0 .. sqrt(-4*y^2+4), y = 0 .. bb/(2)):

(zmin,zmax) := (min,max)(op([1,1],Q)[..,..,3]):

display(PD,PC,BND,
  seq(polygon([[0,0],[0,0]],'thickness'=0.3, 'legend'=sprintf("%5.3f",z),
              'color'=Color("HSV",[0.6666666667*(zmax-z)/(zmax-zmin),1,1])),
      z=zmax..zmin, (zmin-zmax)/(16)),
  'legendstyle'=['location'=':-right'], 'size'=[580,500]);

foo := copy(op([1,4,2],PD)):
foo[..,..,1] := 360*foo[..,..,1]:
foo := ImageTools:-HSVtoRGB(foo):
for i from 1 to N do
  for j from 1 to N do
    if foo[i,j,1]=0 and foo[i,j,2]=0 and foo[i,j,3]=0 then
      foo[i,j,1],foo[i,j,2],foo[i,j,3] := 1,1,1;
    end if;
  end do;
end do:
PDnew := subsop([1,4]=COLOR(RGB,foo),PD):

display(PDnew,PC,BND,
  seq(polygon([[0,0],[0,0]],'legend'=sprintf("%5.3f",z),
              'color'=Color("HSV",[0.6666666667*(zmax-z)/(zmax-zmin),1,1])),
      z=zmax..zmin, (zmin-zmax)/(16)),
  'legendstyle'=['location'=':-right'], 'size'=[580,500]);

Download Usman_contplot3.mw

@Muhammad Usman Why didn't you mention all your requirements when you first posted the Question?

The contourplot command of Maple 2015 doesn't support the colorscheme option, in the way that Maple 2025 does.

So it's a bit tricky to get your [blue,green,yellow,red] color spread for the inner part of the 2D plot.

One way is to use densityplot rather than filled with contourplot.

Also, in Maple 2015 I don't think that a split-by-value is possible for doing densityplot with a colorscheme. (I do some tricks for it below, because I didn't want to mess with a color-function since it'd mean I'd have to force the choices of contour-values...)

But here is one way to get that.

restart;

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Dec 20 2015, Build ID 1097895`

with(plots):

aa := 4: bb := 2:  

f := -(x^2/(aa^2)+y^2/(bb^(2))-1) * aa^2*bb^2/(aa^2+bb^2);

-(1/5)*x^2-(4/5)*y^2+16/5

N := 250:
PD := densityplot(piecewise(x<sqrt(-4*y^2+4),f,undefined),
            x = 0 .. aa/2, y = 0 .. bb/2, style=surface,
            colorscheme=["zgradient",[blue,green,yellow,red]],
            axes = boxed, title = "Contour plot over quarter ellipse",
            grid = [N,N]):

PC := contourplot(f, x = 0 .. sqrt(-4*y^2+4), y = 0 .. bb/(2),
            contours = 15, color=black, axes = boxed, grid = [100,100]):

BND := plot([sqrt(-4*y^2+4),y,y=0..bb/2],thickness=2,color=black):

display(PD,PC,BND);

foo := op([1,4,2],PD):
for i from 1 to N do
  for j from 1 to N do
    if foo[i,j,1]=0 and foo[i,j,2]=0 and foo[i,j,3]=0 then
      foo[i,j,1],foo[i,j,2],foo[i,j,3] := 1,1,1;
    end if;
  end do;
end do:

display(PD,PC,BND);

Download Usman_contplot2.mw

Have you submitted it as a bug report?

While you might not want to use a specific forced method as workaround, I notice that the error gets avoided (due to some internal remember table, I suppose) with a second attempt:

restart;
kernelopts(version);
           Maple 2025.1, X86 64 LINUX, Jun 12 2025, Build ID 1932578

> integrand:=1/2/x^(9/2)*2^(1/2)*Pi^(1/2)/(1/x)^(1/2)*cos(1/x):

> int(integrand,x); int(integrand,x);

Error, (in convert/sincos) too many levels of recursion

                 1/2   1/2                           2
          1/4 I 2    Pi    (4 I x cos(1/x) - 2 I (2 x  - 1) sin(1/x))
          -----------------------------------------------------------
                                  5/2      1/2
                                 x    (1/x)

@WD0HHU On July 18, 2025, Tech Support asked you whether your status bar was enabled, and showed a screenshot of its control in the View ribbon.

@janhardo Btw, it's also bit heavy-handed to do actions like,
    combine(b, trig)
or simplify. The problem is that you don't allow the user an optional way to choose not to do it.

I can make a guess as to one reason why you might have put it in there: Your SplitPDE1 procedure applies expand, to separate the linear/nonlinear terms as separate additive terms (so that they can be select'd). But a factored product, unexpanded, could contain linear terms that get missed, eg.     u_t*(K+u_x)
So, if you inadvertantly expand'd/blew-up some trig terms then you might want to combine to get back politely to the original form. However, one problem is that you cannot know whether that has to be done, and you don't offer a way to not do it.

That's why earlier I mentioned applying normal (or a frontend'd expand) instead of just a plain expand call. That leaves trig stuff like sin(u_t+K+x) alone, and you don't need to take on the burden of possible repair. But normal will still cause the wanted expansion of factored products, so that the linear & nonlinear addends get separated. Eg, compare,
  expr := sin(diff(u(x),x)+K+x) + u(x)*(K+diff(u(x),x));
  expand(expr);
  normal(expr);

Hope that makes sense.

@salim-barzani I've made another edit/correction in my code above (both places).

@janhardo K is constant w.r.t. all of x,y,t, functions u,v of those, and derivatives of those functions.

Being unrelated to all of those, it's importance is not much different from a number here.

1 2 3 4 5 6 7 Last Page 1 of 592