mmcdara

7891 Reputation

22 Badges

9 years, 60 days

MaplePrimes Activity


These are replies submitted by mmcdara

@tomleslie 

Thanks you for this detailed contribution.


You wrote "First time through the "offending" loop the 'sol' procedure will been evaluated from t=0 to t=1, so the event at t=0.5 is available. Second time throught the loop the 'sol' procedure has been evaluated from t=0 to t=0.5, so the event at t=0.3 is available - and so on" : this was also my explanation (and yes, I agree, "when an even doesn't fire! I would have expected some kind of "warning" or NULL answer, rather than returning t=0, x(t)=0" (this situation is always worrying).

I appreciate this explanation of yours
# Technically this is not an "initialization".
# In order to obtain the result at t=1, Maple
# has to "run" the solution module "sol" from
# t=0 to t=1. 


I had the same feeling than you have when you write 
   # Now the 'sol' procedure will be
    # re-initialized and re-evaluated
    # between t=0 and t=te, ie t=0 and
    # t=0.1. Since the stopping point for
    # this evaluation is t=0.1, then any
    # evaents which occur after this time
    # will not be accessible in the next
    # iteration of the loop


But I was lured because I thought that sol(1) was in fact a way to tell Maple that the solution has to be evaluated in the range 0..1.
For my defence the message  Error, (in _dat[1]) 'direction' must be set prior to calling/setting 'eventfired'
you get if you omit sol(1) before extracting the fired events is not very explicit (for me at last) .
Moreover if you replace sol(1) by sol(0.2), only the first event is fired, so the value T you use to "set the direction" is not insignificant.
I thought what I called "initialization" was done only once and that command sol(eventfired=[i]) 
meant "extract, in the range 0..T,  the time when the event "i" is fired"

Your "Since the stopping point for this evaluation is t=0.1, then any events which occur after this time will not be accessible in the next iteration of the loop" gave me the idea to add sol(1) as the first command within the loop.
Then all the events are extracted whatever their order of description.


Pushing further I found this is the  the command xe:=sol(te) which is responsible of the non firing of latter events.
If you write (with no loop here)

sol(eventfired=[1]);
sol(op(%));
sol(eventfired=[2]);
sol(op(%));
sol(eventfired=[3]);
sol(op(%));

You will find that events 2 and 3 are not fired which means that sol(some value) is an "active operation" while sol(eventfired=[i]) is a "passive operation".
The good way to proceed thus is:

  1. set 'direction' by writting sol(T) for any T larger than the suspected time of the latest event
  2. extract all the firing times TF[i]
  3. evaluate the solutions x[i]:=sol(TF[i]) for decreasing values of the TFs



All this exchange gave me precious hints to use events in a more carefully way, even if still do not understand exactly what is the logic that Maple uses behind.

Thanks again for you fruitful help.

@Joe Riel 

This is what I would have done myself (maybe with a scan over relative integers?).
Remark that in this case you can use LinearFit after a propre parameterization:

  • let p = x^n
  • --> y=p/(1+m*p)
  • --> (1+m*p)*y=p
  • let q = p*y= y*x^n
  • --> y+m*q=p
  • --> p-y = m*q
  • let r=p-y = x^n-y
  • --> r = m*q

Then this is trivial ordinary least squares regression of r against q whose solution is 
m_optimal = mean(r*q) / mean(q^2)
No interest in Maple beyond simple curiosity. But this can have some interest if you want to solve the problem with Excel.

Remark : I wasn't able to load the data file to confirm this approach (pb with open office on my Mac)

PS : Fitting y=P(x)/Q(x) where P and Q are polynoms can generally be rewritten as fitting u=R(v) where R is a polynom and u and v a reparameterization u=A(x,y), v=B(x,y) with A and B polynoms depending on P and Q.
But this approach sometimes lead to numerical problems (worse conditioning of the information matrix).

@Carl Love 
@tomleslie
@Kitonum

Sorry for missing the answer, I was a two places at once.
Thank to both of you.

To @Kitonum and @Carl Love : No, I hadn't saw that vertical-axis tickmarks had negative values

Incidentally I had the reverse the list of the events... and this gave good results (see the attached file) ????
More generally, from several other examples, it seems that all goes well if the commands sol(eventfired=[i])are executed in an order such that they return a decreasing values of t.

restart;

interface(version);

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

sys := { diff(x(t), t) = 1, x(0) = 0 }:
evs := ListTools:-Reverse([ [x(t)-0.1, none],  [x(t)-0.3, none], [x(t)-0.5, none] ]):
sol := dsolve(sys, numeric, events=evs):

plots:-odeplot(sol, [t, x(t)], t=0..0.5, gridlines=true);

 

# times that fired the events

sol(1): # initialization

sol(eventfired=[1]);
sol(eventfired=[2]);
sol(eventfired=[3]);
 

[HFloat(0.49999999999999994)]

 

[HFloat(0.3)]

 

[HFloat(0.1)]

(2)

# Same times computed  within a loop

for i from 1 to 3 do
  te := op(sol(eventfired=[i]));
end do;

HFloat(0.49999999999999994)

 

HFloat(0.3)

 

HFloat(0.1)

(3)

# Values of x(t) computed  within a loop

for i from 1 to 3 do
  te := op(sol(eventfired=[i]));
# xe := sol(te);             # this doesn't return the correct result,
# xe := subs(sol(te), x(t)); # this doesn't work neither
  xe := eval(x(t), sol(te)); # this doesn't work neither
end do;

HFloat(0.49999999999999994)

 

HFloat(0.5)

 

HFloat(0.3)

 

HFloat(0.3)

 

HFloat(0.1)

 

HFloat(0.1)

(4)

 


 

Download Incomprehensible_Reversed.mw

 

@Christian Wolinski 

Thanks.

Triangulation is fast and not a problem (I believe).

The procedure "hachurer" which generates the hatches is already slow, and even very slow.
Let's say a hatch is the segment of a straight line wich, its orientation being known, is fully characterized by its intercept.
In order to reduce the computational effort I built, for each triangle in the triangulation, the range of intercepts (p_range_t) of the straight lines that cut these triangles. So I'm optimal from this point of view ... but, if a H hathces ar a in given triangle, I solve H systemes of equations with the command
( solve({La(x, y, A, p), lhs(L_sides[nt][n])} union rel[nt], [x, y]) ) 
and I belive this is where the time is lost.
An improvement could be:

   1/ Precompute 3 formal solutions with the command above where lhs(L_sides[nt][n]) is replaced by 2 relations
       corresponding to each pairs of sides of an abstract triangle.
        Lets callSOL12, SOL23, SOL31 these formal solutions and express each of them as a function 
        SOLij(x1, x2, x3, y1, y2, y3, slope, intercept) of the (abstract) coordinates of the 3 vertices of the triangle and the slope and 
       intercept of a straight line (at this point, as you could have noticed, my code doesn't work for vertical line: this come
       from the slope+intercept parameterization of the hatching straight lines, instead of the more general p*x+q*y+r
)

   2/ for each triangle :
      2.1/  split the intercept ranges for this triangle in 3 ranges : the one oh hatches that cut sides 1 and 2, next sides 2 and
              3, and finally sides 3 and 1. Let R12, R23, R31 these ranges
      2.2/  Now evaluate each SOLij(x1, x2, x3, y1, y2, y3, slope, intercepet) for all the straight lines with intercept in the range Rij
              and angle A
   
With this there is only 3 calls to solve instead of hatch-number.
I think this is the most importent algorithmic improvement.
Beyond that, my programming skills are quite limited and a Maple's expert would surely write something more efficient.

 

For the texting part i think that a prior selection of the elementary patterns that cut a triangle would be a good idea.
Here the slowness comes from the many calls to the procedure plots:-inequal.
An algortihmic improvement would consist to avoid calling this procedure for patterns we know are entirely in a triangle or entirely outside. Easy to say, but more complicated to write.
A possibility could be to use simplex:-convexhull: let T the list of points of a triangle,P the list of points of a pattern (=motif) and W the list of the whole points:
     if simplex:-convexhull(W) =simplex:-convexhull(T) the pattern is entirely within the triangle
     if simplex:-convexhull(W) =simplex:-convexhull(P) the triangle is entirely within the pattern


 

@Christian Wolinski 

Hi, I've just posting my work about hatching and texturing.
You can access it here Triangulation, hatching and texturing of simple...

It's more a bunch of ideas that a ready-to-use application: a lot of work remain to be done to make it computationally efficient. But I'd like to think it would serve as a starting point.

@Christian Wolinski 

I will post today the work I've done on this subject (it's about seven am here and I'm going to work).
I did it gor fun and so it's very imperfect and quite slow. I never had the time to make it better nor faster but I think that some people here could fix that if this stuff presents any interest.
When back home I will add a few comments and deliver the module.

It enables 3 things :
  1/ tesselation of a simple polygon by triangles (Delaunay triangulation triangulates the convex hull og such a polygon)
  2/ hatching
  3/ simple texturing (the pattern nust necessarily be a polygon, not necessarily convex [I was working recently with pentaminos but my original goal was to do Penrose's meshing of arbitrary symple polygons])

See you soon.

@Joe Riel 

Thanks Joe.

@acer 

Your approach is, from far, more intelligent than mine. 

"I'm pretty sure that I've duplicated a similar..." yes, I saw a lot replies you made (even to me) to customize plots by adding special characters (in my case [or maybe it was sand15?] it was to rename vertices of a graph: I saved your mw file but I can't put my finger on your answer )

"However none of this helps..." right, after tour first answer I just went through the help pages and discovered some stuff around MathML: then I began to play.

@Joe Riel  @rlopez


I agree one hundred percent.
Nevertheless, codegen:-MathML or MathML[ExportPresentation] can help to understand how an expression is represented in MathML syntax.
Here is an unpretentious piece of code that allows you to translate a (simple) Maple expression into this "abstruse lingo" rlopez has evoked.


 

restart;


The procedure translate does not account for all the situations, but it
can easily be completed.

translate := proc(f)
local s, t;
uses StringTools:
s := MathML[ExportPresentation](f):
t := SubstituteAll(s, "</mi>", """),"):
t := SubstituteAll(t, "<mi>", "mi("""):
t := SubstituteAll(t, "</mn>", """),"):
t := SubstituteAll(t, "<mn>", "mn("""):
t := SubstituteAll(t, "</mo>", """),"):
t := SubstituteAll(t, "<mo>", "mo("""):
t := SubstituteAll(t, "<mrow>", "mrow("):
t := SubstituteAll(t, "</mrow>", ")"):
t := SubstituteAll(t, "<msub>", "msub("):
t := SubstituteAll(t, "</msub>", ")"):
t := SubstituteAll(t, "<msup>", "msup("):
t := SubstituteAll(t, "</msup>", ")"):
t := SubstituteAll(t, "<mfrac>", "mfrac("):
t := SubstituteAll(t, "</mfrac>", ")"):
t := SubstituteAll(t, "<mfenced>", "mfenced("):
t := SubstituteAll(t, "</mfenced>", ")"):
t := SubstituteAll(t, "),)", "))"):
t := SubstituteAll(t, ")m", "),m");
return t
end proc:


One example

f := Int(1/2*x^2, x):
translate(f)

"<math xmlns='http://www.w3.org/1998/Math/MathML'>mrow(mo("&Integral;"),mrow(mfrac(mn("1"),mn("2")),mo("&InvisibleTimes;"),msup(mi("x"),mn("2"))),mo("&InvisibleTimes;"),mrow(mo("&DifferentialD;"),mi("x")))</math>"

(1)


Now copy the substring between the first > and defore the second < (excluded and
paste it while enclosing it between `# and `
You must obtain this line, next click enter.

`#mrow(mo(\"&Integral;\"),mrow(mfrac(mn(\"1\"),mn(\"2\")),mo(\"&InvisibleTimes;\"),msup(mi(\"x\"),mn(\"2\"))),mo(\"&InvisibleTimes;\"),mrow(mo(\"&DifferentialD;\"),mi(\"x\")))`

`#mrow(mo("&Integral;"),mrow(mfrac(mn("1"),mn("2")),mo("&InvisibleTimes;"),msup(mi("x"),mn("2"))),mo("&InvisibleTimes;"),mrow(mo("&DifferentialD;"),mi("x")))`

(2)


and another one

f := (x__2+y^3)/(1+cos(y)):
translate(f)

"<math xmlns='http://www.w3.org/1998/Math/MathML'>mfrac(mrow(msup(mi("y"),mn("3")),mo("+"),mi("x__2")),mrow(mn("1"),mo("+"),mrow(mi("cos"),mo("&ApplyFunction;"),mfenced(mi("y")))))</math>"

(3)

`#mfrac(mrow(msup(mi(\"y\"),mn(\"3\")),mo(\"+\"),mi(\"x__2\")),mrow(mn(\"1\"),mo(\"+\"),mrow(mi(\"cos\"),mo(\"&ApplyFunction;\"),mfenced(mi(\"y\")))))`;

`#mfrac(mrow(msup(mi("y"),mn("3")),mo("+"),msub(mi("x"),mi("2"))),mrow(mn("1"),mo("+"),mrow(mi("cos"),mo("&ApplyFunction;"),mfenced(mi("y")))))`

(4)

 


 

Download trick.mw

@acer 

I've forgotten the lprint command...

About you programmatic construction: I never suspected it was possible!

Thank you very much

@acer @rlopez

(sent from my personal account)

Such a difference of point of view!
At first sight I would say that acer's answer suits me perfectly: on one side I have rejected the 2D mode long ago, and on the other side acer already familiarised me with MathML through his replies to a lot of questions I already asked concerning labels (graph vertices or plots).
But your trick concerning this conversion to atomic variables that I wasn't aware at all is impressing.
So thumb up for both of you.


By the way, can anyone of you explain me this strange result?

 

restart

e := { `#mrow(msub(mi("H"),mn("2")),msub(mn("0"),mn("2")))`, `#mrow(msub(mi("Fe"),mn("2")),msub(mn("0"),mn("3")))`
}

{`#mrow(msub(mi("Fe"),mn("2")),msub(mn("0"),mn("3")))`, ` #mrow(msub(mi("H"),mn("2")),msub(mn("0"),mn("2")))`}

(1)

e := { NULL, `#mrow(msub(mi("H"),mn("2")),msub(mn("0"),mn("2")))`, `#mrow(msub(mi("Fe"),mn("2")),msub(mn("0"),mn("3")))`
}:
f := e[2..-1]

{`#mrow(msub(mi("Fe"),mn("2")),msub(mn("0"),mn("3")))`, `#mrow(msub(mi("H"),mn("2")),msub(mn("0"),mn("2")))`}

(2)

 


 

Download AtomicVariables.mw

@acer 

Thanks for giving links to all these posts.
I've just given a quick look and all this stuff seems very interesting: I'm going to read them carefully.

@Carl Love @Kitonum

Thank you very much to both of you.

@epostma 

Sorry for this late reply.

Thanks for all these precious informations.
Should I conclude that almost all random variables, with the exception of those in the Gaussian family, have a specific sampler?

First 116 117 118 119 120 121 122 Last Page 118 of 154