Carl Love

Carl Love

28020 Reputation

25 Badges

12 years, 301 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Some of the below may duplicate things that acer already said. I didn't follow the links to worksheets that he gave (can't do it from my phone), so I haven't read the entirety of his Answer.

You wrote:

  •  fsolve and plot only work with the proc if a function was defined in their calling sequence as shown below.

That's totally a red herring reinforced by confirmation bias. (Not your fault: The human mind evolved to fall prey to confirmation bias.) The truth is that in the examples that you show that didn't work, the first argument was not a proc; and in the examples that you show that did work, it wouldn't have mattered whether the proc was defined inside or outside of the command call.

Here are your examples (in 1D input), 1st those that didn't work. Given: is exactly as you've defined above; it is a proc.

fsolve(f(x)[1] - x, w);
plot(f(x)[1] - x, 0. .. 0.02);

The problem is that f(x)[1] - x is not a proc. So the fact that f isproc is made irrelevant in these examples.

Now those that worked:

fsolve(x-> ValveLaminarMassFlowProc(x, P__1, P__2, rho__1, T__1, rho__eta)[1] - x, w);
and likewise for the plot.

You could just as well have defined:

f1:= x-> 
    ValveLaminarMassFlowProc(x, P__1, P__2, rho__1, T__1, rho__eta)[1] - x
:

and then done: 

fsolve(f1, w);
plot(f1, 0. .. 0.02);

and those would've worked also, which proves that having the procedure definition inside the command is irrelevant.

You also asked about trapping the error inside the proc. I suppose that that can be done, but I'd rather show you how to prevent the error inside the proc. This will involve two programming concepts that may be new to you: type-checking and returning unevaluated.

Given: The same and f1 as above, define:

f2:= proc(x::algebraic)
    if x::numeric then f1(x)[1] - x else 'procname'(x) fi
end proc:

The x::numeric is the "type check", and the 'procname'(x) being the last thing evaluated when the else branch is taken is the "returning unevaluated". The keyword procname is used (only inside procedures) to refer to the proc's own name (f2 in this case), i.e., it's a form of indirect self reference. The quotes around the procname are what makes it unevaluated.

That's all there is to it! It's not actually necessary to have a cascade of procedures ff1f2. I just did that to save myself some typing. All this can be done in the original f.

Now you can use f2 like this:

plot(f2(x), x= 0. .. 0.02);
and likewise for fsolve.

My own recent bad restart experience:

What you describe is a serious bug in restart. I also had a recent (about 3 weeks ago) experience that convinced me beyond any reasonable doubt that

  1. Some kernel memory had survived a restart,
  2. That memory was causing an error, and
  3. An expression that had been part of the computation before the restart was being blamed in the error message as the cause of the error even though the computation after the restart did not have the remotest connection to the one before.

When run in a new kernel, the second computation ran without error (as expected). I was able to duplicate this whole scenario multiple times. The expression was a totally ordinary algebraic expression that had been the denominator of an ODE IVP that I had solved numerically (via method= classical[rk4]). Indeed, the expression was d:= (x(t)^2 + y(t)^2)^(3/2). It's unique enough that there's no chance that it would appear anywhere at all in the second computation simply by coincidence. The second computation had nothing to do with differential equations or even algebra. The command that errored was (I think) type(e, with_unit(realcons, 'a', 'b')) where was something very simple like 100*Unit(ohm). (I'm sure about the unit, but not about the coefficient. The coefficient was definitely some positive real number.)

Minutiae regarding my enviroment (in case it may help someone investigating this):

  1. Neither computation involved any user-initiated disk I/O (e.g., savereadsavelib) whatsoever, nor were any directory-checking or directory-changing commands used (e.g., currentdir).
  2. I almost never load packages with with, and I certainly didn't in this situation.
  3. I use Standard interface, 1D Input, 2D output with extended typesetting.
  4. Windows 10; Maple 2020.1 (I don't trust .2).
  5. My initialization file is very basic. My libname includes a few standard 3rd-party packages such as DirectSearch, but none of those were used in either computation.

  6. I don't get the Physics updates because installing them seems like a hassle and I do no symbolic work in physics and very little in special functions or differential equations. (I mostly do numerical analysis, combinatorics, operations research, symbolic linear algebra, logic, and stuff with very complicated types.) However, this new Latex is intriguing, and I hope that it's preinstalled in the next major release.

  7. I used no embedded components in these two computations, not even a plot.

(I have no idea why the MaplePrimes editor double-spaced those last two points, which I've tried to correct several times.)


What restart is supposed to do:

First: (You probably know this already, because I've said it many times over the years (not necessarily to you).) You need to put restart in its own execution group. This is mentioned on its help page. However, it does usually---but not always---work correctly if you don't follow this advice.

AFAIK, restart is supposed to do all these things:

  1. Execute the ModuleUnload procedure of any module in use that has one. (Most modules do not, but I think that Physics has a very long one. You can see what it's doing by setting printlevel:= 10 (or greater) in a separate execution group immediately before doing the restart.)
  2. Reset all variables to their initial state.
  3. Collect the garbage.
  4. Execute your initialization file, if you have one.

And I'm unsure about whether it's supposed to do these:

  1. Close files opened by the user. Surely, it should close them (after posting their write buffers), because their pointers will disappear with the restart, thus making them zombies. But I don't see documentation that it does close them.
  2. Update the status bar to reflect the fact that all the garbage has been collected. (I'm sure that it doesn't do that updating; the part that I'm unsure about is whether it's supposed to. It seems shoddy that it doesn't.)
  3. Things that it's supposed to do, does do, or doesn't do to the GUI are of no concern to me at the moment (except for that status bar point).
     

Something that you can do as a workaround that's a lot faster than what you're doing now:

If you save and reopen your worksheet, that will start a brand-new kernel (a new invocation of mserver.exe). This takes much less time than opening and closing your entire Maple session.

And, as I've said several times before (not necessarily to you), this exact same advice can be used to correct loss-of-kernel errors (rather than closing the entire Maple session, as is incorrectly recommended by the loss-of-kernel error message).

Yes, it's simply indices(node).

P.S.: I believe that an earlier version of this Question referred to node as a "Matrix". I guess that you realized that a matrix can only have two indices and thus edited out that word. So, let me reassure everyone that the command indices will work exactly as shown on any table or rtable (ArrayMatrix, or Vector).

You are ignoring everything that I tried to teach you in my responses to your previous version of this Question, thus causing unnecessarily duplicated effort for other responders. Here's a summary:

  1. Separate the code that numerically evaluates the function from the code that numerically integrates it.
  2. Replace your conditional (if-based) definition of with an unconditional one using abs. You already agreed, after some cajoling, that this is correct, so why go back to the if?
  3. There is some severe numerical instability (or perhaps it's extreme oscillation) in this function that has thwarted every attempt that I've made (using up to 8 million evaluation points) to numerically integrate it (using many different algorithms) to more than 1 digit precision. Your attempts to use the trapezoid rule or Simpson's rule (here and in the prior Question) are worthless if you can't get an upper bound on the error.

Do you think that all the advice that I already gave is worthless just because it didn't finally produce the integral?

Do you think that adding up a bunch of nearly random numbers and calling it an "integral" is worthwhile? 

To proceed in a meaningful way, the numerical error needs to be controlled. I'm hoping that some other respondent can provide some insight into that. Once the error is controlled, one of Maple's standard numerical integration algorithms (which are far more sophisticated than Newton-Cotes (trapezoid, Simpson's, etc.)) should be able to handle it. That makes my point #1 above even more relevant.

The errors and/or warnings that you're getting come when you try to use indexing to add or modify an element in an existing list with more than a certain number of elements (I think the cut-off is 100 elements). In other words, any statement of the form L[k]:= ... will give this message if L is a list and has 100 or more elements. But note that this form of assignment statement is a bad thing to do even in the cases where Maple doesn't obstensibly object to it. It is also just as bad (for efficiency) to add elements to an existing list, set, or sequence without using indexing, although Maple (unfortunately) doesn't complain about it. Here are the 3 most common bad statements of this form:

L:= [op(L), x]; #adding to a list
S:= S, x; #adding to a sequence
S:= S union {x}; #adding to a set

Lists, sets, and sequences are three of the fundamental indexable container data structures of Maple. Those particular three are the immutable containers. (I'm omitting here any further discussion of containers with arithmetic properties such as explicit sums and products, although it's worth keeping in mind that these too are immutable.) All of these statements (including the explicitly indexed assignment) are bad because they don't truly modify the container; instead, they construct a new container from scratch by copying the old container plus the new element. Indeed, immutable literally means "can't be changed", and there's truly no way to change these containers. These statements are one of the main causes of inefficiency in Maple programs.

But Maple also has three fundamental indexable mutable containers: tables, rtables, and records. An rtable is a VectorArray, or Matrix. Tables can be explicitly created by the table command, or they're implicitly created by assignment statements of the form T[k]:= x where is a previously unused name (and and x can be anything, or even sequences of things, or even NULL). Records are created with the Record command; they're fairly specialized, so I won't discuss them further here.

One nice feature of all these mutable containers is that the elements in these containers can be anything at all, even lists, nested lists, sequences, or even NULL  (this underlined part is the direct answer to your specific question). So you can (and probably should) use a one-dimensional Array for the specific situation that you posed:

A:= Array(1..10, i-> [[0,0], 0]);

If you want to extend this Array to 11 elements, you may do this:

A(11):= ...; #note the ( ) instead of [ ]! 

or this is much easier (new to Maple 2018) because you don't need to know the last index:

A,= ...; #that's not a typo!

I recommend that you NOT use a two-dimensional Array for the situation that you described.

Using a table instead of an Array is comparably as efficient and avoids a lot of syntax:

T:= table(): #optional, although recommended.
for k to 10 do T[k]:= [[0,0], 0] od:
T[11]:= [[1,2], 3]: #must be [ ] after T!

There are no restrictions on table indices; unlike Array indices, they don't need to be consecutive integers, and indeed they don't need to be integers at all.

If you come to a point in your code where you no longer need to modify your table or Array (or add elements to it), it's perfectly fine and efficient to convert it to a list at that point. There are several ways. The easiest are these:

#Array A:
L:= [seq](A);

#table T:
L:= [entries](T, 'nolist'):
#or, if T contains "naked" sequences or NULLs this is better:
L:= [entries](T):

If you want your table elements to be listed in a certain order, you'll need to do something different, which isn't difficult, but it may be easier to use an Array in that case.

As of Maple 2019, it's possible to create lists, sets, and sequences of arbitrary non-predetermined length with do loops (forwhile, etc.) without using any of the "bad" statements that I discussed in the first two paragraphs and without any efficiency penalty. For example,

L:= [for k to 10 do [[k+1,k], k^2] od];

The outer square brackets make it a list. They could be changed to { } for a set or ( ) for a sequence. These loops have the full flexibility of any Maple code:

  • There's no limit to how much code is in the do loop.
  • It can be any code, including other loops and if statements. 
  • The length of the container need not be predetermined: The loop can be exited by a break statement or while or until condition.
  • These container-constructing loops can be embedded in other code anywhere the container itself can appear; they needn't be explcitly assigned to anything.

To my mind, these container-constructing do​​​​​​ loops are the greatest innovation in Maple since objects were introduced (in Maple 16, I believe). It's a shame that they're not more widely known, because they provide the easiest solution to the problems that I discussed in the first two paragraphs, which are very serious efficiency burdens.

The following procedure will be significantly faster (than even my modification of Kitonum's procedure) for numbers with many divisors:

SumOddDivisors:= proc(n::integer)
local k, q:= n;
    for k while q::even do q:= iquo(q,2) od;
    2^(k-1)*NumberTheory:-SumOfDivisors(q)
end proc
:  
#Example:
n:= 2^39*3*5*7*11*13*17*19*23*nextprime(2^29);
               n := 32922697295031156088441405440

CodeTools:-Usage(SumOddDivisors(n));
memory used=7.57KiB, alloc change=0 bytes, cpu time=0ns,
real time=0ns, gc time=0ns

                 82255314605128880924739502080

The reason that this is significantly faster is that the sum of the divisors of a number can be computed without computing the divisors themselves; only the factorization into prime powers is needed. Simple example: 81 = 3^4; the sum of its divisors is (3^5 - 1)/(3 - 1) from the well-known formula for a finite geometric sum.

Your equation can be easily converted to a polynomial, and this'll cause fsolve to use its dedicated polynomial algorithm, which often produces better results. However, the fundamental problem is still as acer says: You need to increase Digits because the root is so close to 1. But the polynomial technique may be beneficial to you for some other problems. 

Digits:= 16:
r2poly:= denom(rhs(r2))-numer(rhs(r2))/lhs(r2);
                              7                       2
         r2poly := (ksi + 100)  (ksi - 1) (2 ksi - 99) 

                                  5    5
            + 5.894442105097472 10  ksi 


fsolve(r2poly, ksi);
             -119.4197464595025, 0.9999999999994157


By default, all real roots are returned. The option complex could be added.

Note that xi will prettyprint as the Greek letter that you spell ksi​​​​​​.

So, why did you use option outliers= true?

I think this is what you want:

Statistics[BoxPlot](
    X, orientation= horizontal, deciles= false, mean= false,
    axis[2]= [tickmarks=0], view= [1..80, -.5..1.5]
);

There are two wedges that are bounded by those 3 surfaces: In one case, 0 <= z <= -y; in the other, -y <= z <= 0. Here's one way to plot them:

plot3d([1, theta, z], theta= -Pi..0, z= 0..-sin(theta), coords= cylindrical, scaling= constrained);
plot3d([1, theta, z], theta= 0..Pi, z= 0..-sin(theta), coords= cylindrical, scaling= constrained);

The only difference in those two commands is the range of theta.

If you examine the evaluated result of sqrt(x) with lprint, you'll see x^(1/2). Thus, it doesn't contain the literal expression sqrt. The command has is looking for the literal expression, whereas hastype and types in general allow for more subtlety.

P.S. In the title, I mean that sqrt(x) is not a function in the special sense that that word is used in Maple. Of course it's a function in the ordinary mathematical sense (as long as x is a free variable).

You can use solve to put the initial conditions into standard form:

ics2:= solve({ics}, {a0,a2,b0,b2,a1,b1}(0));
            /        3           3          3           15  
   ics2 := { a0(0) = --, a1(0) = -, a2(0) = --, b0(0) = --, 
            \        16          8          16          16  

             1          -1\ 
     b1(0) = -, b2(0) = -- }
             8          16/ 


dsolve({odes} union ics2, numeric);

It is NOT a solution. Perhaps you mean that z0(t) = t^2*exp(-t) is a solution? This can be verified with odetest:

z0:= t-> t^2*exp(-t):
ode:= diff(z(t),t,t)+2*diff(z(t),t)+z(t)=2*exp(-t):
odetest(z(t)=z0(t), ode);

                               
0

The return value 0 means that it's a solution.

According to its help, the 2nd argument to rgf_findrecur should be a list of 2*consecutive terms of the sequence, where K is the order of the recurrence (as well as the 1st argument). Presumably, extra terms are simply ignored.

Edit: I replaced with 2*K.

There's no point in using CodeTools:-Usage on some code that doesn't even work. To see that it doesn't work, remove the Usage command. The reason that Usage is obscuring that it doesn't work is because you've used output= bytesused, so the output is just that number rather than the result of your computation.

The vast majority of modules should not be packages. What benefit did you think you were getting by using option package? The only good reason (AFAIK) to use option package is if you want to load the module with with.

First 76 77 78 79 80 81 82 Last Page 78 of 395