acer

32490 Reputation

29 Badges

20 years, 8 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Mac Dude Given the subtleties in play it is not prudent to comment on some changes you made to an example, without the complete code-to-reproduce. An upload would do. For example, did you change kernelopts(floatPi=false) earlier?

I suspect that the discrepancy in the floating-point roundoff effect is due to some interplay between the use of an Equation Label and kernelopts(floatPi).

Notice how the argument to sin differs in the last three statements below.

restart;

kernelopts(version);

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

(1)

kernelopts(floatPi=true): # default

Digits:=15;

15

(2)

trace(sin):

x1:=Int(3.52*10^8, ti = 0 .. 1):

R:=value(x1);

352000000.

(3)

sin(2*Pi*(3));

execute sin, args = 2211681228.12721
{--> enter sin, args = 2211681228.12721

<-- exit sin (now at top level) = -0.443987770092724e-5}

 

-0.443987770092724e-5

(4)

sin(2*Pi*R);

execute sin, args = 2211681228.12722

{--> enter sin, args = 2211681228.12722
<-- exit sin (now at top level) = 0.556012229902952e-5}

 

0.556012229902952e-5

(5)

sin('2*Pi'*R);

execute sin, args = 2211681228.12721

{--> enter sin, args = 2211681228.12721
<-- exit sin (now at top level) = -0.443987770092724e-5}

 

-0.443987770092724e-5

(6)

 

Download floatPiEqnLabel.mw

@Earl The frontend command is not necessary here, and I apologize if its inclusion was confusing. (I reach for it as a partial knee-jerk reaction when I see the derivatives, and I was trying some other things first.)

Here is the code, using eliminate without frontend. It is straightforward.

restart;

eq1 := -T*sin(theta(t)) = m*(diff(X(t), t, t)
       +L*(diff(theta(t), t, t))*cos(theta(t))
       -L*(diff(theta(t), t))^2*sin(theta(t))):

eq2 := T*cos(theta(t))-m*g = m*(diff(Y(t), t, t)
       +L*(diff(theta(t), t, t))*sin(theta(t))
       +L*(diff(theta(t), t))^2*cos(theta(t))):

K := eliminate({eq1,eq2},{m}):

conds := simplify(map(`=`,K[2],0)):

conds[1]/T;

L*(diff(diff(theta(t), t), t))+(diff(diff(X(t), t), t))*cos(theta(t))+sin(theta(t))*(diff(diff(Y(t), t), t)+g) = 0

(1)

 

Download eliminate_examp2.mw

There are a variety of ways to eliminate m from this particular pair of equations. Taking a difference of solve(...,m) calls is one alternative, but I find that somewhat ad hoc in the sense that it is not so generally applicable (generally, both equations might not contain or be solvable for m, etc).

To me this seems reasonable as a first approach: If you want to eliminate m then try the eliminate command.  Just my $0.02.

(I like Kitonum's ingenuity, hence my upvote for his answer.)

@vanzzy Post your followups on this topic/thread here. Most people cannot access the very old Maple 12 version. I won't have such access until later in the week at earliest.

I will delete more duplicates of this topic thread, even if you're just running into problems with the code people posted here.

@Mariusz Iwaniuk Calling simplify up front on the integrand is not a whole lot easier than calling convert, IMO.

But sure, it works for the second example above, even in Maple 2016 which is the version with which the OP saved his upload.

Using Maple 2016.2 for 64bit Linux the first example sometimes produces an error "too many levels of recursion", but more often loses the kernel connection, if simplify is first applied to the integrand as Mariusz has shown.

restart

kernelopts(version)

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

(1)

`assuming`([simplify(ln(r)*sin(k1*Pi*ln(r/Rs)/ln(r4/Rs))/r)], [Rs > 0, r4 > 0, r4 > Rs, k1::posint])

ln(r)*sin(k1*Pi*(ln(Rs)-ln(r))/(-ln(r4)+ln(Rs)))/r

(2)

`assuming`([int(ln(r)*sin(k1*Pi*(ln(Rs)-ln(r))/(-ln(r4)+ln(Rs)))/r, r = Rs .. r4)], [Rs > 0, r4 > 0, r4 > Rs, k1::posint])

Error, (in sprintf) too many levels of recursion

 

``

Download 2_integrals_ac.mw

@vanzzy The `size` option was introduced in Maple 18.

Your Maple 12 came ten major releases ago, and is ten years old.

It is unproductive to posts questions for Maple 12 without telling us of that aspect.  There is a drop menu for "Product" available when you submit a question, where you can specify the version.

The good news is that the `size` option is unrelated to your label issue. I used it only to make the plot show a little more clearly.

@vanzzy What version of Maple are you using?

In really old versions you can set _EnvExplicit to true , if the explicit option is not valid when calling solve.

Don't forget the complications about what the "first" root might be, or evaluate to in floating-point (numeric root-finding and convergence), or that the curves may cross, and that there are discontinuties.

restart;

gm := V -> 1/sqrt(1-V^2):
T := w-k*V:
S := w*V-k:

f := unapply(I*B*gm(V)^3*T^3+gm(V)^2*T^2
             -I*gm(V)^3*T*S^2-1/3*I*B*gm(V)^3*S^3*T-1/3*gm(V)^2*S^2,
             w,B,V,k):

Hgen := simplify(rationalize(f(w,B,V,k)));

(-V^2+1)^(1/2)*(((3*k^2-w^2)*V^2-4*V*k*w-k^2+3*w^2)*(-V^2+1)^(1/2)+I*(B*V^3*w^3+(-3*B*k*w^2-3*B*k^2+3*w^2)*V^2+3*k*w*(B*k+2*B-2)*V-B*k^3-3*B*w^2+3*k^2)*(V*k-w))/(3*V^4-6*V^2+3)

altcols := ["Violet","Orange","Navy"]:

Vlist :=  [ 0.8, 0.9, 0.99 ];

[.8, .9, .99]

_EnvExplicit := true;

true

HSols:=[solve(numer(Hgen),w)]:
indets(%,name);

{B, V, k}

#
#
#   HSols[1] is the "first" root
#
#

k := 1:
R1 := simplify(HSols[1], size):
Pk1 := plots:-display(
  seq(plot([Im](eval(R1, V=Vlist[i])), B=0.0 .. 1.0,
           linestyle=[solid], gridlines=false,
           adaptive=false, numpoints=100,
           color=altcols[i], legend=[typeset(V=Vlist[i])]),
      i=1 .. nops(Vlist)),
  title=typeset('k'=k),
  gridlines=false
);
k := 'k':

k := 5:
R1 := simplify(HSols[1], size):
Pk5 := plots:-display(
  seq(plot([Im](eval(R1, V=Vlist[i])), B=0.0 .. 1.0,
           linestyle=[solid], gridlines=false,
           adaptive=false, numpoints=100,
           color=altcols[i], legend=[typeset(V=Vlist[i])]),
      i=1 .. nops(Vlist)),
  title=typeset('k'=k),
  gridlines=false
);
k := 'k':

 

Download rootplot_alt_EnvExplicit.mw

@vanzzy Why on earth did you not tell us what you actually want in the first place? Up until the last Comment you did not specifiy that you wanted only the imaginary part of one "root" (in fact you conveyed something else, IMO). And you did not bother to give the V values before.

You didn't even bother to edit the original Question and supply the equations in copy'able text (or upload in a worksheet). You said you would, but you didn't. Someone else had to do it for you.

restart;

gm := V -> 1/sqrt(1-V^2):
T := w-k*V:
S := w*V-k:

f := unapply(I*B*gm(V)^3*T^3+gm(V)^2*T^2
             -I*gm(V)^3*T*S^2-1/3*I*B*gm(V)^3*S^3*T-1/3*gm(V)^2*S^2,
             w,B,V,k):

Hgen := simplify(rationalize(f(w,B,V,k)));

(((3*k^2-w^2)*V^2-4*V*k*w-k^2+3*w^2)*(-V^2+1)^(1/2)+I*(V*k-w)*(B*V^3*w^3+(-3*B*k*w^2-3*B*k^2+3*w^2)*V^2+3*k*w*(B*k+2*B-2)*V-B*k^3-3*B*w^2+3*k^2))*(-V^2+1)^(1/2)/(3*V^4-6*V^2+3)

altcols := ["Violet","Orange","Navy"]:

Vlist :=  [ 0.8, 0.9, 0.99 ];

[.8, .9, .99]

HSols:=[solve(numer(Hgen),w,explicit)]:
indets(%,name);

{B, V, k}

#
#
#   HSols[1] is the "first" root
#
#

k := 1:
Pk1 := plots:-display(
  seq(plot([Im](eval(HSols[1], V=Vlist[i])), B=0.0 .. 1.0,
           linestyle=[solid], gridlines=false,
           adaptive=false, numpoints=100,
           color=altcols[i], legend=[typeset(V=Vlist[i])]),
      i=1 .. nops(Vlist)),
  title=typeset('k'=k),
  gridlines=false
);
k := 'k':

k := 5:
Pk5 := plots:-display(
  seq(plot([Im](eval(HSols[1], V=Vlist[i])), B=0.0 .. 1.0,
           linestyle=[solid], gridlines=false,
           adaptive=false, numpoints=100,
           color=altcols[i], legend=[typeset(V=Vlist[i])]),
      i=1 .. nops(Vlist)),
  title=typeset('k'=k),
  gridlines=false
);
k := 'k':

 

 

Download rootplot_alt.mw

@erik10 By default dsolve/numeric will return a procedure which takes a numeric value t and returns a list containing t, x(t), and y(t), all estimated for the given numeric value of the argument. You can call this procedurelist output, as it's a single procedure which returns the results in a list. For any given numeric t this single procedure will return all the dependent results computed at that numeric value. If you wanted to use this with, say, fsolve then you'd have to pass fsolve some wrapping procedure which took in numeric t and extracted just a single dependent value from the list. That is doable, but a little awkward.

The alternative suggested below is to specify output=listprocedure which returns a list of procedures, each of which can take a numeric argument and return the estimated scalar result (respectively, for x(t), or y(t), etc). You can extract these procedures from the list, individually, and work with them individually. This is easy to use, and can also bring flexibility, as shown in the Answers to this thead which utilize fsolve.

Both forms of output involve procedures which use the same kinds of mechanisms to numerically solve the DE. The basic mechanism is to construct polynomial interpolants (or a piecewise of such) with coefficients that guarantee the accuracy specified when dsolve was called (at computed/requested data points).

Sometimes (eg, a system with two dependent variables, say x(t) and y(t) and a coupled system) it can be more efficient to use the single procedurelist output, rather than extract two procedures from listprocedure output. It can be beneficial to have the mechanism not duplicate results. But verifying such claims generally, via timing all that, is tricky (as garbage collection and memoization, etc, can come into play).

Now, let's talk about plotting. The plots:-odeplot call you made in a prior Comment computes 201 data points for each of x(t) and y(t). The plot driver will use use piecewise low order interpolation of those data points for the purpose of rendering the curve. It does not make sense to expect that the zoomed plot will be relatively accurate to more than a few decimal places for independent t values lying between those 201 computed data points. When you zoom in on a plot all that happens is that the plot driver recomputes from the previously computed data points; it doesn't request additional data points or additional accuracy. Zooming into a plot just uses the lower order interpolation of the plotter, which is completely separate from any accurate, high order interpolation within the dsolve/numeric engine or its returned procedures.

Now let's talk briefly about using dsolve/events versus the fsolve with the procedures extracted from output=listprocedure. The former has the potential to be quite a bit more efficient, as it allows the same mechanism that computes the dependent results to track itself and, say, detect a zero. It also knows how accurately it is computing the results. In contrast, the fsolve approach can resort to a numeric root-finding method like Newton's or secant, which can bounce all over the domain. And it has to find a suitable initial point that converges. And it has less idea about how accurate the function evaluations are (or could be requested) than does dsolve itself. So here we have a classic dichotomy of performance versus ease-of-use.

By coincidence, there was a recent, related question on stackexchange.com .

@gaurav_rs The numeric quadrature of the summation can be done pretty quickly at Digits=16 too, or at Digits=35. And using the closed form for the sum then for Digits=16 is faster still.

So the slowness you claim may be much more related to your original approach than it is to Digits.

restart;

Digits := 16:

f := (r,t,n) -> add(BesselJ(0, t*r)*r^(2*(i-1)), i=1..n):

CodeTools:-Usage(
  evalf(Int('f'(r,151.5793716314014,100), r=0..1))
);

memory used=172.66MiB, alloc change=38.01MiB, cpu time=1.05s, real time=1.05s, gc time=85.43ms

0.2274316631710838e-1

restart;

Digits := 16:

f := (r,t,n) -> Sum(BesselJ(0, t*r)*r^(2*(i-1)), i=1..n):

CodeTools:-Usage(
  evalf(Int('f'(r,151.5793716314014,100), r=0..1))
);

memory used=156.61MiB, alloc change=126.02MiB, cpu time=900.00ms, real time=871.00ms, gc time=76.16ms

0.2274316631710838e-1

restart;

Digits := 16:

f := unapply(sum(BesselJ(0, t*r)*r^(2*(i-1)), i=1..n), [r,t,n]):

CodeTools:-Usage(
  evalf(Int('f'(r,151.5793716314014,100), r=0..1, method=_Dexp))
);

memory used=52.02MiB, alloc change=4.00MiB, cpu time=354.00ms, real time=355.00ms, gc time=41.19ms

0.2274316631710838e-1

restart;

Digits := 35:

f := (r,t,n) -> add(BesselJ(0, t*r)*r^(2*(i-1)), i=1..n):

CodeTools:-Usage(
  evalf(Int('f'(r,151.5793716314014,100), r=0..1))
);

memory used=227.46MiB, alloc change=38.01MiB, cpu time=1.28s, real time=1.28s, gc time=103.39ms

0.22743166317108384811188537753590361e-1

restart;

Digits := 35:

f := (r,t,n) -> Sum(BesselJ(0, t*r)*r^(2*(i-1)), i=1..n):

CodeTools:-Usage(
  evalf(Int('f'(r,151.5793716314014,100), r=0..1))
);

memory used=211.26MiB, alloc change=126.02MiB, cpu time=1.20s, real time=1.12s, gc time=158.08ms

0.22743166317108384811188537753590361e-1

restart;

Digits := 35:

f := unapply(sum(BesselJ(0, t*r)*r^(2*(i-1)), i=1..n), [r,t,n]):

CodeTools:-Usage(
  evalf(Int('f'(r,151.5793716314014,100), r=0..1))
);

memory used=143.57MiB, alloc change=44.01MiB, cpu time=934.00ms, real time=892.00ms, gc time=118.77ms

0.22743166317108384811188537753590361e-1

 

Download evalfInt_example_16_35.mw

Perhaps only grainy images of solutions or suggestions could be offered in response.

@bliengme I think that it's still unclear what you want, and what is wrong with the following.

Do you accept what Maple's plots:-semilogplot gives you, for the case of a semilog plot with logarithmic scaling on the horizontal axis and linear scaling on the vertical axis? If so, then what's wrong with the reverse of that, as in the plots:-logplot behavior?

(Naturally, the labels could be anything you want.)

Or, perhaps it is the spacing of the subticks (the minor gridlines) that you want different?

Please state exactly what you want, not what some others think of definitions.

restart;

plots:-logplot(10^y, y=0..1,
               gridlines=true, labels=["y","x"]);

plots:-semilogplot(log10(x), x=1..10,
                   gridlines=true, labels=["x","y"]);

 

Download logplot.mw

 

@bliengme In a semilogarithmic plot, one axis has a logarithmic scale and the other axis has a linear scale.

From the description of the logplot command, "semi-logarithmic plot of functions where the vertical axis has a logarithmic scale." And the example shows that the horizontal axis has the usual linear scale.

So, how is the logplot command not supplying what you're asking for?

Anyway, you can also control the mode of either axis, separately. As described in the other link I gave.

Using dualaxisplot for this seems silly.

 

@eslamelidy You try a modestly high-order series solution, for various values of eq(z).

eslam_acc_complex.mw

@xinwmath What special form of generally invalid result from combine(..,symbolic) or simplify(...,symbolic) are you looking for?

restart;

expr3 := 3*(8*L__a+4*lambda+(4*I)*sqrt(3)*lambda)^(2/3)
         *(8*L__a-(4*I)*sqrt(3)*lambda+4*lambda)^(2/3)*lambda^2*(L__a-lambda);

3*(8*L__a+4*lambda+(4*I)*3^(1/2)*lambda)^(2/3)*(8*L__a-(4*I)*3^(1/2)*lambda+4*lambda)^(2/3)*lambda^2*(L__a-lambda)

ans3 := simplify(expand(combine(expr3,symbolic)),symbolic);

48*(L__a^2+L__a*lambda+lambda^2)^(2/3)*lambda^2*(L__a-lambda)

simplify(expand(combine(ans3 - expr3,symbolic)));

48*lambda^2*((L__a^2+L__a*lambda+lambda^2)^(2/3)-((L__a^2+L__a*lambda+lambda^2)^2)^(1/3))*(L__a-lambda)

simplify(expand(combine(ans3 - expr3,symbolic)),symbolic);

0

 

Download symbolic_sick.mw

First 235 236 237 238 239 240 241 Last Page 237 of 594