Carl Love

Carl Love

28130 Reputation

25 Badges

13 years, 310 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@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?

@Mariusz Iwaniuk You need to change method= classical to method= classical[rk4]. Without this, it's using Euler's method.

@taro Here's a worksheet showing a way to do what you asked.

restart:

CC:= Array((1..2)$3, ()-> <args>.10^~<2,1,0>, order= C_order):

C1:= ArrayTools:-Alias(CC, [rtable_num_elems(CC, All)], row);

C1 := Vector[row](8, {(1) = 111, (2) = 112, (3) = 121, (4) = 122, (5) = 211, (6) = 212, (7) = 221, (8) = 222}, order = C_order, attributes = [source_rtable = (Array(1..2, 1..2, 1..2, {(1, 1, 1) = 111, (1, 1, 2) = 111, (1, 2, 1) = 112, (1, 2, 2) = 112, (2, 1, 1) = 121, (2, 1, 2) = 121, (2, 2, 1) = 122, (2, 2, 2) = 122}, order = C_order))])

(1)

CF:= rtable(CC, order= Fortran_order):

C2:= ArrayTools:-Alias(CF, [rtable_num_elems(CF, All)], row);

C2 := Vector[row](8, {(1) = 111, (2) = 211, (3) = 121, (4) = 221, (5) = 112, (6) = 212, (7) = 122, (8) = 222}, attributes = [source_rtable = (Array(1..2, 1..2, 1..2, {(1, 1, 1) = 111, (1, 1, 2) = 111, (1, 2, 1) = 121, (1, 2, 2) = 121, (2, 1, 1) = 211, (2, 1, 2) = 211, (2, 2, 1) = 221, (2, 2, 2) = 221}))])

(2)

 

Download rtable_order.mw

In both cases that it's used, ArrayTools:-Alias produces an inplace copy of a 2x2x2 parent Array as an 8-element row Vector. (ArrayTools:-Alias can only be used to make inplace copies.) The rtable(..., order= ...) command makes a true copy of the parent.

It may be possible to make an inplace copy in the other order by using an indexing function. I've never done that, and I have some doubts that there'd ever be an efficiency benefit from doing that. From reading the Maple help, you may think that there are some commands that can change the order inplace. No, these commands just change the order attribute without changing the actual order. I mean that they simply change---by a trivial relabeling---the value of the order attribute.

Please make sure that you understand the difference between inplace copy and (true) copy. See ?copy, and note that the examples that are shown there for tables apply to rtables also.

Please read these help pages: ?worksheet,managing,exporttolatex and ?latex,viewing. If you then still need help, please upload the relevant files. Use the green up-arrow on the toolbar of the MaplePrimes editor to upload:

  1. the worksheet (either as .mw or .mws);
  2. the LaTeX (you'll need to rename .tex as .txt in order for MaplePrimes to allow it);
  3. the PDF (as .pdf).

And please describe exactly what you did to produce file 3 from file 2.

@taro Maple is designed so that the effect of order isn't usually apparent to the casual user. You may only notice it if you use ArrayTools:-Alias, which'll give very wrong results if you don't respect the order, or convert(..., list), which'll just give different, but not wrong, results.

As you may be able to guess from the names of the option's values, C_order and Fortran_order, this issue is much older than Maple itself.

@Kertab24 Transforming multiple definite integrals to other coordinate systems, such as polar, is not something that's handled well, if at all, by computer algebra systems, such as Maple. Transforming the integrand (the sqrt(x^2+y^2) in your case) is relatively easy, and they can generally handle that. The hard part is transforming the limits of integration, which come from the boundaries of the region of integration (the rectangle in your case). As I think you realize, this requires some visualization skill, or, as you said, "imagination".

So that I don't go too far with explanations that you may not understand, please give me some background on your education:

  1. Do you have any familiarity with (or recollection of) polar coordinates?
  2. How long has it been since you've formally taken a graded course in calculus? I don't mean just reading a refresher.
  3. What about multi-variable calculus, which is where double integrals are covered? For science/engineering/math majors, this would usually be the third-semester course in calculus.
  4. Do you still have the textbook for that course? If so, give me the publication details.

Feel free to send me a private email by selecting "Contact author" from the "More..." pull-down menu at the botton right of this message. I offer private online calculus tutoring.

First 337 338 339 340 341 342 343 Last Page 339 of 711