Carl Love

Carl Love

26827 Reputation

25 Badges

11 years, 278 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Like this:

seq(cat("%15s" $ k), k= 1..4);

To continue, you need to say what mathematically makes the difference between solid and dashed lines.

params:= {
    h= piecewise(
        z <= d+1,   1,
        z <= d+4,   1-delta/2*(1 + cos(2*Pi*(z - 1 - 1/2))),
        z <= d+6,   1 
    ),
    w0= -c*h^2/4 + 3/64*(b*c-4*a)*h^4 + 19/2304*b*(b-4*a)*h^6,
    w1= c/4 + (4*a-b*c)*h^2/16,
    w2= (4*(b*c-4*a)-b*h^2)/256,
    w3= b*(b-4*a)/2304,
    a= x4*S*Gr*sin(alpha)/(4*x1*x5),
    b= 1/Da + x3*M/(x1*(1+m^2)),
    c= Dp/x1,
    Dp= 96*x1/((6-b*h^2)*h^4)*(F + a*h/24 - 11/6144*b*(b-4*a)*h^8),
    x1= ((1-phi1)*(1-phi2))^(-5/2),
    x2= (1-phi2)*(1 - phi1 + phi1*Rs1/Rf) + phi2*Rs2/Rf,
    x3= shnf/sf,
    x4= (1-phi2)*(1 - phi1 + phi1*RBs1/RBf) + phi2*RBs2/RBf,
    x5= khnf/kf,
    shnf= sbf*(ss2+2*sbf-2*phi2*(sbf-ss2))/(ss2+2*sbf+phi2*(sbf-ss2)),
    sbf= sf*(ss1+2*sf-2*phi1*(sf-ss1))/(ss1+2*sf+phi1*(sf-ss1)),
    khnf= kbf*(ks2+2*kbf-2*phi2*(kbf-ks2))/(ks2+2*kbf+phi2*(kbf-ks2)),
    kbf= kf*(ks1+2*kf-2*phi1*(kf-ks1))/(ks1+2*kf+phi1*(kf-ks1)),
    RBs1= 8933*16.7e6,
    RBf= 1063*1.8e6,
    RBs2= 6320*18e6,
    kf= 0.492,
    sbf= 6.67e-1, ss2= 2.7e-8,
    sf= 6.67e-1, ss1= 59.6e6,
    ks2= 76.5, kf= 0.492, ks1= 401,
    phi1= 0.01, phi2= 0.02, alpha= Pi/4, m= 0.5, Da= 0.1, Gr= 5, 
    delta= 1, S= 0.5,  d= 1,
    F= 1, z= 1
}:
                                      
W1:= eval[recurse](w0+w1*r^2+w2*r^4+w3*r^6, params):

plot(
    eval~(W1, M=~ [2, 5, 7]), r= 0..1,
    color= [red, blue, green], axes= boxed,
    legend= typeset~(M^2=~ [2,5,7]),
    labels= [r, w(r,1)], labeldirections= [horizontal, vertical],
    axesfont= [helvetica, 16], labelfont= [helvetica, 16],
    caption= typeset(
        ``(seq(v = eval[recurse](v, params), v= (m, Da, Gr, S, alpha, phi1, phi2)))
    )
);

If I translate your f into 1D-input, it becomes this:

f:= (u[1], u[2])-> u[1]^2*diff(u[2], x);

It's simply a nuance of Maple syntax that indexed variables such as u[1] are not allowed to be formal parameters of procedures. There is a simple fix that allows the variables to still display as subscripted:

f:= (u__1, u__2)-> u__1^2*diff(u__2, x);

c__1 is an alias. Aliases are very different from ordinary variables (and I think it was a bad design choice to use them for dsolve). One of the weird things about aliases is that if they are used in a procedure, they must be defined before the procedure is defined, whereas an ordinary global variable only needs to be defined before the procedure is called. To verify this, re-enter the definition of foo (without restart), and re-call foo (shown below). Also shown below is a safer way to code foo by being agnostic about what the constant will be.

Your Example 2 Question: Solving for _C1 works, even though the ode has c__1 , why?

Because _C1 is an ordinary global variable, not an alias. c__1 is an alias for _C1. To see this:

_C1;
                       c__1

Your Example 3 Question: Forcing arbitraryconstants = subscripted it still does not work inside proc. Why?

c__1 is still an alias. The weird interaction between aliases and procedures supercedes everything else.

 

Example (1) (modified) solving for constant of integration fails inside proc but works outside

 

restart;

foo:=proc(ode::`=`)
local sol,the_constant;
   sol:=dsolve(ode);
   print("sol is ",sol);
   the_constant:=solve(sol,c__1);
   print("the constant is ",the_constant);
end proc
:

#this does not work
ode:=diff(y(x),x) = 3/4*y(x)/x;
foo(ode)

diff(y(x), x) = (3/4)*y(x)/x

"sol is ", y(x) = c__1*x^(3/4)

"the constant is "

#this works
ode:=diff(y(x),x) = 3/4*y(x)/x;
sol:=dsolve(ode);
print("sol is ",sol);
the_constant:=solve(sol,c__1);

diff(y(x), x) = (3/4)*y(x)/x

y(x) = c__1*x^(3/4)

"sol is ", y(x) = c__1*x^(3/4)

y(x)/x^(3/4)

 

####
# All below added by Carl
####
#
#This works (Note: No additional restart!):

foo:= proc(ode::`=`)
local sol,the_constant;
   sol:=dsolve(ode);
   print("sol is ",sol);
   the_constant:=solve(sol,c__1);
   print("the constant is ",the_constant);
end proc
:

foo(ode);

"sol is ", y(x) = c__1*x^(3/4)

"the constant is ", y(x)/x^(3/4)

restart:

#A better way to code foo:
foo:= proc(ode::`=`)
local
    sol:= dsolve(ode),
    C:= indets(sol, And(symbol, Not(constant))) minus indets(ode, And(symbol, Not(constant))),
    the_constant
;
    print("sol is ",sol);
    if nops(C) <> 1 then error "Can't find specific constant:", C fi;
    the_constant:= solve(sol, C[]);
    print("the constant is ", the_constant)
end proc:

ode:=diff(y(x),x) = 3/4*y(x)/x;
foo(ode);

diff(y(x), x) = (3/4)*y(x)/x

"sol is ", y(x) = c__1*x^(3/4)

"the constant is ", y(x)/x^(3/4)

Download AliasedDsolveConstants.mw

When applied to a first-order nonlinear ODE, the term homogeneous has a completely different and standard definition than the one used for linear ODEs and referenced by @C_R . There is no connection between the two definitions. A quick Google search turns up many references to this definition of homogeneous for 1st-order nonlinear ODEs, and it's equivalent to the defintion given on the Maple help page that you referenced. So, I think that that quote from Advanced Engineering Mathematics with Maple is not telling the whole story.

Each of your three algebraic solutions of your ode for diff(y(x), x) can be put in the form F(y(x)/x) (using three slightly different Fs). To see this, after your PDEtools:-Solve command do

simplify(eval(rhs~([sol]), y(x)= x*v), symbolic);

You'll see that none of the three depend on x. (Note that v = y(x)/x.)

What benefit do you think you're getting by using alias? Aliases are quite weird, and if you're thinking of XX as just another variable, you'll likely make mistakes. Why not instead use

export XX:= Sqr1;

or

macro(XX= Sqr1);

?

It seems to me as if you just want an abbreviation for a longer command name. The export above will do that. The only possible reason I know of to use an alias is that it's an abbreviation that gets back-substituted on output. 

Whenever you use is, you need to consider the possibility that it returns FAIL instead of true or false. If all of the is-es return FAIL or false, then ormap will return FAIL.

foo:= proc(a, {S::list:= []})
    S<>[] and
    S[1]::{1, And(specfunc('typeset'), identical(true) &under (s-> ormap(is, s, 1)))}
end proc:

 

1. First, let's focus on setting up your equations; solving them is relatively unimportant at this point. Solving them symbolically may not be possible, but it also may not be necessary or even useful at all. The final mouse driver will use a compiled numeric solver dedicated to this particular set of equations. Solutions accurate to 3 or 4 significant digits will be sufficient. Maple can be used to produce such compiled code, which will ultimately run independently of Maple (not even requiring the existence of Maple).

From the code that you showed, I'm not convinced that your equations are set up correctly. In my plain human reading of your code, I've noticed several potential problems (some mentioned below). I cannot really check your code until you upload a worksheet.

2. You can upload a Maple worksheet using the green up-arrow on the toolbar of the MaplePrimes editor. If you get a message about the upload having failed, ignore it. That message only means that MaplePrimes cannot display the uploaded worksheet. Almost always the file link is correctly uploaded anyway.

3. Forget about your linearization T; only work with L. An "exact" symbolic solution of a crude symbolic approximation is of very dubious value.

4. This is a description of the main potential error that I see in your code. I will call a variable indexed if it ends with a pair of square brackets ( [...). (Whether there's anything in the square brackets is irrelevant.) So, in your code, p[ni]m[m]s[n]phi[1], ..., phi[3]theta[1], ..., theta[6], etc., are all indexed variables. It is highly error prone to use variables in both indexed and unindxexed forms in the same problem. For example, you use unindexed theta as a variable (it's one of the solution variables) and you assign values to the indexed forms theta[1], etc.

5. I will call a variable subscripted if it contains a double underscore in its name. For example, your p__n is subscripted. Subscripted variables are immune to any of the potential errors that indexed variables can lead to. Note well that I'm not saying that there's anything wrong with using indexed variables; they can be quite useful. The problems may occur when you try to use them simultaneoulsy in both indexed and unindexed form, and especially when some of those have values assigned to them.

Assuming that cornerPoints is an array or matrix with dim = 3 rows and 4 columns, you could also do this, which requires no indexing or dimension references:

plots:-textplot3d(
    [op,`<,>`]~(convert(cornerPoints^%T, listlist)), 
    align= {above, right}, font= [Courier, bold, 20]
);

The absolute value of your final computation is independent of the order of the 4 vectors, so it doesn't matter which one you subtract from the others.

V:= (p,q,r,s)-> abs(LinearAlgebra:-`&x`(p-s, q-s) . (r-s));

proc (p, q, r, s) options operator, arrow; abs(LinearAlgebra:-`&x`(p-s, q-s).(r-s)) end proc

Pts:= ['LinearAlgebra:-RandomVector(3)' $ 4];

[Vector(3, {(1) = -59, (2) = 12, (3) = -62}), Vector(3, {(1) = 72, (2) = 42, (3) = 18}), Vector(3, {(1) = 52, (2) = -13, (3) = 82}), Vector(3, {(1) = 70, (2) = -32, (3) = -1})]

for pts in combinat:-permute(Pts) do V(pts[]) od;

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

851671

 

Download BoxProduct.mw

Your final formula can be simply proved by mathematical induction, without any need for a closed-form algebraic expression for the Fibonacci sequence.

restart:

interface(prompt= ""):

The prroposition to prove:

eval(G(2*n) = G(2*n+1) + G(2*n+2), G= arctan@`/`@F);

arctan(1/F(2*n)) = arctan(1/F(2*n+1))+arctan(1/F(2*n+2))

for positive integers n where F = combinat:-fibonacci.

 

We start by having Maple confirm a fairly well-known identity for a sum of arctangents:

combine(arctan(1/a) + arctan(1/b)) assuming a > 1, b >= 1;

arctan((b+a)/(a*b-1))

So we just need to verify this Fibonacci identity:

1/F(2*n) = subs({a= F(2*n+1), b= F(2*n+2)}, op(%));

1/F(2*n) = (F(2*n+2)+F(2*n+1))/(F(2*n+1)*F(2*n+2)-1)

It'll be easier to work with in polynomial form. This will be the verification goal of mathematical induction:

Goal:= numer((lhs-rhs)(%)) = 0;

-F(2*n)*F(2*n+2)-F(2*n)*F(2*n+1)+F(2*n+1)*F(2*n+2)-1 = 0

To form the induction hypothesis, change n to n-1:

IH:= eval(Goal, n= n-1);

-F(2*n-2)*F(2*n)-F(2*n-2)*F(2*n-1)+F(2*n-1)*F(2*n)-1 = 0

Verify the induction base case:

is(eval(Goal, {n=1, F= combinat:-fibonacci}));

true

And verify the Goal by simplifying with respect to the induction hypothesis and the defining recursive formula for the Fibonacci sequence:

is(simplify(Goal, {IH, seq(F(2*n+k)+F(2*n+k-1) = F(2*n+k+1), k= -1..1)}));

true

 

Download FibonacciAtanIdentity.mw

@C_R Your understanding of (...)(...) is totally correct. The reason that the extra parentheses (as shown by @nm) are necessary is the precedence of Maple's binary infix operators: -> has lower precedence than = (just like has lower precedence than *). See ?operators,precedence.

You can avoid the need for the extra parentheses by replacing x->x with x@@0, which is the identity function. `@@` has higher precedence than =. You could also use 1@@0_@@0, etc.

The multiplication symbol can be easily suppressed:

make_nice:= e-> subs(
    "&lowast;"= "&InvisibleTimes;",   
    InertForm:-Display(evalindets(e, And(`*`, fraction &under curry(op,1)), `%*`@op))
): 

Also, by changing @nm's typespec from `&*`(fraction, anything) to And(`*`, fraction &under curry(op,1)), my procedure can handle products of any number of factors, such as 3/8*ln(x/y)*exp(z).

You just need to add the arrows option to your DEplot3d command. For example, by adding arrows= cheap I get this:

In your display command, remove the options gridlines= true and axes= none, and add the options

axis[1]= [gridlines=17, tickmarks=0], axis[2]= [gridlines= 8, tickmarks=0], view= [0..16, 0..7]

You should get this:

1 2 3 4 5 6 7 Last Page 1 of 385