acer

32313 Reputation

29 Badges

19 years, 314 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Scot Gould I suspect that you're talking about a problem that affects all array-plots done using plots:-display, and not about the current query about how to deal with an odd number of plots.

There is a long-standing quirk in the GUI, such that when plots:-display is used for array-plotting the plots are all rendered left-aligned in their respective Table cells. There are two manual ways to correct for that: 1) using the mouse-pointer, select an outer (say, right) exterior border of the Table and adjust its position slightly, and 2) right-click in the Table and adjust its width property slightly. Those both cause the plots to properly realign at the Cell's centers.

This is one of several reasons why I prefer DocumentTools:-Tabulate over plots:-display, for handling Arrays/lists/Matrices/Vectors of plots, ie. the former doesn't suffer from this bug. And customized solutions using DocumentTools primitives even offers programmatic control of several choices of alignment.

I reported this problem (with intial Table rendering after calling plots:-display) years ago.

@Scot Gould Thank you. It was a typo in the accompanying text (now corrected). The code was already correct.

At the default settings (Digits=10<=15 and UsehardwareFloats=deduced) the Eigenvalues command will use parallelism within LAPACK.

The timings you describe hint that you are using a higher working precision, with software (rather than hardware double precision floats). We cannot tell because you haven't bothered to provide the code of your computation.

If the two attachments are the dotm (.m) format of your Matrices then you should inform us of that.

If this is a continuation of your very recent Question then you should provide the code you're using. You should provide the code anyway. It's OK if that means using read on .m files. But you really ought to provide the calls you make to do the computation as well as the Digits, etc.

There seems to be a problem with getting to the particular LAPACK function that benefits most from symmetric positive-definite A and symmetric B for the generalized eigen-problem. Instead it is calling the general (nonsymmetric routine) which is slower. But that is secondary to the slowdown you'll see moving from hardware Digits<=15 to higher precision software floats. But if you really need the higher precision then this might well be important. (I have not looked at the conditioning...yet.)

It's difficult to say more without your providing adequate details.

The fact that you appear to have unnecessarily split this into a separate Question thread is not very helpful.

note: The "generic" CLAPACK external code for higher (software float) precision is not parallelized. So it's slower for two reasons: 1) it's using software Maple-floats, and 2) it's not got the the same LAPACK parallelism as the hardware double-precision Intel MKL implementation. I don't know whether it is possible to beat the software float, generic, unparallelized eigenvector QRQZ compiled code with a Maple Threads/Grid parallelized implementation that does, say, separate nullspace (eigenspace basis) computations in parallel after finding the eigenvalues.

@mmcdara You may also declare names local at the top-level.

restart;

local Temperature, gamma:

Temperature[1] := 10;

10

gamma[r] := 10;

10

 

Then it becomes as straightforward to encapsulate such code within procedures. (Note that :- colon-minus syntax allows the original global name to still be referenced, if needed.)

IMO being able to declare the names (including a collection of nonpredetermined number) as local within a procedure is important. So often it becomes desirable to turn one-off code into reusable procedure(s). I reserve using double-underscore for base names that I'll never want to use to create several names (especially if the number may vary).

@mmcdara It's a personal preference, but I'd prefer to concatenate as little as possible. (The aliasing below re-uses the assigned displacements.)

Also, in general there might be other aliases in play, so you could target the ones you want here more precisely, again re-using displacements in the lookup.

restart:

mass_names := [A, B, C]

[A, B, C]

masses        := seq(M__||m, m in mass_names);
displacements := seq(x__||m, m in mass_names);

stiffnesses   := Matrix(3$2, (i,j) -> `if`(i=j, 0, K__||(cat(mass_names[i],mass_names[j]))));
dampings      := seq(xi__||m, m in mass_names);

masses := `#msub(mi("M"),mi("A"))`, `#msub(mi("M"),mi("B"))`, `#msub(mi("M"),mi("C"))`

displacements := `#msub(mi("x"),mi("A"))`, `#msub(mi("x"),mi("B"))`, `#msub(mi("x"),mi("C"))`

stiffnesses := Matrix(3, 3, {(1, 1) = 0, (1, 2) = `#msub(mi("K"),mi("AB"))`, (1, 3) = `#msub(mi("K"),mi("AC"))`, (2, 1) = `#msub(mi("K"),mi("BA"))`, (2, 2) = 0, (2, 3) = `#msub(mi("K"),mi("BC"))`, (3, 1) = `#msub(mi("K"),mi("CA"))`, (3, 2) = `#msub(mi("K"),mi("CB"))`, (3, 3) = 0})

xi__A, xi__B, xi__C

map(f->alias(f=f(t)), [displacements]):

a:=map(f->alias:-GlobalToContent[f],[displacements]):

velocities := diff(a, t):

accelerations := diff(a, t$2):

masses *~ accelerations =~ stiffnesses . `<,>`(displacements)
                           +~ dampings *~ velocities^~2;

Vector[column]([[M__A*(diff(x__A(t), `$`(t, 2))) = K__AB*x__B+K__AC*x__C+xi__A*(diff(x__A(t), t))^2], [M__B*(diff(x__B(t), `$`(t, 2))) = K__BA*x__A+K__BC*x__C+xi__B*(diff(x__B(t), t))^2], [M__C*(diff(x__C(t), `$`(t, 2))) = K__CA*x__A+K__CB*x__B+xi__C*(diff(x__C(t), t))^2]])

Download purpose_ac.mw

Alternatively, if you create your new list (ie, your `a`) prior to making the aliases then you don't have to resort to concoctions to form it afterwards. Eg,

restart:

mass_names := [A, B, C]:

masses        := seq(M__||m, m in mass_names):
displacements := seq(x__||m, m in mass_names):

stiffnesses   := Matrix(3$2, (i,j) -> `if`(i=j, 0, K__||(cat(mass_names[i],mass_names[j])))):
dampings      := seq(xi__||m, m in mass_names):

a := apply~([displacements],t):

alias~(Equate([displacements],a)):

velocities := diff(a, t):

accelerations := diff(a, t$2):

masses *~ accelerations =~ stiffnesses . `<,>`(displacements)
                           +~ dampings *~ velocities^~2;

Vector[column]([[M__A*(diff(x__A(t), `$`(t, 2))) = K__AB*x__B+K__AC*x__C+xi__A*(diff(x__A(t), t))^2], [M__B*(diff(x__B(t), `$`(t, 2))) = K__BA*x__A+K__BC*x__C+xi__B*(diff(x__B(t), t))^2], [M__C*(diff(x__C(t), `$`(t, 2))) = K__CA*x__A+K__CB*x__B+xi__C*(diff(x__C(t), t))^2]])

Download purpose_acc.mw

ps. You've been showing quite a bit of code using concatenations over the past year. But that doesn't allow safe local declaration within procedures, and so hinders programmatic encapsulation via procedures. Why do you use concatenated names instead of indexed names? (Both pretty-print with subscripts...) For example,

restart:

mass_names := [A, B, C]:

masses        := seq(M[m], m in mass_names):
displacements := seq(x[m], m in mass_names):

stiffnesses   := Matrix(3$2, (i,j) -> `if`(i=j, 0,
                                           K[mass_names[i],
                                             mass_names[j]])):
dampings      := seq(xi[m], m in mass_names):

a := apply~([displacements],t):

alias~(Equate([displacements],a)):

velocities := diff(a, t):

accelerations := diff(a, t$2):

masses *~ accelerations =~ stiffnesses . `<,>`(displacements)
                           +~ dampings *~ velocities^~2;

Vector[column]([[M[A]*(diff(x[A](t), `$`(t, 2))) = K[A, B]*x[B]+K[A, C]*x[C]+xi[A]*(diff(x[A](t), t))^2], [M[B]*(diff(x[B](t), `$`(t, 2))) = K[B, A]*x[A]+K[B, C]*x[C]+xi[B]*(diff(x[B](t), t))^2], [M[C]*(diff(x[C](t), `$`(t, 2))) = K[C, A]*x[A]+K[C, B]*x[B]+xi[C]*(diff(x[C](t), t))^2]])

Download purpose_accc.mw

@sand15 The addresses of the names returned by alias() are not the same as the addresses of the original global names. That is why your attempt did not work. (This relates to the magic way that the kernel builtin procedure alias:-BuiltinAlias works.)

An alternative might be,

restart;

alias(f=f(t)):
alias(g=g(t)):

a:=[indices(alias:-ContentToGlobal,':-nolist')];

            a := [f, g]

diff(a,t);

            d     d
           [-- f, -- g]
            dt    dt

I will file a report.

It returned both (as lists) in my Maple 16.02, but not in 17.02.

@John2020 

frontend(taylor,[cos(theta(t)),theta(t)=0,7],
         [{`+`,`*`,`=`,specfunc(cos)},{}]);

1-(1/2)*theta(t)^2+(1/24)*theta(t)^4-(1/720)*theta(t)^6+O(theta(t)^8)

subsindets(evalindets('taylor(cos(theta(t)),theta(t)=0,7)',
                      specfunc(theta),freeze),name,thaw);

1-(1/2)*theta(t)^2+(1/24)*theta(t)^4-(1/720)*theta(t)^6+O(theta(t)^8)

Download taylor_frontended.mw

And, as alternate, using frontend rather than explicitly freezing,

restart;
with(Physics):

r := x*(diff(theta(t), t))^2+y*(diff(varphi(t), t))^2:
g := (4*(f+T))*(diff(theta(t), t))^2+u*(diff(varphi(t), t))^2:

frontend(solve,[identity(r=g,diff(theta(t),t)),[x,y]],
         [{`+`,`*`,`=`,list,specfunc({identity})},{}]);

           [[x = 4 T + 4 f, y = u]]

@John2020 If you load the Physics package then the name diff is rebound.

I wasn't suggesting that you instead try to freeze calls to Physics:-diff. I was trying to suggest that you'd need to distinguish between Physics:-diff and the original, global name (and try to ensure that it is the latter which gets frozen). I posted my earlier reply from a phone, which is why it is terse. Sorry.

The code I originally provided freezes calls to the name diff, without forcibly referring to the original global name :-diff rather than the current binding of that name. But it will continue to work -- even with Physics loaded -- provided that you refer explicitly to that global name :-diff instead of the rebound diff.

restart;

with(Physics):

r := x*(diff(theta(t), t))^2+y*(diff(varphi(t), t))^2:
g := (4*(f+T))*(diff(theta(t), t))^2+u*(diff(varphi(t), t))^2:

solve(identity(subsindets(r=g,specfunc(:-diff),freeze),
               freeze(:-diff(theta(t),t))),
      [x,y]);

             [[x = 4 T + 4 f, y = u]]

@John2020 So why not provide the additional problematic example!?

It might be that you need to adjust what is frozen, eg. Physic:-diff versus diff or some such thing.

[edit. The above was not adequately clear, sorry. I meant that one might have to ensure that it is not Physics:-diff which is frozen, but the original global name :-diff. I meant that one might distinguish one versus the other.]

As the OP has learned an exact expression can be achieved here, though some of the methodology is not (as yet) explicitly given.

One way in which it can be achieved is through a change of variables under some reasonable assumptions that agree with the supplied plotting range.

This form of exact result does not suffer from the same numeric difficulties under fsolve. This particular form of the exact result can also be evaluated numerically without undesirable small (or float 0.0) imaginary components.

H := Int(1/sqrt(sin(x0)-sin(x)),x=0..x0);

Int(1/(sin(x0)-sin(x))^(1/2), x = 0 .. x0)

R := simplify(value(IntegrationTools:-Change(H,s=sin(x),s)))
     assuming s>0, s<1, x0>0, x0<Pi/2;

2^(1/2)*(EllipticK((1/2)*(2*sin(x0)+2)^(1/2))-EllipticF(1/(sin(x0)+1)^(1/2), (1/2)*2^(1/2)*(sin(x0)+1)^(1/2)))

f := unapply(R, x0):

f(0.1);

.448035371*2^(1/2)

g :=  alpha -> fsolve(x0 -> f(x0) - 2*sqrt(alpha), 0 .. Pi/2):

g(eval(R, x0=0.1));

.5614162073

plot(f, 0 .. Pi/2, labels = [x0, alpha])

Download inverse_function_with_fsolve_ac2.mw

[edit] In fact the same well-behaved exact form can be obtained even more directly. (The key here is to also pass the lower-bound assumption x0>0.)

int(1/sqrt(sin(x0)-sin(x)),x=0..x0) assuming x0>0, x0<Pi/2;

2^(1/2)*EllipticK((1/2)*(2*sin(x0)+2)^(1/2))-2^(1/2)*EllipticF(1/(sin(x0)+1)^(1/2), (1/2)*(2*sin(x0)+2)^(1/2))

 

@tomleslie The edit I suggested is not too hard to guess, and it causes the code to "run away" for me.

But I agree with you that a worksheet would be better. The OP's last Question was posted with a crucial detail missing (ie. Digits value, which made a significant difference).

@tomleslie I don't think that the OP intended

    if root_r != 0 then

as 1D input, since then it tests a factorial against zero (with unfortunate consequences on whether equations get augmented or not).

It's more likely that the OP intended,

   if root_r <> 0 then

It may even be that the OP used the != syntax in 2D Input in Maple itself, but made the mistake of thinking that literal could be sensibly shown here as plaintext code.

@Jaime_mc2 Why are you using,

   evalf(allvalues(RootOf(HermiteH(7, x), x)));

instead of, say,

   [fsolve(simplify(HermiteH(7, x)))];

 

Also, would a single solution suffice (if there were infinitely many)? I am wondering also why you are using solve rather than fsolve, when dealing with the equations.

First 115 116 117 118 119 120 121 Last Page 117 of 591