acer

32303 Reputation

29 Badges

19 years, 309 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

I am unsure of the root cause of the discrepancy. You could resolve that particular inconsistency with,

   int(..., numeric, method=_Dexp)

or,

   evalf(Int(..., method=_Dexp))

Also, your integrands may be complex-valued due to floating-point (ie. approximate) coefficients along with the presence of I, the imaginary unit.

You may get rid of small floating-point imaginary artefacts in results by wrapping the results in call to fnormal.

You may get rid of 0.0*I imaginary components of a result by wrapping in a call like simplify(...,':-zero').

Eg,
   simplify( fnormal( evalf(Int(..., method=_Dexp) ), ':-zero');

I don't understand why you apply evalf to the integrands passed to the quadrature call.

Btw, while you might be deliberately trying numeric quadrature for speed (say), much of your results can be obtained using exact symbolic computations (including erf calls in results) if you use exact rationals instead of floating-point values for your parameters, etc. In at least one way that can also symbolically show that the imaginary part of your integral is exactly zero.

The following is not real-valued for purely real x. But it's not awful.

restart;

kernelopts(version);

`Maple 2022.1, X86 64 LINUX, May 26 2022, Build ID 1619613`

with(IntegrationTools):

 

f := sin(x)/sqrt(2+sin(2*x));

sin(x)/(2+sin(2*x))^(1/2)

raw := eval(value(Change(expand(Int(f,x)), t=tan(x), t)),
            t=tan(x)):

Q := combine(simplify(convert(combine(simplify(raw,tan),
                                      radical,symbolic),
                              sincos),
                      symbolic));

-(1/2)*arctanh((2+sin(2*x))^(1/2)/(sin(x)+cos(x)))-(1/2)*arctan((-sin(x)+cos(x))/(2+sin(2*x))^(1/2))

plot([Re,Im](Q), x=-2*Pi..2*Pi);

Download intfun1.mw

I don't know how to get it as H or F in the Comments on the Question, which are real-valued for real x.

I'm not sure exactly how the names you used figure in the process.

But perhaps something like the following can be turned to your purpose.

restart;
randomize():

 

f := x -> x^2;

proc (x) options operator, arrow; x^2 end proc

ff := proc(x) option remember; f(x); end proc:

 

F := proc(f) local r;
  r := rand(1..100);
  f(r());
  f(r());
  NULL;
end proc:

 

F(ff);

 

T := op(4,eval(ff));

table( [( 46 ) = 2116, ( 62 ) = 3844 ] )

F(ff);

 

T := op(4,eval(ff));

table( [( 46 ) = 2116, ( 41 ) = 1681, ( 42 ) = 1764, ( 62 ) = 3844 ] )

ff := subsop(4=NULL, eval(ff)): # clear ff's remember-table

F(ff);

 

T := op(4,eval(ff));

table( [( 36 ) = 1296, ( 98 ) = 9604 ] )

Download memoi.mw

Here is one way:

restart;

kernelopts(version);

`Maple 2016.2, X86 64 LINUX, Jan 13 2017, Build ID 1194701`

interface(typesetting=extended):

A := RealRange(-2,Open(5));

RealRange(-2, Open(5))

B := RealRange(Open(0),7);

RealRange(Open(0), 7)

solve(convert(Or(x::A,x::B),relation));

RealRange(-2, 7)

solve(convert(And(x::A,x::B),relation));

RealRange(Open(0), Open(5))

solve(convert(And(x::A,Not(x::B)),relation));

RealRange(-2, 0)

Download some_set_stuff.mw

If you really want you could devise an overloaded scheme where those names union, minus, and intersection act similarly (say, by performing the operations such as I showed, underneath).

You may pass the option explicit (or explicit=true) to the solve command and obtain the explicit roots of a 4th degree polynomial directly, with a single command.

For your example, say,

    solve(Opt_effort_FOC = 0, e[n], explicit)

Inside the procedure p2 you need to use the global name :-debug as the name of the keyword option when it calls procedure p1.

p1 := proc(a, b, {debug::truefalse := false})
   print("p1", debug);
   a + b;
end proc:

p2 := proc(a, b, {debug::truefalse := false})
   print("p2", debug);
   p1(a, b, ':-debug' = true);
end proc:

p2(1, 2, 'debug' = true);

              "p2", true
              "p1", true
                  3

Otherwise your p2 is calling p1 with true=true, since inside p2
    debug = true
evaluates using the value of the p2 parameter that is also named debug.

Here is what's happening in your original. Notice the arguments that p2 passes to p1 include
   true = true
and not
   debug=true
which was what you'd expected.

restart;

p1 := proc(a, b, {debug::truefalse := false})
   print("p1", debug);
   a + b;
end proc:

p2 := proc(a, b, {debug::truefalse := false})
   print("p2", debug);
   p1(a, b, debug = true);
end proc:

trace(p1):
trace(p2):

p2(1, 2, 'debug' = true);
{--> enter p2, args = 1, 2, debug = true
                           "p2", true

{--> enter p1, args = 1, 2, true = true
                          "p1", false

                               3

<-- exit p1 (now in p2) = 3}
                               3

<-- exit p2 (now at top level) = 3}
                               3

Now observe what happens in the corrected version. Notice what arguments p2 passed to p1.

restart;

p1 := proc(a, b, {debug::truefalse := false})
   print("p1", debug);
   a + b;
end proc:

p2 := proc(a, b, {debug::truefalse := false})
   print("p2", debug);
   p1(a, b, ':-debug' = true);
end proc:

trace(p1):
trace(p2):

p2(1, 2, 'debug' = true);
{--> enter p2, args = 1, 2, debug = true
                           "p2", true

{--> enter p1, args = 1, 2, debug = true
                           "p1", true

                               3

<-- exit p1 (now in p2) = 3}
                               3

<-- exit p2 (now at top level) = 3}
                               3

By the way, I wrote as ':-debug' = true so that it would work even if the global name :-debug had also been assigned a value (at the top-level, say). Strictly speaking I didn't have to do that, since :-debug is a protected name. But for the same situation involving other (unprotected) names you could have that situation. So I wanted to illustrate to you the more robust way.

[note: See here for details I gave in your earlier Question thread, for plotting the 0-contour using interpolation of 0-height data that is computed more robustly. That seems to make this Question thread less relevant...]

You appear to be using Maple 2021, so you can use the newer Interpolation package (new in Maple 2018), which makes the task of interpolation "on-the-fly" easier.

The black-box that it returns can be subsequently polled at arbitrary x-y pairs, so you can easily use it with the contourplot command.

The problem is that you original data is quite flat for higher y-values, and the data you generated for it seems to lack enough precision to allow the 0-contour to be traced out nicely. (I alluded to this in your earlier Question thread, since I had already tried what I show below...)

In the plot below you can see that it has difficulty following the 0-contour, and slight inaccuracy in the original data can cause it to veer off the (exact) course.

Using ArrayInterpolation would not fix this issue, in and of itself. That command merely provides a mechanism to compute all the interpolated x-y target points in batch. I suspect that a good fix would require greater accuracy in the originally generated data -- which would be slower, and you earlier asked for it to be faster. It's possible that a scaling of the data to give it greater slope might be a kludge. But since the curve represents a function (in the mathematical sense) then it be more robustly found using fsolve (to find the alpha corresponding to each beta); again, see here.

restart;

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

with(plots):

excelfile:= FileTools:-JoinPath([kernelopts(':-homedir'),
                                "mapleprimes","N_data.xlsx"]):

NN := ImportMatrix(excelfile,'source'=':-Excel'):

(m,n) := op(1,NN);

51, 51

alpha := Vector([seq(0..10.0,'numelems'=n)],'datatype'=':-float[8]'):

beta := Vector([seq(0..10.0,'numelems'=m)],'datatype'=':-float[8]'):

 

func := Interpolation:-CubicInterpolation([beta,alpha],NN):

 

contourplot(func(x,y), x=0..10, y=0..10, 'contours'=[-1e-8,0,1e-8],
            'grid'=[m*2,n*2], coloring=[red,blue],
            'view'=[0..10,0..10], 'size'=[400,400]);

Download test1_interp_acc.mw

It's not very helpful to split these detail and highly related queries across multiple Question threads. At this time the .xlsx data file is at attachment in your earlier thread, but not in this one. And the quantitative nature of that data bears on which are decent courses of action.

Please try to keep your related queries in the same thread, even if nobody yet answers all to your satisfaction.

There is a special substructure of a 2D plot which the GUI recognizes as instructions for filling a popup box when the mouse-pointer's focus lies over that particular plot feature.

You can build that into your individual plot features (curves, points,...) by using the annotation option.

See the Help page for Topic plot,annotation .

The names generated by solve come in only a few flavours, so you might be able to get by without need for consideration of the names in the original expression f. For example,

restart;
assume(M::nonnegative):
f := 10*cos((-1+t)/sqrt(1+M))-10*cos(t/sqrt(1+M)):
g := solve({diff(f, t), t>0 }, t, allsolutions)[1][1]:

indets(rhs(g),And(name,Not(constant),
                  ':-suffixed(:-_)'));

           {_Z2}

indets(rhs(g),And(name,Not(constant),
                  ':-suffixed({:-_B,:-N,:-_Z})'));

           {_Z2}

I don't see a natural reason to attempt a heavy handed approach like as follows. (The original expression might have already had such generated names, in which case a combined approach might be needed.)

convert(indets(f,And(name,Not(constant))),`global`):
convert(indets(rhs(g),And(name,Not(constant))),`global`) minus %;

           {_Z2}

Simpler results can be had in Maple 2022.1, using the method below.

Perhaps someone can come up with something more direct. I programmatically determine the second root, to be less ad hoc about passing a range to Calculus1:-Roots.

The Calculus1:-Roots and maximize commands have a special (peppering) technique for determining values of the _Zxx assumed-integer parameter that solve returns here when passed its allsolutions option.

The solve command has improved in the past few releases, as far as returning a finite number of explicit roots when passed inequalities that denote a finite range. But it still seems not generally as strong as that internal technique of Roots and maximize. For this example solve returns a general formula when supplied its allsolutions option, but seems to return NULL when instead passed various reasonable inequalities, e.g t>0,t<10.

Here it is first with Maple 2015.2,

restart;

kernelopts(version);

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

with(plots):

shock := piecewise(t <0, 0, t < 1, 10, 0):
sys   := {(M__p+M__a)*diff(x(t), t$2)=M__p*shock-x(t), x(0)=0, D(x)(0)=0}

sys := {(`#msub(mi("M"),mi("p"))`+`#msub(mi("M"),mi("a"))`)*(diff(x(t), t, t)) = `#msub(mi("M"),mi("p"))`*piecewise(t < 0, 0, t < 1, 10, 0)-x(t), x(0) = 0, (D(x))(0) = 0}

sol := unapply(rhs(dsolve(sys)), (M__p,M__a))

proc (M__p, M__a) options operator, arrow; piecewise(t < 0, 0, t < 1, -10*M__p*cos(t/(M__p+M__a)^(1/2))+10*M__p, 1 <= t, 10*M__p*cos((-1+t)/(M__p+M__a)^(1/2))-10*M__p*cos(t/(M__p+M__a)^(1/2))) end proc

# Computation of tend for M__p=M__a=2
#
# Im expecting to get a value around 6.8

gensol := solve({sol(2,2),t>0},t,allsolutions);

{t = -2*arctan((cos(1/2)-1)/sin(1/2))+2*Pi*_Z2}

RootFinding:-NextZero(unapply(sol(2,2),t),0);
secondroot := RootFinding:-NextZero(unapply(sol(2,2),t),%);

6.783185307

13.06637061

# Maple 2022.1 gives 2*Pi+1/2 directly here, without using
# 'expand' with 'Roots', and without even any simplification
# done to the result
#
cands := {Student:-Calculus1:-Roots(expand(sol(2,2)),t=0..secondroot)[]};
tend := [op(cands minus {0})][1];
evalf(tend);

{0, -2*arctan((cos(1/2)-1)/sin(1/2))+2*Pi}

-2*arctan((cos(1/2)-1)/sin(1/2))+2*Pi

6.783185308

# Maple 2022.1 gives 0 here, without calling 'expand'.
#
xtend := simplify(expand(eval(sol(2,2),t=tend)));

0

# Maple 2022.1 gives 40*sin(1/4) directly here, without even calling 'simplify'
#
xmax := maximize(20*cos(((t - 1)*sqrt(4))/4) - 20*cos(sqrt(4)*t/4));
xmax := simplify(xmax);
evalf(xmax);

40/(1+cos(1/2)^2/sin(1/2)^2+2*cos(1/2)/sin(1/2)^2+1/sin(1/2)^2)^(1/2)

20*2^(1/2)*(-cos(1/2)+1)^(1/2)

9.896158366

# Maple 2022.1 gives Pi+1/2 here, without even using 'simplify'.
#
tmax := Student:-Calculus1:-Roots(expand(sol(2,2))=xmax,t=0..10)[1];
tmax := simplify(tmax);
evalf(tmax);

-2*arctan((cos(1/2)+1)/sin(1/2))+2*Pi

-2*arctan((cos(1/2)+1)/sin(1/2))+2*Pi

3.641592654

display(
  pointplot([[tmax,xmax], [tend,xtend]],
            symbolsize=20, symbol=solidcircle),
  plot(sol(2,2),t=0..10, size=[400,200])
);

 

Download ToyProblem_ac2015.2.mw

And here that is in Maple 2022.1,

restart;

kernelopts(version);

`Maple 2022.1, X86 64 LINUX, May 26 2022, Build ID 1619613`

with(plots):

shock := piecewise(t <0, 0, t < 1, 10, 0):
sys   := {(M__p+M__a)*diff(x(t), t$2)=M__p*shock-x(t), x(0)=0, D(x)(0)=0}

sys := {(M__p+M__a)*(diff(x(t), t, t)) = M__p*piecewise(t < 0, 0, t < 1, 10, 0)-x(t), x(0) = 0, (D(x))(0) = 0}

sol := unapply(rhs(dsolve(sys)), (M__p,M__a))

sol := proc (M__p, M__a) options operator, arrow; piecewise(t < 0, 0, t < 1, -10*M__p*cos(t/sqrt(M__p+M__a))+10*M__p, 1 <= t, -10*M__p*cos(t/sqrt(M__p+M__a))+10*M__p*cos((-1+t)/sqrt(M__p+M__a))) end proc

# Computation of tend for M__p=M__a=2
#
# Im expecting to get a value around 6.8

gensol := solve({sol(2,2),t>0},t,allsolutions);

{t = -2*arctan((-1+cos(1/2))/sin(1/2))+2*Pi*_Z2}

RootFinding:-NextZero(unapply(sol(2,2),t),0);
secondroot := RootFinding:-NextZero(unapply(sol(2,2),t),%);

6.783185307

13.06637061

cands := {Student:-Calculus1:-Roots(sol(2,2),t=0..secondroot)[]};
tend := [op(cands minus {0})][1];
evalf(tend);

{0, 1/2+2*Pi}

1/2+2*Pi

6.783185308

xtend := simplify(eval(sol(2,2),t=tend));

0

xmax := maximize(20*cos(((t - 1)*sqrt(4))/4) - 20*cos(sqrt(4)*t/4));
evalf(xmax);

40*sin(1/4)

9.896158372

tmax := Student:-Calculus1:-Roots(sol(2,2)=xmax,t=0..10)[1];
evalf(tmax);

1/2+Pi

3.641592654

display(
  pointplot([[tmax,xmax], [tend,xtend]],
            symbolsize=20, symbol=solidcircle),
  plot(sol(2,2),t=0..10, size=[400,200])
);

 

Download ToyProblem_ac2022.1.mw

Is this the kind of thing you were trying to accomplish?

I wasn't sure whether you wanted only the floating-point approximate root where f(n)=g(n), or an exact result for that, or the next highest integer above that where f(n)>g(n).

restart;

f := n -> 8*n^2:

g := n -> 64*n*log[2](n):

fsolve(f(n)=g(n),n=5..infinity);

43.55926044

ceil(%);

44

plot([f(n),g(n)],n=5..50, size=[500,300]);

Download ex1.mw

An exact result is also possible, in terms of the LambertW special-function.

solve({f(n)=g(n),n>5},n):

    { n = -8/ln(2)*LambertW(-1,-1/8*ln(2)) }

evalf(%);

   {n = 43.55926044}

Your Question shows an unevaluated call to intervalsolve, as a returned value. There is no such command in stock Maple. There is such a command in the Gym package (used by Danish schools). If using that command was your intention then you made two mistakes:

1) Your unevaluated return value shows that you failed to refer to its name properly. You could either use the command's fully qualified name (by including its package name in the call), or you could first load the package and then use the command's short name. See below, for examples of both.

2) The calling sequences of arguments that you passed is also incorrect for that Gym package's command. See below.

Here below is that command working in Maple 2021.2, called in a few alternate ways.

(These produce the same results in Maple 2022.1 and Maple 2020.2. You didn't tag your Question with your Maple version number; that can be useful information for us if it's an old version.)

restart;

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

f := n -> 8*n^2:

g := n -> 64*n*log[2](n):

Gym:-intervalsolve(f(n)=g(n),n=5..infinity);

[-8*LambertW(-1, -(1/8)*ln(2))/ln(2)]

Gym:-intervalsolve(f(n)=g(n),n=5..50);

[-8*LambertW(-1, -(1/8)*ln(2))/ln(2)]

with(Gym):

intervalsolve(f(n)=g(n),n=5..infinity);

[-8*LambertW(-1, -(1/8)*ln(2))/ln(2)]

intervalsolve(f(n)=g(n),n=5..50);

[-8*LambertW(-1, -(1/8)*ln(2))/ln(2)]

evalf(%);

[43.55926044]

Download intervalsolve_Gym.mw

plot( x, labels = [ `&Aring;`, y ],
         labelfont = [ Times, 14 ] );

Alternatives include the following (the Mapleprimes backend doesn't render them all, but my Maple 2022.1 GUI does):

`&Aring;`
`&#8491;`
`&#x212B;`

There is also,
   `#mo("&Aring;");`
which renders in upright Roman, instead of italic. It's possible to get alternate fonts with this upright variant, but it's messier and in that case you might well find Tom's Unit(angstrom) suggestion easier. For the italic, the above seems to work as well with alternate font families.

See also here, if you're more generally interested in fantastic beasts and where to find them.

unapply(x*cos(x), x)

or, in your example,

unapply(f,x)

where x is the intended symbol that appears in the expression to which f evaluates.

Below, CL is a list of the coefficients of powers of e^(-gamma*tau).

Is that what you're after?

restart

f := f__0(eta)+e^(-gamma*tau)*F(eta, tau)

f__0(eta)+e^(-gamma*tau)*F(eta, tau)

(1)

pde := diff(f, eta, eta, eta)+6*k*(diff(f, eta, eta, eta))*(diff(f, eta, eta))^2+(2/3)*f*(diff(f, eta, eta))-(1/3)*(diff(f, eta))^2-(diff(f, eta, tau))+(2/3)*tau*((diff(f, eta))*(diff(f, eta, tau))-(diff(f, eta, eta))*(diff(f, tau)))+1/3 = 0

diff(diff(diff(f__0(eta), eta), eta), eta)+e^(-gamma*tau)*(diff(diff(diff(F(eta, tau), eta), eta), eta))+6*k*(diff(diff(diff(f__0(eta), eta), eta), eta)+e^(-gamma*tau)*(diff(diff(diff(F(eta, tau), eta), eta), eta)))*(diff(diff(f__0(eta), eta), eta)+e^(-gamma*tau)*(diff(diff(F(eta, tau), eta), eta)))^2+(2/3)*(f__0(eta)+e^(-gamma*tau)*F(eta, tau))*(diff(diff(f__0(eta), eta), eta)+e^(-gamma*tau)*(diff(diff(F(eta, tau), eta), eta)))-(1/3)*(diff(f__0(eta), eta)+e^(-gamma*tau)*(diff(F(eta, tau), eta)))^2+e^(-gamma*tau)*gamma*ln(e)*(diff(F(eta, tau), eta))-e^(-gamma*tau)*(diff(diff(F(eta, tau), eta), tau))+(2/3)*tau*((diff(f__0(eta), eta)+e^(-gamma*tau)*(diff(F(eta, tau), eta)))*(-e^(-gamma*tau)*gamma*ln(e)*(diff(F(eta, tau), eta))+e^(-gamma*tau)*(diff(diff(F(eta, tau), eta), tau)))-(diff(diff(f__0(eta), eta), eta)+e^(-gamma*tau)*(diff(diff(F(eta, tau), eta), eta)))*(-e^(-gamma*tau)*gamma*ln(e)*F(eta, tau)+e^(-gamma*tau)*(diff(F(eta, tau), tau))))+1/3 = 0

(2)

CL := frontend(PolynomialTools:-CoefficientList, [lhs(pde), e^(-gamma*tau)])

[diff(diff(diff(f__0(eta), eta), eta), eta)+6*k*(diff(diff(diff(f__0(eta), eta), eta), eta))*(diff(diff(f__0(eta), eta), eta))^2+(2/3)*f__0(eta)*(diff(diff(f__0(eta), eta), eta))-(1/3)*(diff(f__0(eta), eta))^2+1/3, diff(diff(diff(F(eta, tau), eta), eta), eta)+12*k*(diff(diff(diff(f__0(eta), eta), eta), eta))*(diff(diff(f__0(eta), eta), eta))*(diff(diff(F(eta, tau), eta), eta))+6*k*(diff(diff(diff(F(eta, tau), eta), eta), eta))*(diff(diff(f__0(eta), eta), eta))^2+(2/3)*F(eta, tau)*(diff(diff(f__0(eta), eta), eta))+(2/3)*f__0(eta)*(diff(diff(F(eta, tau), eta), eta))-(2/3)*(diff(f__0(eta), eta))*(diff(F(eta, tau), eta))+gamma*ln(e)*(diff(F(eta, tau), eta))-(diff(diff(F(eta, tau), eta), tau))+(2/3)*tau*((diff(f__0(eta), eta))*(-gamma*ln(e)*(diff(F(eta, tau), eta))+diff(diff(F(eta, tau), eta), tau))-(diff(diff(f__0(eta), eta), eta))*(-F(eta, tau)*gamma*ln(e)+diff(F(eta, tau), tau))), 6*k*(diff(diff(diff(f__0(eta), eta), eta), eta))*(diff(diff(F(eta, tau), eta), eta))^2+12*k*(diff(diff(diff(F(eta, tau), eta), eta), eta))*(diff(diff(f__0(eta), eta), eta))*(diff(diff(F(eta, tau), eta), eta))+(2/3)*F(eta, tau)*(diff(diff(F(eta, tau), eta), eta))-(1/3)*(diff(F(eta, tau), eta))^2+(2/3)*tau*((diff(F(eta, tau), eta))*(-gamma*ln(e)*(diff(F(eta, tau), eta))+diff(diff(F(eta, tau), eta), tau))-(diff(diff(F(eta, tau), eta), eta))*(-F(eta, tau)*gamma*ln(e)+diff(F(eta, tau), tau))), 6*k*(diff(diff(diff(F(eta, tau), eta), eta), eta))*(diff(diff(F(eta, tau), eta), eta))^2]

(3)

simplify(add(CL[i+1]*(e^(-gamma*tau))^i, i = 0 .. 3)-lhs(pde))

0

(4)

simplify(sum(CL[i+1]*(e^(-gamma*tau))^i, i = 0 .. 3)-lhs(pde))

0

(5)

NULL

Download How_to_collect_Coefficent_ac.mw

ps. Your expression contain e^(-gamma*tau) which is not the same as exp(-gamma*tau). Is that intentional?

First 61 62 63 64 65 66 67 Last Page 63 of 336