Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

I agree with everything that Acer said above. I just thought that I'd state it more concisely and also list most of the low-level commands that print things as side effects.

Ordinarily (*footnote 1), what you see displayed are only the return values (*footnote 2) of the top-level procedures that you call. (This is even true for procedures that produce plots.) If a procedure has an untrapped error, then it has no return value, so all you see is the error message.

On the other hand, the following commands display things as a side effect, so that that displaying happens regardless of where in the procedure the command is executed: print, lprint, printf, and userinfo (that list does NOT include sprintf or nprintf). And many (most? all?) file output commands allow you to specify terminal as the file, which'll ordinarily (*footnote 3) cause display to the screen.

(*1) "Ordinarily" can be changed by printlevel, trace, or other debugging commands.

(*2) Return values of type function can be trapped and have their display format changed by procedures named `print/...`. In several very common cases (such as PLOT, Matrix, and Typesetting:-mrow) this trapping and reformatting is handled by an internal procedure. If you're debugging and suspect that this is happening, look at the output with lprint.

(*3) This "ordinarily" can be changed by writeto, i.e., things intended for the screen can be sent to a file instead.

Summarizing what happens in TSprintf: My procedure returns an object of type function whose function name is Typesetting:-mrow. That function name is trapped and the object is reformatted by the Standard GUI using internal code. So, my procedure actually returns a data structure rather than displaying something by side effect. Since the display is not done by side effect, it can be blocked by error.

Like this:

Sols3:= (H::algebraic, F::list(algebraic), i::posint, j::posint)->
   nops(
      select(
         s-> andmap(hastype, remove(evalb, s), hname),
         remove(
            s-> true in map(e-> (is(rhs(e) <= 0) assuming positive), s),
            [GTS(H, F, i, j)]
         )
      )
   )
:
(n,m):= (4,4); #Change to whatever size you want
Statistics:-HeatMap(Matrix((n,m), curry(Sols3, H, F)));

 

I do not know the formula for the polynomial. I guess that it depends on the function that's being interpolated? So, in the code below, you just need to fill that formula in where I put randpoly in the first line.

L:= (k::nonnegint, x::name)-> randpoly(x, degree= k):

make_pw:= proc(n::posint, interval::range(algebraic), x::name)
local 
   a:= op(1, interval), 
   h:= (op(2, interval) - a)/n, 
   k,
   x_k:= proc(k) option remember; a+k*h end proc
;
   unapply(piecewise(seq([x < x_k(k), L(k,x)/x_k(k)][], k= 1..n)), x)
end proc:

#applying it to your example:
f:= make_pw(10, 0..1, x);

#or, the same thing with decimal coefficients:
f:= make_pw(10, 0. .. 1., x);

The fs constructed above are procedures which when expanded, as in f(x), will be algebraic piecewise expressions if x is symbolic or just numbers if x is numeric.

Like this:

restart:
betas:= [2, 2.5, 5]:
colors:= table(betas =~ [black, red, blue]):
Eq1:= diff(y(x),x$2) + 'beta'*y(x) = 0, y(1) = sin(1), D(y)(1) = cos(1):
Eq2:= y = 'beta'*exp(x) + sin(2*x) + x^2:
dsol:= dsolve({Eq1}, numeric, parameters= ['beta']):
for beta in betas do
   dsol(parameters= [beta]);
   DP[beta]:= plots:-display(
      [plots:-odeplot(dsol, [x,y(x)], x= 0..0.1, legend= ['beta'(An)=beta]),
       plot(rhs(Eq2), x= 0..0.1, linestyle= dashdot, legend= ['beta'(Num)=beta])
      ],
      color= colors[beta]
   )
od:      
plots:-display(
   [entries(DP, nolist)],
   title= typeset("Variation of velocity for different ", 'beta')
);

 

As far as I know, the option implicit= true has nothing to do with the numeric use of dsolve; it only applies to symbolic solutions. What was your reason for using it? I wouldn't be surprised at all if some change made to dsolve's input processor between the two versions of Maple inadvertently caused it object to this input. Try removing it and see what happens. I can't test because you didn't provide your input or worksheet.

The symbol args has special meaning when it's used inside a procedure; it's essentially a reserved word. (Anything constructed with the arrow operator -> is a procedure, just as if it were constructed with proc() ... end proc.) So, you need to change args to something like Args.

Like this:

restart:
F:= (k::posint)-> assign(m||(1..2*k) =~ seq([dx||i, dy||i][], i= 1..k)):
F(3), m||(1..6);
                  dx1, dy1, dx2, dy2, dx3, dy3

 

A variable is only implicitly made local to a Maple procedure if an explicit assignment (using :=) is made to that variable within the procedure.

There are three (and only three) Maple commands that will accept a bound variable in their arguments and give it a kind of local status: seq, add, and mul. These are all built-in, not written in Maple. There's no way for someone writing a Maple procedure to give bound variables that will appear in their procedure's arguments that same status. This would be a great addition to the Maple language. Then the bound variables in numerous commands such as plot, int, diff, limit, etc. would be easier to use. As it stands, it's up to the procedure's user to make sure that bound variables are unassigned. This is a serious shortcoming of the whole Maple programming language.

Even the bound variables that are used with seq, add, and mul are not 100% local: If those variables are protected, it will cause an error. This is also a serious shortcoming, and one for which I see no excuse.

So, the answer to your final question is unfortunately Yes: If you're going to use a bound variable in a procedure and you want to be a careful programmer, then you need to declare that variable local, even if it's the bound variable for seq, add, or mul. An alternative that is practical and safe in the vast, vast majority of cases (but not 100% perfectly safe) is to begin the variable's name with an underscore.

The first method is overwhelmingly more efficient even if you don't have foreknowledge of nFrames. In Maple, there's no such thing as dynamically growing a list or set. If any change whatsoever is made to a list or set, then the entire list or set is reconstructed from scratch regardless of how large it is. But you can dynamically grow a table or Array.

Furthermore, there's no need to preallocate memory, although it doesn't hurt either. However, creating a list of zeros to initialize an Array in a complete waste. Arrays are essentially initialized to 0 anyway.

The issues in the first paragraph above are extremely important. Attempting to dynamically grow lists or sets may be the number-one source of inefficiency in Maple programs. The issues in the second paragraph are minor.

If you understand the notation of asymptotic time complexity: The first method is O(nFrames); the second method is O(nFrames^2). This is easy to demonstrate experimentally and easy to explain theoretically.

This issue has nothing to do with animation in particular. It is pervasive throughout Maple.

Upto my understanding of this subject, which is admittedly weak, the likelihood function is only properly defined for distributions with at least one symbolic parameter. That seems to be confirmed by both the Wikipedia article "Likelihood function" and the fact that your workaround using symbolic m worked. Given this, it makes sense that the result is always 1 when there are no parameters. In other words, What is the likelihood, or probability, that you've correctly estimated the parameters when there are no parameters? Of course it's 1.

What I'm wondering is What happened to the factors of 1/sqrt(2*Pi) that usually appear in the Normal PDF?

There's only one practical way to limit the time that a Maple command uses: Wrap the command in the appropriately named timelimit command. See ?timelimit. Using kernelopts(cpulimit= ...) for this is brutal, crude, and leads to lost work.

Further finesse is achieved by wrapping the timelimit command in a try-catch block. See the examples on the referenced help page.

Since you are a student of number theory, you should immediately prove the following three elementary propositions about divisors, which'll lead to one efficient method of solving your posted problem (efficient for the vast, vast majority of practical cases). I both came up with and proved these propositions as my first formal written proof as a 16-year-old college freshman.

1. The number of divisors of a positive integer that are less than its square root exactly equals the number that are greater.

2. The number of divisors of a positive integer n is odd iff n is a perfect square.

3. If a positive integer n has no divisors d such that 1 < d <= sqrt(n), then n is prime.

(The only criticism that I received on that freshman work was due to the fact that in my naivete, I called them "factors" rather than "divisors".)

I'll let that sit for a few hours, then I'll post two algorithms, one based on sort that uses proposition 1 above and one that uses min and select, akin to Mr Heinz's but much faster. Now, based on my theoretic knowledge of algorithmic asymptotic time complexity, I know that there's guaranteed to be some unimaginably HUGE number n such that the min/select method is faster than the sort method, but I'd guess that the smallest such n has millions of divisors or more. That why I said "practical cases" above.

I've read the paper that you linked the PDF of. The ODE system in the paper is a two-point BVP of four second-order ODEs. Maple's dsolve comes with four solvers dedicated to two-points BVPs. RK4 is an IVP solver (and, as I said above, it's been long ago superseded by several better IVP solvers, the most notable being RKF45). To use an IVP solver to solve a BVP, one must construct an elaborate "shooting method":

  1. let some of the initial conditions be chosen arbitrarily
  2. solve an IVP
  3. evaluate the solution at the other boundary point. If they're close enough, you're done; if not, refine the initial conditions and go back to step 2.

The authors only give the barest hint of this process by mentioning Newton's method. Since BVP solvers are much flakier than IVP solvers, the above process is occasionally necessary. In the case of this BVP, it most definitely is not. So, it seems to me that in lieu of any explanation for why they used RK4 + shooting + Newton---and they give none---the authors are foolish to use this elaborate method. They should provide you with their Maple code.

In the following worksheet, I easily and systematically generate the first 56 plots from the paper, that being all the plots in figures 2 - 5. I could've made the rest of the plots, but I got tired of it, and I think that you can and should follow my example, which'll easily get you the remaining umpteen 2D plots. If you care to generate the 3D plots also, I can help you with that just as systematically.

Note that my plots are much smoother than those shown in the paper, especially when you view them in Maple rather than on MaplePrimes. That's because the BVP solver is automatically choosing the interval width needed to make them smooth. It's all automatic! All the default dsolve settings worked perfectly (although that's often not the case with BVPs).

This BVP is homogenous (i.e., the independent variable does not appear) and nearly linear. The only nonlinearities are two quadratic terms neither of which involve the highest-order derivatives. That's about as easy as it gets for practical numeric BVPs that we see on MaplePrimes.
 

restart:

ODEs:= <
   diff(U(Y),Y$2) = -Ha*diff(B(Y),Y) - theta(Y) - Br*phi(Y) - V0*diff(U(Y),Y), #Eq. 8
   diff(B(Y),Y$2) = -(V0*diff(B(Y),Y) + Ha*diff(U(Y),Y))*Pm,                   #Eq. 9
   diff(theta(Y),Y$2) = -Nt*diff(theta(Y),Y)^2 -
      (V0*Pr + Nb*diff(phi(Y),Y))*diff(theta(Y),Y),                            #Eq. 10
   diff(phi(Y),Y$2) = -Nt/Nb*diff(theta(Y),Y$2) - V0*Sc*diff(phi(Y),Y)         #Eq. 11
>;
ICs:= <
   phi(0) = 1, B(0) = 0, U(0) = 0, D(theta)(0) = -1,                           #Eq. 12
   phi(1) = 0, D(B)(1) = 0, U(1) = 0, theta(1) = 0                             #Eq. 12
>;
param_names:= [V0, Ha, Pm, Nt, Nb, Br, Sc, Pr];
J:= -D(B):                                                                     #Eq. 13

ODEs := Matrix(4, 1, {(1, 1) = diff(diff(U(Y), Y), Y) = -Ha*(diff(B(Y), Y))-theta(Y)-Br*phi(Y)-V0*(diff(U(Y), Y)), (2, 1) = diff(diff(B(Y), Y), Y) = -(V0*(diff(B(Y), Y))+Ha*(diff(U(Y), Y)))*Pm, (3, 1) = diff(diff(theta(Y), Y), Y) = -Nt*(diff(theta(Y), Y))^2-(V0*Pr+Nb*(diff(phi(Y), Y)))*(diff(theta(Y), Y)), (4, 1) = diff(diff(phi(Y), Y), Y) = -Nt*(diff(diff(theta(Y), Y), Y))/Nb-V0*Sc*(diff(phi(Y), Y))})

 

ICs := Matrix(8, 1, {(1, 1) = phi(0) = 1, (2, 1) = B(0) = 0, (3, 1) = U(0) = 0, (4, 1) = (D(theta))(0) = -1, (5, 1) = phi(1) = 0, (6, 1) = (D(B))(1) = 0, (7, 1) = U(1) = 0, (8, 1) = theta(1) = 0})

 

[V0, Ha, Pm, Nt, Nb, Br, Sc, Pr]

(1)


Maple BVPs do not accept parameters in the same way that I showed you before, which was for IVPs. You need to dsolve for each new set of parameters. Here's a procedure Solve to do that. I've taken the default parameter values from the captions of Figs. 2-4.

Solve:=
   subs(
      _P= param_names,
      proc({
         V0::realcons:= 1,
         Ha::realcons:= 5,
         Pm::realcons:= 0.8,
         Nt::realcons:= 0.1,
         Nb::realcons:= 0.1,
         Br::realcons:= 1,
         Sc::realcons:= 1,
         Pr::realcons:= 10
      })
         userinfo(1, Solve, param_names=~ _P);
         dsolve(
            eval(
               convert(ODEs, set) union convert(ICs, set),
               param_names=~ _P
            ),
            numeric
         )
      end proc
   )
;
   

proc ({ Br::realcons := 1, Ha::realcons := 5, Nb::realcons := .1, Nt::realcons := .1, Pm::realcons := .8, Pr::realcons := 10, Sc::realcons := 1, V0::realcons := 1 }) userinfo(1, Solve, `~`[:-`=`](param_names, ` $`, [V0, Ha, Pm, Nt, Nb, Br, Sc, Pr])); dsolve(eval(`union`(convert(ODEs, set), convert(ICs, set)), `~`[:-`=`](param_names, ` $`, [V0, Ha, Pm, Nt, Nb, Br, Sc, Pr])), numeric) end proc

(2)

infolevel[Solve]:= 1:

Now I generate the plots in Fig. 2 (changing values of V0):

P:= V0:
vals:= [0.2, 0.4, 0.6, 1.0]:
sols:= [seq(Solve(P= v), v= vals)]:
colors:= [red, green, blue, purple]:
for F in [U,B,J,theta,phi](Y) do
   print(plots:-display(
      [seq(
         plots:-odeplot(sols[k], [Y,F], color= colors[k], legend= [P= vals[k]]),
         k= 1..nops(vals)
      )],
      labeldirections= [horizontal,vertical]
   ))
od:

Solve: [V0 = .2 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = .4 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = .6 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = 1.0 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]

 

 

 

 

 

 

Fig. 3 (changing values of Ha):

P:= Ha:
vals:= [1, 5, 10, 20]:
sols:= [seq(Solve(P= v), v= vals)]:
colors:= [red, green, blue, purple]:
for F in [U,B,J](Y) do
   print(plots:-display(
      [seq(
         plots:-odeplot(sols[k], [Y,F], color= colors[k], legend= [P= vals[k]]),
         k= 1..nops(vals)
      )],
      labeldirections= [horizontal,vertical]
   ))
od:

Solve: [V0 = 1 Ha = 1 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]

Solve: [V0 = 1 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 10 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 20 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]

 

 

 

 

Fig. 4 (changing values of Pm):

P:= Pm:
vals:= [0.02, 0.2, 0.4, 0.8]:
sols:= [seq(Solve(P= v), v= vals)]:
colors:= [red, green, blue, purple]:
for F in [U,B,J](Y) do
   print(plots:-display(
      [seq(
         plots:-odeplot(sols[k], [Y,F], color= colors[k], legend= [P= vals[k]]),
         k= 1..nops(vals)
      )],
      labeldirections= [horizontal,vertical]
   ))
od:

Solve: [V0 = 1 Ha = 5 Pm = 0.2e-1 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 5 Pm = .2 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 5 Pm = .4 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1 Sc = 1 Pr = 10]

 

 

 

 

Fig 5 (changing values of Br):

P:= Br:
vals:= [0.1, 0.5, 1.0, 2.0]:
sols:= [seq(Solve(P= v), v= vals)]:
colors:= [red, green, blue, purple]:
for F in [U,B,J](Y) do
   print(plots:-display(
      [seq(
         plots:-odeplot(sols[k], [Y,F], color= colors[k], legend= [P= vals[k]]),
         k= 1..nops(vals)
      )],
      labeldirections= [horizontal,vertical]
   ))
od:

Solve: [V0 = 1 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = .1 Sc = 1 Pr = 10]

Solve: [V0 = 1 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = .5 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 1.0 Sc = 1 Pr = 10]
Solve: [V0 = 1 Ha = 5 Pm = .8 Nt = .1 Nb = .1 Br = 2.0 Sc = 1 Pr = 10]

 

 

 

 

 


 

Download nanofluid_BVP.mw

[Before you download the above, note that there's an updated version in a Reply below. The update includes everything that's in the original.]

Vote up for a great Question.

The symbols &under and `&+` are defined on the page ?type,structured, which is one of deepest pages in whole help system and one of my most commonly reread (the other being ?operators,precedence). I must add that both &-operators and their use in structured types has been around for a long time---at least 18 years, which covers the extent of my usage.

The & is being used in two different ways above. Let me first explain &under. Any symbol (with a few exceptions to disallow ambiguity or excessive weirdness) that begins with & can be used as a binary infix operator, regardless of whether the symbol's meaning has been predefined. If I type A &foo B, then the function `&foo`(A,B) is constructed. If a procedure `&foo` has been defined, then it is called; otherwise the construction remains an unevaluated function like any other. Note that the quotes are required for the prefix form `&foo`(A,B) and that they are not used in the infix form A &foo B, and these two expressions mean exactly the same thing.

So, the call type(e, T &under F) becomes type(e, `&under`(T,F))---a parameterized type. Like any old-style (i.e., pre-TypeTools) parameterized type that doesn't have a built-in handler, this invokes `type/&under`(e, T, F). Use showstat to see the one-line code of `type/&under`. Then you'll see that the original type command is equivalent to type(F(e), T). So why would one use the former? Because now T &under F is itself a type that can be used with the numerous fundamental commands that require a type argument, such as indets, hastype, membertype, procedure parameter restrictions, etc.

(Typing up explanation for `&+` now. Need to post so the above doesn't get lost. One wrong keypress in MaplePrimes and your whole post is gone...

Okay, here it is....)

Note very carefully that & is not itself an operator! It's just a character and has no semantics by itself. Any characters---no matter how weird (including control characters, spaces, escape characters, unicode, quotation marks, etc.)---can be part of a Maple symbol as long as the symbol is enclosed by ` `. I've worked on a large system where all the internally generated symbols (mostly dummy variables of integration) were Chinese characters. If the symbol's first character is &, and the rest are normal printing characters, then the symbol can used infix (there are also some ambiguous constructions that are not allowed).

The purpose of `&+` is to have something that looks kind of like `+` but has no procedure assigned to it. If you were to instead use `+`(...), then the addition would immediately evaluate before the type command saw it. The command type(e, `&+`(T1, T2)) returns true iff e is the addition of something whose first operand is type T1 and whose second operand is type T2 (and they must be in that order). The type `&+` is built-in to the type command; you can't see its code.

The type system is the essence of symbolic programming.

I think that "unknown" is a poor and confusing choice of words for this. I think that the term anonymous procedure is the most-used term for this concept in computer science, and even in the Maple documentation.

Of course the system knows what procedure produced the error. The issue is that in Maple it's easy and common to create and use a procedure without ever giving it a name. It's only the procedure's name that's unknown, not the procedure itself.

For example:

map(x-> sin(1/x), [0, 1]);

Here, the procedure x-> sin(1/x) is anonymous, and the error message is like what you show.

To make matters worse, Maple actually gives names to these anonymous procedures. They're all of the form `unknown/H` where H is some string of hexadecimal digits. If this name were shown in the error message, it would make things so much easier! Then you'd be able to view the errant procedure by applying eval to its name. I'm try to figure out a way to get the name. I know it'll involve using debugopts(callstack), but I haven't worked out the details yet.

First 176 177 178 179 180 181 182 Last Page 178 of 395