3728 Reputation

17 Badges

6 years, 21 days

MaplePrimes Activity

These are replies submitted by mmcdara


Just a remark about your procedure EQSYS 

  • it returns an error when the two planes are parallel and unconfounded
  • it responds "false" if the planes are confounded and the line belongs to each of them
    (assuming that the intersection of two confounded planes is these same planes and that each line in these latter can thus be considered as a particular intersection among many; thus "true" seems to be a better answer, even if geom3d syays that the planes do not intersect ... maybe a definition of the intersection of two planes?)
sys1 :={16*x-2*y-11*z=0, 14*x-y-10*z=3}:
sys1a:={16*x-2*y-11*z=0, 32*x-4*y-22*z=0}: #confounded planes
sys1b:={16*x-2*y-11*z=0, 32*x-4*y-22*z=1}: #parallel uncounfounded planes 

sys2:={(x-2)/3=(y-5)/2, (y-5)/2=(z-2)/4}:

EQSYS(sys1 , sys2);  # answer is true
EQSYS(sys1a, sys2);  # answer is false but the line is in the plane
EQSYS(sys1b, sys2);  # error

In cases sys1a and sys1b geom3d would respond

plane(P1, ...., [x,y,z]);
plane(P2, ...., [x,y,z]);
line(L1, [P1,P2]);

P1 = sys1a --> intersection: the two planes P1 and P2 are the same
               Error, (in geom3d:-line) the two given planes are parallel

P2 = sys1b --> intersection: the two planes P1 and P2 are parallel
               intersection: the given objects do not intersect



Of course...
But in a general situation you must distinguish several cases for the two planes (which may intersect or not):

  • Let a*x+b*y+c*z+d the equation of a plane and a'*x+b'*y+cC'*z+d' the equation of the second one
  • if a'/a=b'/b=c'/c=d'/d then P and P' are identical
  • if a'/aA=b'/b=c'/c then Pand P' are parallel
  • if at least one of the three equalities a'/a=b'/b, b'/b=c'/c, a'/a=c'/c  is not true then P and P' intersect along a line D
    In this case the elementary way to obtain the solution is:
    • to consider that one of the unknown (x or y or z) is a "parameter" t (for instance x=t) and solve a system of two equations in the two remaining unknowns.
    • to solve a system of two equations in the two remaining unknowns, for instance 
      solve({a*tx+b*y+c*z+d=0,  a'*t+b'*y+cC'*z+d'}, {y, z})


For your case this gives (for instance)

e1 := 16*s-2*y-11*z:
e2 := 14*s-y-10*z-3:

# Set x as a parameter "s"
E0 := [x=s, solve({e1, e2}, {y, z})[]]
               [           2     11      4     2]
               [x = s, y = - s + --, z = - s - -]
               [           3     3       3     3]

# verify that this line is the one geom3d finds
plane(P1, 16*x-2*y-11*z, [x,y,z]):
plane(P2, 14*x-y-10*z-3, [x,y,z]):
line(L1, [P1,P2]):
E1 := [x, y, z] =~ Equation(L1,'t');
              [    1                             ]
              [x = - + 9 t, y = 4 + 6 t, z = 12 t]
              [    2                             ]

# as the parameterization can be different, express s as a finction of t
# and use this to re-parameterize E0 in terms of E1 to verify that E0 and 
# E1 are indeed identical

solve(rhs(E0[1])=rhs(E1[1]), s);
eval(E0, s=%);
                            - + 9 t
              [    1                             ]
              [x = - + 9 t, y = 4 + 6 t, z = 12 t]
              [    2                             ]

So to answer your question "Is it possible showing without geom3d ?", just do 

E0 := [x=s, solve({e1, e2}, {y, z})[]]

And now to verify if your two lines are the same :

X := s:
Y := solve((y-5)/2=(s-2)/3, y):
Z := solve((z-2)/4=(s-2)/3, z):
[x, y, z] =~ [X, Y, Z]

You might think it's shorter than using geom3d, but for general planes this doesn't work (look at the different situations in the beginning of this reply).
Why did I use geom3d? Just because I didn't know if your two planes where confounded, parallel or secant.

      It contained an error 

@Preben Alsholm  @Michael_Watson

Right, but I think the original mistake could be a typo:

ex := L/(z+sqrt(z^2 + L^2));  # + or - instead of *
eval(ex, L=alpha*z);
simplify(%) assuming z > 0
                     /     2    \         
                     \alpha  + 1/      + 1


eval(ex, z=beta*L):
simplify(%) assuming L > 0
                    /    2    \            
                    \beta  + 1/      + beta



I naively tried to use rsolve for the first problem but I got no answer. 
I am wrong writting this

rsolve({a(n+1)=M(a(n), b(n)), b(n+1)=a(n)*b(n)/M(a(n), b(n))}, {a, b})

or am I facing some limitation of rsolve?

Concerning problem b:

  • I have been stucked for a while trying to find an invariant or some map from (a, b) to (u, v) whoch could transform the problem in something simpler.
    Do you have some more information about the proof that "The limit M⊗N is the logarithmic mean"?
  • About the convergence speed :
    cat(`b__`, n+1)-cat(`a__`, n+1) = N(a__n,b__n)-M(a__n,b__n);
                                     1        1     
                   b__n+1 - a__n+1 = - b__n - - a__n
                                     2        2     
     cat(`b__`, n+1)-cat(`a__`, n+1) = (b__0-a__0)/2^(n+1)
                                       b__0 - a__0
                     b__n+1 - a__n+1 = -----------
                                         (n + 1)  


Point 1
I had interpreted your claim "Using Maple 2021.1, we can work around that by using..." as a way (not the way) to handle the problem in Maple 202.
Does your last reply mean that 

plot('funktion(r)', r=r1 .. r2);

won't work in Maple 2021 in the case we use Units?

Point 2
There is a basic problem in the definition of the funktion the OP uses

funktion := x -> if 2*Unit('m') < x then 1/x; else 2*Unit('m'); end if:

For the comparison to be possible x must ne e Length (m, cm, inches, ...)
Assuming this is the case the first part of the test (2*Unit('m') < x)  returns a quantity homogeneous to 1/Length, as the else part of the test funktion returns a quantity homogeneous to a Length.
Think to this: what label would you put on the vertical axis if the plot was to be from 1.5m to 4m?
Would it be m or 1/m ?

Point 3
As you probably remember I use Maple 2015
Could you please tell me if newer versions give two identical plots or if, as with Maple 2015, the second plot is incorrect.

funktion := x -> piecewise(2*op(2, x) < x, 1/x, 1/(2*op(2, x))):
a := 1.0:
b := 4.0:

correct := plot(
  [seq([r, op(1, funktion(r*Unit('m')))], r in [seq(a..b, (b-a)/100)])],
  labels=[typeset(r*[Unit('m')]), typeset(f*[Unit('m^(-1)')])],
  color=blue, thickness=7, transparency=0.5, legend="correct"

wrong := plot(funktion, a*Unit('m')..b*Unit('m'), color=red, legend="wrong"):

plots:-display(correct, wrong)

Said otherwise, what plot do you obtain with Maple 2021 if you run your code with this slight modification?

funktion := x -> if 2*Unit('m') < x then 1/x; else 2*Unit('m'); end if:
r1 := 1.0*Unit('m'):
r2 := 6*Unit('m'):

plot(funktion, r1 .. r2);

I hope you get this

Thanks in advance

@Thomas Richard 

Is there any inconvenience to use 

plot('funktion(r)', r=r1 .. r2);

instead of

plot(funktion, r1 .. r2);

The fact that we have to change the syntax of the plot procedure seems to me to be a weakness of Maple, a lack of genericity of this procedure in the sense that Maple should (IMO) understand that it does not have to evaluate the function prematurely.
This is however what it does here

plot([seq([r, funktion(r)], r=3..6)])

It is funny how the same people can sometimes understand the meaning of the word "extract" while categorically refusing it at other times


Tgank you acer


To understand without any procedure why A->B->D->E->C->A (or any cyclic permutation of it) is the shortest cycle:

  • introduce these 3 fictitious points P=(0, 5), Q=(5, 5), R=(5, 0)
  • the cycle A->P->C->Q->E->R->A has length 20 (obvious)
  • if you decide once arrived at E to take the interiorreturn path E->D->B->A instead of E->R->A you won't reduce the length of the path from E to A.
  • thus the shortest cycle (a shortest cycle among many others) has length 20

@acer @tomleslie

I agree that "The plots:-surfdata command can produce either a 3D surface plot or a 2D density-style plot", but if one does this

cosdata := [seq([ seq([i,j, evalf(cos((i+j)/5))], i=-5..5)], j=-5..5)]:
a := surfdata(cosdata, dimension=2, axes=none):

one sees that op(a) is a MESH structure.
In the plot,structure help page its written
The MESH structure represents surfaces in 3-D space defined by a grid of values

Thus, I would be tempted to say that although surfdata can produce a 2D density graph, it seems essentially to be a 3D graph and that dimension=2 "simply makes a projection" of the 3D surface onto the [z,y] plane?

So, is Tom's statement really that wrong? 



It was not a suggestion I made to the OP but to you, you were free to take it into account or not. 
If I had wanted to write to the OP I would have done so.


"By the way it is alway helpfull..."
It was just a suggestion for a tiny improvement of your own code, more a private communication than a real answer.
So I felt that publishing a code 99% identical to yours would be of no use, not to mention plagiarism.
And between you and me, I don't think you spent that much time putting this suggestion into your code.


For the phase portrait you can try this
(3 types of solutions are displayed if you play the animation)


  {diff(x(t), t)-v(t)=0, diff(v(t), t)-2*Heaviside(x(t)^2-1)=0},
  [x(t), v(t)],
  [[v(0)=-2, x(0)=-4], [v(0)=-1, x(0)=2], [v(0)=1, x(0)=-2]]
  ,linecolor=[gold, green, cyan]
  ,arrows=curve        # try comet also
  #,dirfield=[30, 30]  # try other values


Why not define Cnames this way;

Vnames := map(u -> sprintf("%a", u), faces);

then the names of the faces directly refer to the vertices of the primal graph.


An event-based solution (same system as in my previous reply)

sol := dsolve(
           [[x(t)= 1, And(v(t) > 0)], A(t)=1], 
           [[x(t)= 1, And(v(t) < 0)], A(t)=0], 
           [[x(t)=-1, And(v(t) > 0)], A(t)=0], 
           [[x(t)=-1, And(v(t) < 0)], A(t)=1]

# to verify A(t) behaves correctly plot A(t) versus x(t)

plots:-odeplot(sol, [x(t), A(t)], t=0..2)

An Heaviside-based solution

ode := diff(x(t), t)=v(t), diff(v(t), t) - 2*Heaviside(x(t)^2-1) = 0:
ic  := x(0)=-2, v(0)=1:                  
sys := {ode, ic}:
sol := dsolve(sys, numeric);

PS: I guess someone here will be able to present a way to solve formally this "system" (I can't).


First 10 11 12 13 14 15 16 Last Page 12 of 86