sand15

1628 Reputation

16 Badges

11 years, 187 days

MaplePrimes Activity


These are replies submitted by sand15

@nm 

which avoids constructing  H_lines and V_lines:

restart
max_X:=8: 
max_Y:=5:

grid_points:= [seq(seq([n,m],n=0..max_X),m=0..max_Y)]:
the_points := plots:-pointplot(grid_points, symbol=solidcircle, symbolsize=12, color=red):
plots:-display(
  the_points
  , title="my grid"
  , axis[1] = [gridlines = [max_X, color = blue]]
  , axis[2] = [gridlines = [max_Y, color = blue]]
  , scaling=constrained
);

@Ali Guzel  

restart;

with(plots):

T:=60:
Digits:=10:
epsilon:=20:

F:=(x,y)->epsilon*(1-x^2)*y-x:

sys:=diff(X(t),t)=Y(t),diff(Y(t),t)=F(X(t),Y(t)):

vars:={X(t),Y(t)}: ic:=X(0)=0.2,Y(0)=0.1:

pow := [$-4..-2]:

for p in pow do
  sol[p] := dsolve(
              [sys,ic]
              , numeric
              , stepsize=10^p
              , method=classical[foreuler]
              , output=listprocedure
              , maxfun=0
            ):
end do:

col := [red, green, blue]:

display(
  seq(
     odeplot(
       sol[p]
       , [t,X(t)]
       , 0..T
       , axes=normal
       , style=line
       , numpoints=n
       , labels=["t","x"]
       , color=col[abs(p)-1]
       , legend=typeset(h = cat(`#msup(mo("10"),mo("`, p, `"))`))
     )
     , p in pow
  )
)

 

default_sol := dsolve(
                 [sys,ic]
                 , numeric
                 , output=listprocedure
                 , maxfun=0
               ):

display(
  odeplot(
     sol[pow[1]]
     , [t,X(t)]
     , 0..T
     , axes=normal
     , style=line
     , numpoints=n
     , labels=["t","x"]
     , color=col[abs(pow[1])-1]
     , legend=typeset(h = cat(`#msup(mo("10"),mo("`, pow[1], `"))`))
   )
   ,
   odeplot(
     default_sol
     , [t,X(t)]
     , 0..T
     , axes=normal
     , style=line
     , labels=["t","x"]
     , color=gold, thickness=7, transparency=0.7
     , legend="default solver"
   )
  )

 
 

 

Download Are_you_aware_of_convergence_issues.mw

@Alfred_F 

Years ago I bought the first French translation (2015) of this famous book "What is mathematics? An elemenrtary approach to ideas and methods" by Richard Courant and Herbert Robbins (1941).

Its 1st chapter concerns the natural numbers and on line 9 of the the introductory paragraph (page 7) it is written nombres naturels1✻ ... (The footnote  says that the notes tagged with an asterisk are from the book's translator and are thus absent in the vernacular version of the book).
This note (1) can be read page 27 and I reproduce it here (I hope I'm not infringing copyright in doing so).
Note_1.pdf
There is also a discussion about positive, strictly positive and non negative numbers.
 

@Mac Dude 

About the density plot: you could of use the orientation option in the Histogram3D plot command to get this density plot.
But you would have to adjust the lightning model and remove cell boundaries to get a pretty display.
I believed it's better to have a specific procedure dedicated to density plotting.

@emendes 

If it's not too much to ask, could you explain to me what your code is meant to calculate?
(This is more of a curiosity than something that could help me to answer your question.)


 

restart

ic := 1, 2, 3:
rec := rsolve({f(n) = f(n-1) + f(n-2) + f(n-3), seq(f(i)=ic[i], i=1..3)}, {f});

{f(n) = Sum((_R+1)*(1/_R)^n/(3*_R^2+2*_R+1), _R = RootOf(_Z^3+_Z^2+_Z-1))}

(1)

F := unapply(value(rhs(rec[1])), n)

proc (n) options operator, arrow; sum((_R+1)*(1/_R)^n/(3*_R^2+2*_R+1), _R = RootOf(_Z^3+_Z^2+_Z-1)) end proc

(2)

print~([ seq([i, ic[i]], i=1..3), seq([n, F(n)], n=4..10) ]):

[1, 1]

 

[2, 2]

 

[3, 3]

 

[4, 6]

 

[5, 11]

 

[6, 20]

 

[7, 37]

 

[8, 68]

 

[9, 125]

 

[10, 230]

(3)

ic := 'ic';  # unassign the name ic

ic

(4)

tri := proc(ic)
  local rec, F:
  rec := rsolve({f(n) = f(n-1) + f(n-2) + f(n-3), seq(f(i)=ic[i], i=1..3)}, {f});
  F   := unapply(value(rhs(rec[1])), n);
end proc:

ic := [u, v, w];
tri(ic)(n)

[u, v, w]

 

sum(-(-2*_R^2*v+_R^2*w-2*_R*u-_R*v+_R*w+u+v-w)*(1/_R)^n/((3*_R^2+2*_R+1)*_R), _R = RootOf(_Z^3+_Z^2+_Z-1))

(5)

ic := [1, 2, 3];
tri(ic)(n);
map(n -> tri(ic)(n), [$1..10])

[1, 2, 3]

 

sum((_R+1)*(1/_R)^n/(3*_R^2+2*_R+1), _R = RootOf(_Z^3+_Z^2+_Z-1))

 

[1, 2, 3, 6, 11, 20, 37, 68, 125, 230]

(6)

ic := [1$3];
tri(ic)(n);
map(n -> tri(ic)(n), [$1..10])

[1, 1, 1]

 

sum(-(-_R^2-2*_R+1)*(1/_R)^n/((3*_R^2+2*_R+1)*_R), _R = RootOf(_Z^3+_Z^2+_Z-1))

 

[1, 1, 1, 3, 5, 9, 17, 31, 57, 105]

(7)

 


 

Download tripatouille.mw

It would be much simpler to help you if you upload your worksheet (use the big green up-arrow in the menubar, watchout: the file name must not contain special characters).

An advice: verify if the commands you use are all threadsafe (for instance zip is not)

@dharr 

You're right, I have been careless in the conclusion I drew from Descarte's rule of signs.

I got back to what I did to see that in fact I don't use it explicitely (I don't even know why I mentioned it: probably I thought by the time it would help me, but it didn't).
Anyway, I found my error.
At some point I solve a 3r degree equation and Maple, as usual, gives me the roots r1, r2, r3 where r2 and r3 appears in conjugate complex forms, let's say r2 = a+Ib and r3 = a-Ib.
Given these 3 rootsshould be real be real I conclude that b = 0 and thus r2 = r3, so there is a root of multiplicity 2 (for instance y = z).

Replacing z by y in the three symetric functions of (x, y, z) led me to w = 0 and x = 0.
I vaw to quick: in the expression r2 = a+Ib nothing tels me that 'a', or 'b', or both are not imaginary themselves, meaning the shotcut "b must be null" is an unforgiveable error.

I should have taken the time to analyze my work instead of assuming I was on the right track.
As they say, forewarned is forearmed; let’s hope this serves me well as a lesson.

See you next time!

@dharr Interesting...

From wiki "entiers naturels" (which translates into "natural integer", eventually English is not that difficult) you read

whose translation is "Thus, in textbooks that include 0 among the natural numbers—which is the case in most school textbooks in France—we find the notations ...".

I was always taught 0 was the smallest natural integer and never realized that other cultures could suse other conventions.
As @Alfred_F seems to be French-speaking (maybe French) I implicitely understood 0 was a natural number... which is why I was quite surprised by the "big" solution he gives at the end of its post.

Never mind:  I "demonstrated" in my reply that there was no other solution but 'w = 0'... which is obviously wrong because this "big" solution is indeed a solution.
This prove I did some mistake in my reasoning. I will try to track it if I have some time.

Isn't {u = 2, v = 2, w = 0, x = 0, y = 2, z = 2} the smallest solution (up to an arbitrary permutation on x, y, z)? 

My point of view:

  • The set E := {x*y*z = w^2, x+y+z = u^2, x*y+x*z+y*z = v^2} of conditions can be interpreted this way:
    • Let Z(t) = t^3 + b*t^2 + c*t + d and let (x, y, z) the rootsof the equation Z(t) = 0.
    • The symetric functions of (x, y, z) are
      x+y+z = -b
      x*y+y*z+z*x = c
      x*y*z = -d
    • Thus the equation Z(t) = 0 writes 
      t^3 - u^2*t^2 + v^2*t - w^2 = 0
  • Note that by Descarte's rule of signs this 3rd degree polynom with indeterminate 't' has 3 sign chsnges and thus 3 positive roots (x, y, z)
     
  • After some algebra one can show (if i did no mistake) that one root, let's say 'x' has multiplicity 1 and that y = z, meaning the second root (let's say 'y') has multiplicity 2.
    This means Z(t)  = 0 writes
    (t-x)*(t-y)^2 = t^3 - u^2*t^2 + v^2*t - w^2
  • Let EQ this equation.
    Here are all the (u, v, w, x, y, z} solutions we get
    sol := { solve(identity(EQ, t), {x, y, u, v, w}) }:
    print~(sol):
  • Let us take w = 0 (∎).
    We get
  • The first solution is basically the trivial one you got (take u = 0)/
    Among the four other ones, only one contains no minus sign and is thus the "generator" for other solutions based on the choice w = 0.
    It is  easy to see that y = 2 (and thus  z = 2) gives the smallest, non identically null, solution I presented at the top of this reply.
    All other solutions are given by y = 2*p^2 where p is any positive number.

(∎) One might think the choice  w = 0 is arbitrary. But for other integer positive values of 'w', we need that the system

sys := {2*y^3+w^2 = p^2, y*(y^3+2*w^2) = q^2}

have integer solutions (p::posint, q::posint) in order that positive integer (x, y, u, v, w) could exist.
Multiplying sys[1] by 'y' and substractind twice sys[2] gives

-3*w^2*y = p^2*y-2*q^2

whose the only integer solution is


As you see we return to the 'w = 0' initial choice and noother value of 'w' can lead to a solution.

@acer 

Thanks for all these details.
Funny to see that you were able to see "that the conditional block of the original code... the inert representation of that original procedure ..". This doesn't seem that obvious to me at first sight, but I guess I will undoubtedly always be a bit of a DIYer when it comes to using Maple.

And to end this reply, yes, you're right writing "that if this kind of manipulation gets too involved then one might as well print the proc, copy, and edit the new version as desired".
I already did that sometimes and, after all, sparsematrixplot is quite easy rewrite or customize.

By the way, my own question was motivated by this question from  @Mac Dude.

@acer 

Great, and thanks a lot.

I use to use the gap option with matrixplot and I was disappointred not to find it for sparsematrixplot.
Your first proposal is quite “technical,” but very effective, because, but if I understand correctly, it's a fairly general method for overriding the default options of any (?) tracing procedure.

Your second one is based on the use of the sparse option in matrixplot and looks more like a trick, but a very astute one (and indeed Maple 2015 does not support this option, so it's a dead line for me).

Given I'm stuck with Maple 2015 right now I will then use your first approach.

A question: how did you end writting this

`plots/sparsematrixplot` := FromInert(
    subsop(
      [5,3]=[
              op([5,3], ToInert(eval(`plots/sparsematrixplot`)))
              ,_Inert_STATSEQ(_Inert_ASSIGN(_Inert_LOCAL(8), _Inert_FLOAT(_Inert_INTPOS(0), _Inert_INTPOS(0))))
            ][]
            ,ToInert(eval(`plots/sparsematrixplot`))
   )
):

?
This is all but trivial to find and writting

showstat(`plots/sparsematrixplot`)

didn't help me to understand your  intellectual journey.

Thanks again

@Mac Dude 

Meanwhile here is slight variant of my Histogram3D procedure described in my previous reply (usable for both dependent and independent random variables).

I also included a new feature to draw the equivalent of what plots:-densityplot does for functions. The code is very simple but strongly benefits from @acer's help, see HERE.

NOTE: In my code (I forgot to mention it uses Maple 2015) I used the first method acer proposed (in the code below his piece of code is reproduced in green highlighted b
lue font: it should be easy for you to locate it and replace it by one of the other solutions acer proposed, depending on your Maple version).

3DHist_3DDensity_2.mw



As I told you before one can get far better rendering for both the histogram and, in case this would interest you, the density plot.
But doing so requires some time to develop a proper code based upon PLOT/POLYGONS and/or plottools built-in functions. The plots:-matrixplot I used is a quick shortcut to drax the histogram you have in mind, but it is versatile enough to produce pretty outputs.
The same observationholds for plots:-sparsematrixplot as a workaround to a "discrete" plots:-densityplot.

To jump on your initial remark "I could of course program one myself, but as I am lazy and this seems like a common kind of plot (at least in statistics and in physics), someone may well have done such a thing. DDG search did not find anything, at least not wrt. Maple."
I'm sure you can and I understand your reluctance to do it.
And yes, it's a very common plot in Statistics and that's unfortunate Maple doesn't offer any feature to draw 3D histograms and density plotof a cloud of points (as a lot of other languages do).


This is the kind of result you can obtain using an Histogram3D procedure based upon plottools:
3DHist_using_plottools.mw

@Mac Dude

Code: 3DHist_dependent.mw

Scatterplot of a sample drawn from a couple ofdependent random variables:
 
3D histogram of this same sample

Watchout: I assumed xvec and yvec are samples of two independent random variables X and Y

In my initial answer the procedure Histogram3D might produce an error when xtally or ytally were sparse (meaning for instance that some indices did not existed in one or two of these tables).
The new Histogram3D procedure in worksheet 3DHist_robustified.mw accounts for this situation (see text in [courier black] font).

Example 
             3D histogram over a 30-by-30 grid                                     sparse structure of the 30-by-30 grid

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