Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

@Ida2018 Sorry, I should've been more precise. I meant replace O__L by O__L*Unit(kN) in the equation. The other O__L, which is the second argument to solve, should remain as is.

@Joe Riel I think that time would be better spent by first doing obvious optimizations (obvious to you and me, that is) than in analyzing a profile or dump of a program that is obviously so raw. But first, I want to make sure that it works.

@kfli Read? Ha! I found it out by experimenting with your code for about 90 minutes. You can confirm what I said by doing this little experiment:

A:= Matrix():
whattype(A);

                             Matrix
print(whattype(A));
                      proc()  ...  end;

This is the same difference (and for exactly the same reason) as can be observed in

P:= proc() end proc:
P;
print(P);

That difference is described on help page ?last_name_eval.

Your procedure is great. I see a great many ways that it can be improved, and I look forward to working on that. But before I work on anything other than stylistic issues, I want to make sure that it works. Why spend time optimizing something that doesn't even work? You'd be the best person to say whether it works.

This little code fragment seems problematic to me:

R := norm(df,2);
Q := R/~df;
gamma := R/~Transpose(Q).fval;

Please tell me in purely algebraic terms (by which I mean without making any reference to the rest of the algorithm or to what df or fval represent) what you think it does. All that it actually does is gamma:= Transpose(df).fval (possibly the round-off errors are different).
 

 

@waseem I accept your apology. Thank you.

Thank you for the new paper. Your code shows that you have obviously learned a tremendous amount from reading my MaplePrimes posts. It's quite impressive.

Do you have any information at all that could help me understand equations 23 - 31? Obviously, the authors (Gupta, et al) have written some code. Why don't they supply it? Given any code in any language, I could translate it to Maple and improve it tremendously.

The notions of error control mentioned in the Gupta paper are ridiculous. Basically, they have said something akin to "For one set of parameter values, we got convergence to 4 decimal places using 160 nodes; therefore, we use 160 nodes for all sets of parameter values." and "For one set of parameter values, we got convergence to 4 decimal places by setting infinity to 8; therefore, we use infinity = 8 for all sets of parameter values." No self-respecting mathematican would say those things.

The Maple BVP solver has very fine error control. If it detects a failure to converge, it will report that. If you ask for 4 digits and it thinks that it can only get 3, it will report that. So, if it produces a numeric result, I trust it. At this point, I have matched the skin friction numbers in the first column of Table 3 to 1-2 decimal places. Thus, I assume that the numbers reported by Gupta, et al, in that table are not accurate to 4 decimal places and the Maple numbers are correct.

I am interested in coding their FEM, but I can't understand Eqs 23-31. When I code it, I will dynamically adjust the mesh size on every run and make sure that it converges.

 

@666 basha You say that "the problem can be solved very easily using RK method dsolve." Well, I'd certainly like to see that. Please show me a solution technique (of either the corrected or uncorrected system) by whatever method you choose to use. Having any solution (for comparison) will help with coding and testing the authors' FEM method.

Before the problem can be solved by any method, it must be input correctly into that method. At this point, I am absolutely certain that the system as typed into the paper is incorrect (in both places that it appears in the paper), and I believe that Preben is also certain of that. I think that the only error is that the sign of lambda needs to be changed, but I am not certain of that.

If I had plunged ahead with coding the paper's FEM, and I'd assumed that the system was correct as published, then there would have been enormous and lengthy frustration trying to find a non-existent error in the FEM code. Please tell me: Can you understand such a fundamental principle of engineering? 

If you have some understanding of equations 23 - 31, please share that. I am still trying to understand them. To begin with, how are the mesh nodes eta[e] (which appear as the limits of integration) chosen? Are they simply evenly spaced over the eta interval (0 to 8 in this case)? 

@Ritti For real x, it is obvious that abs(cosh(x)) = cosh(x). If this is not obvious to you, then you should take a step back from the problem that you're currently working on and spend a few minutes reviewing the definition of cosh. Any decent calculus textbook will help, as will the Wikipedia article "Hyperbolic function".

@vv By "sequential order", the OP is referring to tree depth within the expression, not tree breadth. In other words, they want the innermost function (or perhaps functions) with the desired type of arguments, not the rightmost.

@Ritti Just ask Maple again:

is(abs(cosh(x)) <= exp(x^2/2)) assuming x::real;
      true

@zhuxian This is also easy with basic Maple commands. Suppose that is the expression that contains the constants and _C7 is the only one that shouldn't be set to 0. Then do

eval(S, (indets(S, suffixed(_C)) minus {_C7})=~ 0);

When you enter a Question, isn't there a box labelled Tags immediately below the box where you enter the text of the Question?

I reposted your actual math/Maple question as an official Question.

I used my moderator privileges to turn your question into an official Question (which a moderator can do even without tags). Do you see a box +Add Tags  under the word Posted on the left side of the header, to the lower right of your avatar? Try clicking on that.

@Rohith The new example can be handled with a relatively small change to the filter:

value(
   indets(
      InertForm:-Parse("a+(c/d)*log(1+b)+sin(a-b)+a^c"), #string!
      #selects operators with exactly two operands, both of which are names,
      #negated names, integers, or floats:
      And(specfunc({name, (-1) &* name, numeric}, {`%+`, `%*`, `%/`, `%^`}), 2 &under nops) 
   )                                            
);

So, the only changes that I made were adding (-1) &* name to the allowed operands and `%^` to the allowed operators. So, the items selected by this new filter are necessarily a superset of those selected by the previous one since both the operand set and the operator set have been enlarged.

Regarding your terminology: You should say that you're extracting "operators acting on exactly two atomic operands". What you mean will then be much clearer (even though this doesn't conform precisely with the technical definition of atomic used by Maple), because operand alone can apply to subexpressions of any complexity. For example, in a + b/c the `+` has exactly two operands: a and b/c.

@Preben Alsholm If an operator is restricted to having exactly two operands, then sometimes lhs and rhs are used to refer to those operands, even when the operator is not =, <<=<>. This is the sense by which the OP is using the terms "lhs" and "rhs"

@waseem Do not repeat the same inane question over and over again. It makes me very angry. Are you less than 10 years old? Next time, I'll just delete it. Be patient: It's not like I get paid to answer questions on MaplePrimes. It is just my hobby. As you can see, I've been working on the problem, but I got stuck. That's in addition to the numerous other MaplePrimes problems that I'm working on.

I'd suggest that if you're asking for help from someone whom you're not paying, that you limit any prodding to at most once per week, and that you learn to do it in a polite respectful way.

 

I am trying to solve the system with Maple's built-in numeric BVP solver (which is using a nonlinear FEM with dynamic mesh size anyway). The system solves very easily, but the values that I am getting do not agree with the paper at all, not even to one significant digit! I feel that I must have entered something wrong, like a minus sign perhaps. But I have gone over what I typed "with a fine-tooth comb" at least 10 times, comparing my worksheet with the paper character by character, using the 3 places in the paper where the ODEs are given and the two where the BCs are given.

So, I'd like some assistance to find the mistake. Or, at least I'd like someone to tell me that I've entered the system and parameter values exactly as they're given in the paper.
 

restart:

Digits:= 15:

ODEs:= [
   #Eq 9/18/24:
   (1+K)*diff(f(eta),eta$3) + f(eta)*diff(f(eta),eta$2) - 2*n/(n+1)*diff(f(eta),eta)^2 -
   2/(n+1)*M*diff(f(eta),eta) + K*diff(g(eta),eta) + 2/(n+1)*sigma*theta(eta)  =  0,
   #Eq 10/19/25:
   (1+K/2)*diff(g(eta),eta$2) + f(eta)*diff(g(eta),eta) -
   (3*n-1)/(n+1)*diff(f(eta),eta)*g(eta) - 2/(n+1)*K*(2*g(eta)+diff(f(eta),eta$2))  =  0,
   #Eq 11/20/26:
   diff(theta(eta),eta$2) + Pr*f(eta)*diff(theta(eta),eta)  =  0
]:   

<ODEs[]>;

Vector(3, {(1) = (1+K)*(diff(diff(diff(f(eta), eta), eta), eta))+f(eta)*(diff(diff(f(eta), eta), eta))-2*n*(diff(f(eta), eta))^2/(n+1)-2*M*(diff(f(eta), eta))/(n+1)+K*(diff(g(eta), eta))+2*sigma*theta(eta)/(n+1) = 0, (2) = (1+(1/2)*K)*(diff(diff(g(eta), eta), eta))+f(eta)*(diff(g(eta), eta))-(3*n-1)*(diff(f(eta), eta))*g(eta)/(n+1)-2*K*(2*g(eta)+diff(diff(f(eta), eta), eta))/(n+1) = 0, (3) = diff(diff(theta(eta), eta), eta)+Pr*f(eta)*(diff(theta(eta), eta)) = 0})

BCs:= [
   #Eq 12/21:
   f(0) = -lambda, D(f)(0) = -1, g(0) = -1/2*(D@@2)(f)(0), D(theta)(0) = -c*(1-theta(0)),
   #Eq 13/22:
   D(f)(%infinity) = 0, g(%infinity) = 0, theta(%infinity) = 0
   #I changed infinity to %infinity to make it a parameter.
];

[f(0) = -lambda, (D(f))(0) = -1, g(0) = -(1/2)*((D@@2)(f))(0), (D(theta))(0) = -c*(1-theta(0)), (D(f))(%infinity) = 0, g(%infinity) = 0, theta(%infinity) = 0]

Solve:= subs(
   _BVP= {ODEs[],  BCs[]},
   proc(
      #The values presented here are simply defaults. They can be changed when the
      #procedure is called.
      {K::realcons:= 0, Pr::realcons:= 0.733, sigma::realcons:= 0, lambda::realcons:= 3,
       c::realcons:= 1, M::realcons:= 1, n::realcons:= 1, %infinity::positive:= 8
      }
   )
      dsolve(_BVP, numeric, _rest)      
   end proc
):
   

Try to duplicate Table 2,  second column:

Table2:= Solve(K= 0, Pr= 0.733, sigma= 0, lambda= 3, c= 1, M= 1, n= 1):

seq(eval(diff(f(eta),eta), Table2(Eta)), Eta= 0..8);

HFloat(-0.9999999999999999), HFloat(-1.0), HFloat(-1.0000000000000004), HFloat(-1.0000000000000002), HFloat(-0.9999999999999994), HFloat(-0.9999999999991915), HFloat(-0.9999999968802413), HFloat(-0.9999664933149132), HFloat(0.0)

Totally wrong!

 

Use known exact solution (Eq 32) for this special case to try to duplicate Table 2, first column:

#Eq 32:
F:= eval[recurse](lambda - (1-exp(-eta*z))/z, [z= (lambda+sqrt(lambda^2+4*M-4))/2, lambda= 3, M= 1]);

3-(1-exp(-eta*(3/2+(1/2)*9^(1/2))))/(3/2+(1/2)*9^(1/2))

evalf(map2(eval,  diff(F,eta), eta=~ [$0..8]));

[-1.00000000000000, -0.497870683678640e-1, -0.247875217666636e-2, -0.123409804086680e-3, -0.614421235332820e-5, -0.305902320501826e-6, -0.152299797447126e-7, -0.758256042791190e-9, -0.377513454427910e-10]

It agrees with the table!

 

Try Table 3, column 1:

seq(eval(diff(f(eta),eta$2), Solve(K=1, Pr=0.733, sigma=5, c= 1, M= 3, n= 0.5, lambda= L)(0)), L= 2.0..4.0, 0.5);

HFloat(2.6198066501250747), HFloat(2.4131979709211966), HFloat(2.216852164633253), HFloat(2.040203526244527), HFloat(1.8844737238149782)

 

Totally wrong. It's even moving in the wrong direction!


 

Download nanofluidBVP.mw

 

First 308 309 310 311 312 313 314 Last Page 310 of 709