Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

@Adam Ledger The issue that leads to the error is that the bound (or "dummy") variables of integration are being considered by NonlinearFit as if they were model parameters to be fitted. Hence it tries substituting numeric values for them. This issue is easy to fix by turning the model into an arrow procedure whose procedure parameters are the independent variable, x, followed by the true model parameters, beta and `ΔA`. This removes the error, but actually getting an answer is another matter.

@nm I'm sorry, but you've mistaken my ad hoc workaround Answer above as being more general than I intended. By ad hoc, I mean that it was intended to apply to the specific plot that you originally posted, on which the original scale for all three axes was 0..10. It should be obvious that that particular 0.4-factor adjustment cannot work for a general plot. And by workaround, I mean something that's meant to bypass a shortcoming or bug of Maple until better library code is written; I don't mean something that's intended to give Maple functionality equivalent to Mathematica in some area.

It'd be fairly easy for me to write a procedure to rescale the lengths of all the axes of a general plot to visually conform to any "box ratios" that you wanted. But rescaling the tickmarks of an unseen Maple plot is hellishly complex because a Maple programmer has no access to the tickmarks of the original plot when they are generated by the plot renderer in the GUI, which is the case for the vast majority of plots. In other words, a plot that a user sees on their screen has more information, more detail, than a Maple programmer has access to.

This is a shortcoming of Maple. Mathematica handles it much better. You'll get no argument from me about that.

@Carl Love You also asked if there was a way that was even more efficient than your first. Yes:

nFrames:= 10:
plots:-display(
   [seq(plot(sin(i*t), t= -2*Pi..2*Pi), i= 0..nFrames-1)],
   insequence
);

There are two separate simplifications that I've done: I replaced for with seq, and I eliminated the temporary variables frames and w (which become garbage that needs to be collected). When creating sequences, lists, or sets iteratively, for can be replaced by seq whenever

  1. the element being created doesn't depend on any of the previous elements, and
  2. the condition for ending the sequence doesn't depend on any of the previous elements.

The efficiency difference between your first and second methods is huge for large nFrames. The additional efficiency achieved by using my variation above is tiny by comparison.

@Adam Ledger Yes, your observation about the implication of the coprimality of j and k is correct. But even thinking about the simple concept of coprimality is already getting more complicated than is necessary for the problem.

@vv Yes, of course, his algorithm is fine (and perhaps asymptotically the best). It's only its implementation in Maple that's bad. In Maple, one shouldn't use sqrt on integers with a huge number of divisors, use is on large-number exact inequalities, or use arrow expressions containing function calls that are constant with respect to the arrow's parameters (unless one is sure that there's a remember table in effect). The first two of these are Maple-specific pitfalls; the third is something that a competent programmer in any language should foresee (though it's easy in Maple to overlook the fact that an arrow expression is a procedure like any other).

I was quite surprised at how brutal on my system and unstoppable his function was when applied to my data. And I tried several times, each time using smaller and fewer data.

@Carl Love 

Here are the two algorithms that I was talking about, plus Heinz's:
 

restart:

I want to measure the times of the algorithms without measuring the time required by numtheory:-divisors. Since its time is common to both, it's irrelevant to comparing these two algorithms. To do this, I use a remember table.

Divisors:= proc(n::posint)
option remember;
   numtheory:-divisors(n)
end proc:

Algorithm using sort:

MidDivisor_sort:= (n::posint)-> (D-> D[1+iquo(nops(D),2)])(sort([Divisors(n)[]])):

Algorithm using select:

MidDivisor_select:= (n::posint)-> min(select(`>`, Divisors(n), isqrt(n))):

Heinz's algorithm is the same as the second, but his implementation of it is horrible because he uses the highly symbolic and very slow commands is and sqrt (without evalf).

Heinz:= n-> min(select(d-> is(d=n or d>sqrt(n)), Divisors(n))):

Generate some random test data: A list of numbers that have a lot of divisors.

E:= rand(0..19):

N:= ['rand()*mul([2, 3, 5, 7]^~['E()'$4])' $ 64]:

Load the remember table:

Divisors~(N):

Let them run!

gc(); L1:= CodeTools:-Usage([MidDivisor_sort~(N), gc()]):

memory used=71.81MiB, alloc change=0 bytes, cpu time=1.36s, real time=1.25s, gc time=625.00ms

 

gc(); L2:= CodeTools:-Usage([MidDivisor_select~(N), gc()]):

memory used=463.00MiB, alloc change=0 bytes, cpu time=8.48s, real time=5.62s, gc time=4.78s

 

Warning: Do not try to test Heinz's code on this data! You will crash Maple. His code is only good for baby numbers.

 

Verify that the results are the same:

evalb(L1 = L2);

true

(1)

 

Download MidDivisor.mw

Now, that stuff about asymptotic time complexity: Let ND be the number of divisors of n. It's clear that the select does a linear scan of length ND, and the min does a linear scan of length ND/2. The isqrt is sublinear (that's not so clear, but it's clear that it's done only once for each n). So the second algorithm is O(ND). Other than hand, the first algorithm does a sort, so it must be at least O(ND*log(ND)). That's why I know that it's guaranteed that there's some n (indeed, a whole infinitude of them) such that the select algorithm is faster than the sort algorithm. But those n are likely incredibly large, since Maple's sort is very fast.

Some readers may ask Why use sort at all? Aren't sets of integers sorted always already in Maple? Yes, that's certainly true for word-sized integers at this point, and perhaps all integers. But it's bad programming practice to rely on it. So, that's a problem (admittedly, a very small problem), with Tom's "break" algorithm.

@Gabriel samaila 

It's easy to include extra boundary values with exactly one new parameter each in the dsolve command (only for a BVP). The command will solve these for you, filling in the numeric values of the parameters. Since the skin friction = C__f and Nussel = Nu are defined this way (Eq. 14), I included them and produced Tables 1-3. For Tables 2 and 3, I got a perfect match with the paper down to the last digit. My Table 1 is incorrect because they only specify 4 of the 8 required parameters. I couldn't find what values they used for the others! Once you find these 4 values (by email maybe), you can easily fill them into my code. Here's the updated worksheet. It includes everything from the previous worksheet also, so you might as well overwrite your existing copy.

nanofluid_BVP.mw

@Adam Ledger No references are needed; this can be proved from first principles. But some concrete visualization helps:

Imagine a long high school corridor lined with lockers along the left wall for as far as the eye can see. The lockers are consecutively numbered starting with 1. The numbers are clearly displayed on the doors. They go up to some n >= 999. The exact n doesn't matter; pick whatever reasonable number helps you visualize the situation. If you can't decide, use 1729 (It doesn't have any special properties related to this problem. I just like it for unrelated reasons.) It is the last day of the school year. All of the locker doors start out open. The n students are exiting, single file, walking past the lockers. The first student closes the door of every locker. The second student then opens the door of every second locker starting with locker 2. The third student changes the state of every third door starting with locker 3. The kth student changes the state of every kth door starting with locker k. Finally, the last student just changes the state of the last door. Now, which locker doors are closed?

@Adam Ledger Oh, I didn't mean that I'd post the solutions in a few hours. I meant that I'd post the algorithms, only one of which uses one of the propositions. In this case, seeing the algorithm won't provide any hint towards proving the proposition. It would give a hint towards formulating the proposition, but I already did that completely above.

@Gabriel samaila So why did my the vote up on my original Answer disappear, and none has appeared on this one? I'm not saying that you necessarily had anything to do with that.

So, does the paper pass or fail your review? I'd fail it unless they make a few changes:

Provide the Maple code

and either

Justify the method that they use and explain the "Newton" part OR Use a more-modern more-automatic method.

@vv Yes, I wondered that myself! It'd be nice to see an example where it's used in library code.

@Kertab24 

The Edwards & Penney calculus textbook is good. I've used it, and I have a copy (different edition). Do you still have that textbook? So you've finished a semester of calculus using this book. That would be roughly the first 5 chapters, right? So, you've barely scratched the surface of integration, right? You've covered u-substitution, right? But that's about where you stopped, right?

So, summarizing and interpolating a bit from what you said: You're several years into an undergraduate chemical engineering program, which requires three semesters of calculus and likely also one semester of differential equations (essentially a fourth semester of calculus), yet you only have one semester of calculus, which you did well in. How did that happen? I mean, why did they let you bypass the prerequisites? Have you had two semesters of physics (one of mechanics and one of electricity and magnetism)?

Maple is very good at both symbolic and numeric integration, and I think that it can be a great help to you in your program. But I think that you'll need to know at some point soon how to at least set up double and triple integrals in polar, cylindrical, and spherical coordinates. Maple should be able to evaluate the vast majority of them if they're set up correctly. And as the Answer by VV shows, Maple can sometimes do, without coordinate changes, integrals that almost anyone doing by hand would use a coordinate change for. But my two-year-old Maple couldn't do it.

@Gabriel samaila 

There's no need to manually convert the system to first-order equations. The dsolve command handles that automatically.

Here's a complete solution to your problem with a comparison between RK4 and rkf45 (the default method).
 

restart:

dsys:= {
   diff(x[1](y),y) = 1,
   diff(x[2](y),y) = x[3](y),
   diff(x[3](y),y) = -x[6](y)-H[a]*x[3](y)-V[0]*x[3](y)-B[r]*x[8](y),
   diff(x[4](y),y) = x[5](y),
   diff(x[5](y),y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y),
   diff(x[6](y),y) = x[7](y),
   diff(x[7](y),y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y),
   diff(x[8](y),y) = x[9](y),
   diff(x[9](y),y) =
      -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b] +
      N[t]*P[r]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y)
}:

ICs:= {
   x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0, x[5](0) = 0,
   x[6](0) = 0, x[7](0) = -1, x[8](0) = 1, x[9](0) = -.68
}:

Params:= [
   B[r] = 1, H[a] = 5, N[b] = .1, N[t] = .1,
   P[m] = .8, P[r] = 10, S[c] = 1, V[0] = 1
]:

#Method 1: Evaluate the system with respect to the parameters
#before calling dsolve.

Sol1:= dsolve(eval(dsys union ICs, Params), numeric):

MyPlot:= Sol->
   plots:-odeplot(
      Sol1, [seq([y, x[k](y)], k= 1..9)], y= 0..2,
      linestyle= [seq([solid, dash, dot][], k= 1..3)],
      color= [red $ 3, green $ 3, blue $ 3],
      legend= [seq(x[k], k= 1..9)],
      size= [2000,1000]
   )
:

MyPlot(Sol1);

 

#Method 2: Specify the numeric value of the parameters after calling dsolve.
Sol2:= dsolve(dsys union ICs, numeric, parameters= lhs~(Params)):

Sol2(parameters= Params):

MyPlot(Sol2);
#The plot is necessarily identical to the above, and I omitted it here.

#Compare with RK4:
Sol3:= dsolve(
   dsys union ICs, numeric, parameters= lhs~(Params),
   method= classical[rk4], stepsize= 0.1
):

Sol3(parameters= Params):

MyPlot(Sol3);

 

#Numeric comparison at y=1:
RKF45:= eval[recurse](Sol2(y), [y= 1])[2..];

[x[1](1) = HFloat(1.0), x[2](1) = HFloat(-0.058264834925667204), x[3](1) = HFloat(-0.0010873159899245078), x[4](1) = HFloat(0.10808962402573118), x[5](1) = HFloat(0.14658764048208384), x[6](1) = HFloat(-0.10159131480245949), x[7](1) = HFloat(-5.077149160781932e-5), x[8](1) = HFloat(-0.020372829806255117), x[9](1) = HFloat(-0.6595763987021374)]

(1)

RK4:= eval[recurse](Sol3(y), [y= 1])[2..];

[x[1](1) = HFloat(1.0), x[2](1) = HFloat(-0.05826577858436468), x[3](1) = HFloat(-0.0010813843558394268), x[4](1) = HFloat(0.10808999748780099), x[5](1) = HFloat(0.14659111634721797), x[6](1) = HFloat(-0.10158856147456753), x[7](1) = HFloat(-6.078236913125798e-5), x[8](1) = HFloat(-0.02037430234974795), x[9](1) = HFloat(-0.6595649152811208)]

(2)

#The difference:
lhs~(RK4) =~ rhs~(RK4) -~ rhs~(RKF45);

[x[1](1) = HFloat(0.0), x[2](1) = HFloat(-9.43658697474814e-7), x[3](1) = HFloat(5.931634085080918e-6), x[4](1) = HFloat(3.734620698109259e-7), x[5](1) = HFloat(3.475865134133782e-6), x[6](1) = HFloat(2.753327891957813e-6), x[7](1) = HFloat(-1.001087752343866e-5), x[8](1) = HFloat(-1.4725434928329617e-6), x[9](1) = HFloat(1.1483421016533768e-5)]

(3)

#It looks like the RK4 results are pretty good.


 

Download How_to_dsolve.mw

@tomleslie Is there any other way to export a worksheet file as LaTeX other than to use the File -> Export As menu?

First 335 336 337 338 339 340 341 Last Page 337 of 709