mmcdara

7891 Reputation

22 Badges

9 years, 60 days

MaplePrimes Activity


These are replies submitted by mmcdara

I've written this piece of code xhich seems to work correctly.
Nevertheless I'm not one hundred percent sure it still works correctly for other expressions of the same kind.
Probably this could be written more efficiently...


 

restart


The black code above comes from this Mapleprimes' Post
(personal modifications are pink highlighted)

with(DocumentTools):

with(DocumentTools:-Layout):

F := Font("Liebniz notation", size=12, color=blue, style=:-Hyperlink):

H := DocumentTools:-Layout:-Hyperlink(F,linktarget="https://www.mapleprimes.com/maplesoftblog/210289-Writing-Derivatives-At-A-Point-Using"):

InsertContent(Worksheet(Group(Input(Textfield("Original post from Technical Support ",H))))):

 

# Type to check if an expression is a derivative using 'D', e.g. D(f)(a) and D[1,2](f)(a,b).

 

TypeTools:-AddType(   

 

        'Dexpr',      

 

        proc( f )     

 

               if op( [0,0], f ) <> D and op( [0,0,0], f ) <> D then

 

                       return false;

 

               end if;       

 

               if not type( op( [0,1], f ), 'name' ) or not type( { op( f ) }, 'set(algebraic)' ) then

 

                       return false;

 

               end if;       

 

               if op( [0,0,0], f ) = D and not type( { op( [0,0,..], f ) }, 'set(posint)' ) then

 

                       return false;

 

               end if;       

 

               return true;          

 

        end proc      

 

):

 

 

# Create a local version of 'print', which will print expressions like D[1,2](f)(a,b) in a custom way,

 

# but otherwise print in the usual fashion.

 

Dprint := proc()

 

 

        local A, B, f, g, L, X, Y, Z;

 

 

        # Check that a valid expression involving 'D' is passed, along with a variable name or list of variable names.

 

        if ( _npassed < 2 ) or ( not _passed[1] :: 'Dexpr' ) or ( not passed[2] :: 'Or'('name','list'('name')) ) then

 

               #return :-print( _passed );

               return  _passed;

 

        end if;

 

 

        # Extract important variables from the input.

 

        g := _passed[1]; # expression

 

        X := _passed[2]; # variable name(s)

 

        f := op( [0,1], g ); # function name in expression

 

        A := op( g ); # point(s) of evaluation

 

 

        # Check that the number of variables is the same as the number of evaluation points.

 

        if nops( X ) <> nops( [A] ) then

 

               #return :-print( _passed );

               return  _passed;

        

        end if;

 

 

        # The differential operator.

 

        L := op( [0,0], g );

 

 

        # Find the variable (univariate) or indices (multivariate) for the derivative(s).

 

        B := `if`( L = D, X, [ op( L ) ] );

 

 

        # Variable name(s) as expression sequence.

 

        Y := op( X );

 

 

        # Check that the point(s) of evaluation is/are distinct from the variable name(s).

 

        if numelems( {Y} intersect {A} ) > 0 then

 

               #return :-print( _passed );

               return  _passed;

        

        end if;

 

 

        # Find the expression sequence of the variable names.

 

        Z := `if`( L = D, X, X[B] );

 

        #return print( Eval( Diff( f(Y), Z ), (Y) = (A) ) );

        return Eval( Diff( f(Y), Z ), (Y) = (A) );

 

 

 

end proc:

# Example

expr := -(Diff(f(mu__1, mu__2), mu__1))^2*mu__1^2-(Diff(f(mu__1, mu__2), mu__2))^2*mu__2^2+2*(Diff(f(mu__1, mu__2), mu__1))^2*lambda__1^2+2*(Diff(f(mu__1, mu__2), mu__2))^2*lambda__2^2+2*(Diff(f(mu__1, mu__2), mu__1))*(Diff(f(mu__1, mu__2), mu__2))*lambda__1*lambda__2-2*(Diff(f(mu__1, mu__2), mu__1))*mu__1*(Diff(f(mu__1, mu__2), mu__2))*mu__2 + anything_else;

-(Diff(f(mu__1, mu__2), mu__1))^2*mu__1^2-(Diff(f(mu__1, mu__2), mu__2))^2*mu__2^2+2*(Diff(f(mu__1, mu__2), mu__1))^2*lambda__1^2+2*(Diff(f(mu__1, mu__2), mu__2))^2*lambda__2^2+2*(Diff(f(mu__1, mu__2), mu__1))*(Diff(f(mu__1, mu__2), mu__2))*lambda__1*lambda__2-2*(Diff(f(mu__1, mu__2), mu__1))*mu__1*(Diff(f(mu__1, mu__2), mu__2))*mu__2+anything_else

(1)

LiebnizNotation := proc(expr, vars)
  local DiffTerms, Remainder, opDiffTerms, S, A, u, v:
  DiffTerms   := map(u -> if has(u, 'Diff') then u end if, [op(expr)]):
  Remainder   := expr-add(DiffTerms);
  opDiffTerms := map(u -> [op(u)], DiffTerms):
  S := 0:
  for u in opDiffTerms do
    A := 1:
    for v in u do
      if has(v, 'Diff') then
        patmatch(v, (d::algebraic)^n::integer,'la');
        A := A * Dprint(convert(eval(d, la), 'D'), [x__1, x__2])^eval(n, la);
      else
        A := A*v
      end if:
    end do:
    S := S+A:
  end do:
  return S + Remainder
end proc:

LiebnizNotation(expr, [u__1, u__2])

-(Eval(Diff(f(x__1, x__2), [x__1]), (x__1, x__2) = (mu__1, mu__2)))^2*mu__1^2-(Eval(Diff(f(x__1, x__2), [x__2]), (x__1, x__2) = (mu__1, mu__2)))^2*mu__2^2+2*(Eval(Diff(f(x__1, x__2), [x__1]), (x__1, x__2) = (mu__1, mu__2)))^2*lambda__1^2+2*(Eval(Diff(f(x__1, x__2), [x__2]), (x__1, x__2) = (mu__1, mu__2)))^2*lambda__2^2+2*(Eval(Diff(f(x__1, x__2), [x__1]), (x__1, x__2) = (mu__1, mu__2)))*(Eval(Diff(f(x__1, x__2), [x__2]), (x__1, x__2) = (mu__1, mu__2)))*lambda__1*lambda__2-2*(Eval(Diff(f(x__1, x__2), [x__1]), (x__1, x__2) = (mu__1, mu__2)))*mu__1*(Eval(Diff(f(x__1, x__2), [x__2]), (x__1, x__2) = (mu__1, mu__2)))*mu__2+anything_else

(2)

 

 


 

Download Liebniz_Notation.mw

@mmcdara 

To complete Carl's reply look here  https://www.random.org/randomness/

(and possibly here for more details on PRNG https://en.wikipedia.org/wiki/Pseudorandom_number_generator)

@digerdiga 


randomize() is just a command that "randomly" generates a seed and has nothing to do with rand(...)
randomize() is useful  if you do not want to generate ALWAYS the same random numbers  (note that Seed:=randomize() is useless,  unless you want to reuse the same seed "Seed" to generate the same random numbers again)

To see this execute many times this block of instruction  

restart;
Seed := randomize();
rfloat := rand(0. .. 1.0);
rfloat();


and compare the results to many executions of this one

restart;
rfloat := rand(0. .. 1.0);
rfloat();



 

@Mariusz Iwaniuk 

Good, thank you Mariusz

@acer 

Great... but I would understand that this holiday season, you have better things to do than spend time answering me. 

Best regards

@acer 

Thanks for your quick reply.
It's a very interesting feature for which I will find to find a lot of uses.

If it's not too much to ask, how can I define now  a link to the help pages?

My purpose is to display a table which presents a few elementary characteristics of the different continuous random variables implemented in Maple. Its first column is aimed to contain the name of the distribution and I'd like this name to appear as an hyperlink to the corresponding help page.

Here is a file which will give you a better idea of what I try to do.

ContinuousRandomVariables.mw

Thanks in advance

Merry Christmas

@Carl Love 

When I was a kid I was fascinated by a living room lamp in a shop window, made of a long transparent tank that contained a green liquid and contiously oscillated around a horizontal axis.
Years later, as I was designing front tracking methods, I was led to simulate this kind of phenomenon (2D), called sloshing.
One way to simulate sloshing is to solve Navier-Stokes equations, but is extremely computational demanding, specially to simulate 3D sloshing.
Probably the most commonly used approach today (at least to create virtual reality films, which is quite different than the computation, for instance, of resonance frequencies, or damping times) is SPH (Smooth Particle Hydrodynamics)

https://www.youtube.com/watch?v=Lo3H7Ut6iyM
https://www.youtube.com/watch?v=iHACAlfYeiQ

I thing it would be a really challenging, but motivating,  problem to simulate sloshing, even in 2D, with MAPLE.

 

@Rouben Rostamian  

In case the file os not correctly loaded here is the corresponding code 

restart;
with(plots):

de := diff(x(t),t) = v(t),
      diff(v(t),t) = 1/m*F(x(t), v(t), u);

FrictionCoefficient := c -> (mu__s-(mu__s-mu__k)*(1-exp(-abs(xi*c)))) * tanh(c/epsilon) ;
plot(eval(FrictionCoefficient(c), [mu__s=0.8, mu__k=0.5, epsilon=0.001, xi=50]), c=-0.1..0.1);


F := (x, v, u) -> -k*x - m*g*FrictionCoefficient(v-u):

ic := x(0)=2, v(0)=0;

params := m=1, k=1, g=1, mu__k=0.5, mu__s=0.8, u=0.25, epsilon=0.001, xi=20;

tmax := 40;
%dsolve({de,ic}, numeric, output=operator, range=0..tmax):
eval(%, {params}):
dsol := value(%);

myx := eval(x, dsol):
myv := eval(v, dsol):

%odeplot(dsol, [[t, x(t)],
                      [t, v(t)],
                      [t, F(x(t),v(t),u)],
                      [t, mu__s*m*g/k],
                      [t,-mu__s*m*g/k]
                     ], color=[red,"Green", blue, gray, gray],
                     refine=1, size=[800,300],
                     legend=[displacement,velocity,force,"",""]):
subs(params, %):
value(%);

@Rouben Rostamian  

Hi Rouben, 

Here is the implementation of the smoothed friction model.
It depends on two extra parameters :

  • epsilon smoothes the transition between negative and positive values of v-u
  • xi smoothes the transition between pu__s and mu__k (you could find info on wikipedia about smoothed Coulomb's model frictions)

With this model events are no longer necessary which (for my experience events are sometimes the cause of earlier stopping of the simulation).

You will see that m=1 works correctly.
Of course, in a real application, epsilon and xi would have to be adjusted correctly... here it is a simple illustrative choice.

Now an important point: smooth solid friction models are not just mathematical tricks to remove discontiniuties, they are models that mimic the reality more than models with discontinuitues that, in fact, are over simplified approximations of the reality.
It's normal that, as you write  "The discrepancy [between the unsmoothed and the smoothed friction models] is most noticeable in the plot of v(t).  The original v(t) has flat parts, indicating a motion with constant velocity (sticking to the belt)..." but the question we hsould ask to ourselves is "what is the more realistic behaviour (let's say the one we could observe during an experiment): the one with an idealized friction model or the smoth one (after a suitable tunning of its coefficients)".
If we adopt the purely mathematical point of view of solving an ode with discontinuous coefficient, the question of using smooth friction models is stupid.
But if want to solve an ode which is more representative of some observed phenomenon, then maybe it worths using a smooth model (or more complicated ones, see for instance http://www.dymoresolutions.com/resume/publications/DoFerriBauchau07.pdf
).

--------------------------------------------------------------------
Let's share experiences
To give you a quick idea of the problems  I use to solve, think to a closed box moving under accelerations conditions (the problem is more easily picture in 1D) with a mass inside which can have displacements relative to the box.
A prestressed spring maintains the mass against a wall of the box until the inertia force is large enough its effects. Then the mass moves submitted to the effect of the spring, friction effects and hydraulic damping.
The mass can reach the opposite wall, return to its original rest position, and when this happens, it can stay here until the inertia force enables it to move again.
Finally not so different to your own problem.
To solve this here are a few things I use :

  • option parameters in dsolve
  • method = rosenbrock
  • two events  [ [x(t)-RightWall, v(t) > 0], halt],  [ [x(t)-LeftWall, v(t) < 0], halt] (among others)
  • as soon as one of these two events is fired I enter a specific procedure where:
    1. Let FLY the name of this solution.
    2. Let Tstop the time the event fired.
    3. I compute the time when the mass will move again (this time is the solution of a static equilibriium relation obtained by equating v(t) to 0 and x(t) to the abscissa of the wall the contact to fired the event).
      Let Tstart this time
    4. In the range [Tstop, Tstart] the solution is [x(t) = RightWall or LeftWall, v(t)=0]
      Let REST the name of this solution
    5. The solution is piecewise(t < Tstop, FLY(t) , if t < Tstart, REST(t))
    6. Go to point 1 while the maximum time has not been reached
  • At the end the whole solution is a sequence i-of alternate FLY and REST phases

Building it doesn't take more time than building it in one single shot. The advantage is that you can do a fine tunning of the events.
Maybe this will give you some ideas to improve your own approach?
--------------------------------------------------------------------


Here I used a smooth model as a simple trick to help Maple solving the ode an doing this, I agree, I do not solve the same problem than you did.


Unfortunately I can't load the content of the file, just its link : given the message I got, please tell me is the file is readable or if it is corrupted.

Improvement.mw
Failed to load the worksheet /maplenet/convert/Improvement.mw .

Nice.

In your mw file you wrote "Note:  This demo is not without its flaws.  Depending on the choice of the parameter values....  I don't know how to get around that problem"

I use to simulate mechanical systems with friction models of Coulomb's type or of more complex types (models that generate stick-slip phenomena, see for instance https://en.wikipedia.org/wiki/Stick-slip_phenomenon if you want to know more about them).

A key trick in the numerical treatment of all these friction models is to smooth the discontinuity at v(t)=0, for instance by writing 
friction = mu * tanh( v(t) / eps )
instead of 
friction = piecewise(v(t) > 0, mu*signum(v(t)), 0)

The higher eps (homogeneous to a velocity) the smoother the transition between v(t) < 0 and v(t) > 0.
I don't know if this trick will help you but, I repeat, it's practically used in any computatoinal code involving non-reversible solid friction models.


 

 

I repeated what I did in my original post, now to the case the dice is modelled with a ProbabilityTable random variable.
You will see that we still have incorrect results for ProbabilityFunction and Support (even if, rather strangely, their expressions are different for those for EmpiricalDistribution and DiscreteUniform ... which could denote a coding from different persons at different times).

The funny thing is that probably someone on the development team must have noticed  that something went wrong because the plots of the probability (mass) function are correct in the MathApps in help(Statistics[StatisticsAndProbability) (select next Discrete random Variables).
Maybe he thought it was a simple problem of graphical representation and he didn't report the information that something was wrong?

restart:


A quick reference to the terminology (in particular for discrete RVx)

https://en.wikipedia.org/wiki/Probability_distribution#Functions_for_Discrete_Variables

with(Statistics):

P := [(1/6)$6];
AnotherDice := RandomVariable(ProbabilityTable(P)):

[1/6, 1/6, 1/6, 1/6, 1/6, 1/6]

(1)

exports([attributes(AnotherDice)][3]);

Conditions, ParentName, Parameters, CDF, CentralMoment, CharacteristicFunction, CGF, Mean, Median, Mode, MGF, ProbabilityFunction, Support, RandomSample, RandomVariate

(2)

# Given the result of the incorrect plot I guess ProbabilityFunction refers to the Probability Mass Function
# defined in the reference above..
#
# Its definition is incorrect
# Its plot is incorrect too

ProbabilityFunction(AnotherDice, x);
plot(ProbabilityFunction(AnotherDice, x), x=0..7);

piecewise(x < 1, 0, x <= 6, [1/6, 1/6, 1/6, 1/6, 1/6, 1/6][floor(x)], 0)

 

 

# Wrong results.
# The correct result is [ 0, 0, 1/6, 0, 1/6, 0, 1/6, 0, 1/6, 0, 1/6, 0, 1/6, 0 , 0 ]

f := unapply(ProbabilityFunction(AnotherDice, x), x):
map(f, [seq(0..7, 1/2)])

[0, 0, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, 0, 0]

(3)

# Wrong result.
# The support is the set {1, 2, 3, 4, 5, 6};

Support(AnotherDice)

RealRange(1, 6)

(4)

 


 

Download ProbabilityTable_RV.mw

@Carl Love 

Thanks Carl, you did well (for displacing it to a post of course).
 

@acer 

I spend a lot of time working with the Statistics package and I had never seen this ExploreRV command before. (Why is it not referenced? It would at least deserve to be put in the applications of the start page) 

Damn me, I didn't notice ExploreRV is in the Student package!
It seems very fine and it could suit me well.
Thanks for the tip.

I often encounter a lot of problems with Explore, a powerful by rather delicate tool to use.I'm sure your advices will be extremely useful to me.

Thanks for all, see you next time.

@Carl Love 

Ok, thanks.
The good point is that works with Maple 2019. 
I will test this tomorrow at work.

If it's not too much to ask, the attached file contains an example of what I really wanted to do.
This example works perfectly but I need to change several little things as soon as I change the random variable VA.
Explore_Initial_Attempt.mw

So I had the idea to embed this piece of code within a procedure OBSERVE(RV, SliderRanges) where RV is a random variable and SliderRanges a list of ranges that must contain as many elements as RV has parameters.
In my mind this procedure would be called by writting for instance
OBSERVE(RandomVariable(Normal(a, b)), [-2.0..2.0, 1.0..3.0] ) where -2.0..2.0 denotes the range of "a" and 1.0..3.0 the range of '"b".

Here is a sketch of the procedure I wrote and from which I extracted the piece of code of my initial question:
Explore_Procedure.mw

Do you think this could work ?

 

@Carl Love 

Sorry Carl, I've just seen I forgot to define f.
Your reply and my correction crossed each other

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