acer

32395 Reputation

29 Badges

19 years, 343 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

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.

@janhardo Your approach seems to need some more corrections. What do you think?

restart

with(PDEtools)

undeclare(prime, quiet); declare(u(x, y, t), quiet); declare(v(x, y, t), quiet); declare(f(x, y, t), quiet)

 

pde1 := diff(u(x, y, t), t)+diff(u(x, y, t), `$`(x, 3))+diff(u(x, y, t)*v(x, y, t), x)+t*(diff(u(x, y, t), t))+K

diff(u(x, y, t), t)+diff(diff(diff(u(x, y, t), x), x), x)+(diff(u(x, y, t), x))*v(x, y, t)+u(x, y, t)*(diff(v(x, y, t), x))+t*(diff(u(x, y, t), t))+K

pde2 := diff(u(x, y, t), x)-(diff(v(x, y, t), y))-y*(diff(u(x, y, t), y))-t

diff(u(x, y, t), x)-(diff(v(x, y, t), y))-y*(diff(u(x, y, t), y))-t

T := 'T'; linpred := proc (term) options operator, arrow; not has(simplify((eval(term, {u(x, y, t) = T*u(x, y, t), v(x, y, t) = T*v(x, y, t)}))/T), T) end proc

pde1_lin, pde1_nonlin := selectremove(linpred, expand(pde1)); pde2_lin, pde2_nonlin := selectremove(linpred, expand(pde2))

pde1_lin; pde1_nonlin

diff(u(x, y, t), t)+diff(diff(diff(u(x, y, t), x), x), x)+t*(diff(u(x, y, t), t))

(diff(u(x, y, t), x))*v(x, y, t)+u(x, y, t)*(diff(v(x, y, t), x))+K

pde2_lin; pde2_nonlin

diff(u(x, y, t), x)-(diff(v(x, y, t), y))-y*(diff(u(x, y, t), y))

-t

Download splitsen_in_lineair_en_niet_linear_werkend_def_12-8-2025_ac.mw

@salim-barzani 

Your approach with,
    eval(term, u(x, y, t) = T*u(x, y, t))/T
or,
    eval(term, {u(x, y, t) = T*u(x, y, t), v(x, y, t) = T*v(x, y, t)})/T
(and variants on that idea) seems to be poor. It doesn't work properly for some very simple examples.

note: I made an edit to fix a typo, where I accidentally used Z1 instead of Z2 in a step of dealing with pde2. I've corrected that now, in the code above. It happened to not affect your original example.

You could alternatively make a reusable procedure,
[edited]

restart

with(PDEtools, declare); with(PDEtools, undeclare)

undeclare(prime, quiet)

declare(u(x, y, t), quiet); declare(v(x, y, t), quiet)

 

Split := proc (expr::algebraic, F::(list(name)), L::(list(name))) local Z, dum, res; Z := map(proc (v) options operator, arrow; v = freeze(v) end proc, [indets(expr, And(specfunc(diff), satisfies(proc (fn) options operator, arrow; depends(fn, L) end proc)))[], `~`[apply](F, L[])[]]); res := thaw([selectremove(type, eval(expr, Z)+dum, Or(linear(`~`[rhs](Z)), satisfies(proc (f) options operator, arrow; not depends(f, `~`[rhs](Z)) end proc)))])[]; res[1]-dum, res[2] end proc

 

pde1 := diff(u(x, y, t), t)+diff(u(x, y, t), `$`(x, 3))+diff(u(x, y, t)*v(x, y, t), x)

diff(u(x, y, t), t)+diff(diff(diff(u(x, y, t), x), x), x)+(diff(u(x, y, t), x))*v(x, y, t)+u(x, y, t)*(diff(v(x, y, t), x))

pde2 := diff(u(x, y, t), x)-(diff(v(x, y, t), y))

diff(u(x, y, t), x)-(diff(v(x, y, t), y))

Split(pde1, [u, v], [x, y, t])

diff(diff(diff(u(x, y, t), x), x), x)+diff(u(x, y, t), t), (diff(u(x, y, t), x))*v(x, y, t)+u(x, y, t)*(diff(v(x, y, t), x))

Split(pde2, [u, v], [x, y, t])

diff(u(x, y, t), x)-(diff(v(x, y, t), y)), 0

Split(pde1, [u, v], [x, y, t]), Split(pde2, [u, v], [x, y, t])

diff(diff(diff(u(x, y, t), x), x), x)+diff(u(x, y, t), t), (diff(u(x, y, t), x))*v(x, y, t)+u(x, y, t)*(diff(v(x, y, t), x)), diff(u(x, y, t), x)-(diff(v(x, y, t), y)), 0

Download linear_ac2.mw

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