mmcdara

6650 Reputation

18 Badges

8 years, 132 days

MaplePrimes Activity


These are replies submitted by mmcdara

@acer 

IEEE 75says that elementary operations, sqrt included, must be correctly rounded.
In my case I have a compound expression involving two sqrt and a multiplication.
Each of them has to be correctly rounded.

Ina "classical" programming language (Fortran, C, ...) all this stuff is implicitely done by the compiler and the precessor behind.
In such a lanquage something like a = sqrt(2) * sqrt(3) is decomposed i, a syntactic graph and each meave or branch is correctly rounded in the IEEE sense.
At the very end the rounding error is propagated along the graph up to the root "a".
Then a could be written as  round( round(sqrt(round(2)) *  round(sqrt(round(3)) ).

What puzzled me is that a := evalf[8](sqrt(2) * sqrt(3)) seems to hide this kind of recursive rounding and looks like some composite expression wher evalf[8] is the last function to apply.

When I write  a := sqrt(2.0)*sqrt(3.0), a is not a 8-digits evaluation of sqrt(6) but a 10-digits one. Then taking evalf[8](a) operates on a sinle number and rounds it correctly (=2.4494897).
But if I consider that MAPLE acts as Fortran or C (for instance), then it applies the rounding on each leave and branch of the syntactic graph of a = sqrt(2) * sqrt(3). Then a =evalf[8]( evalf[8](sqrt(evalf[8](2.0))) * evalf[8](sqrt(evalf[8](2.0)))  ) = 2.4494898

evalf[8](sqrt(evalf[8](2.0))):
evalf[8](sqrt(evalf[8](3.0))):
evalf[8](%*%%) ;
              2.4494898

 

So, ok you get the point and I was wrong.

Just a last question (put to the trash the question I asked in may last mail [evaluation.mw]) :  to understand how MAPLE actually works on this problem, I did this 

f := proc(x, y)
evalf[8](sqrt(x)*sqrt(y))
end proc;
proc(x, y)  ...  end;
dismantle(op(f));

PROC(11)
   EXPSEQ(3)
      NAME(4): x
      NAME(4): y
   EXPSEQ(1)
   EXPSEQ(1)
   EXPSEQ(1)
   FUNCTION(3)
      TABLEREF(3)
         NAME(4): evalf #[protected]
         EXPSEQ(2)
            INTPOS(2): 8
      EXPSEQ(2)
         PROD(5)
            FUNCTION(3)
               NAME(4): sqrt #[protected, _syslib]
               EXPSEQ(2)
                  PARAM(2): [1]
            INTPOS(2): 1
            FUNCTION(3)
               NAME(4): sqrt #[protected, _syslib]
               EXPSEQ(2)
                  PARAM(2): [2]
            INTPOS(2): 1
   EXPSEQ(1)
   EXPSEQ(1)
   EXPSEQ(1)
   BINARY(2)
      2
   EXPSEQ(3)
      LIST(2)
         EXPSEQ(5)
            STRING(4): ""
            INTNEG(2): -1
            INTPOS(2): 4
            INTPOS(2): 6
      LIST(2)
         EXPSEQ(5)
            STRING(4): ""
            INTNEG(2): -1
            INTNEG(2): -1
            INTPOS(2): 23

I can understand some part of the output "dismantle" returns but not the full output : where can I find informations to analyse the output.

Thanks for all

 

 

@acer 

I agree sqrt(2.0)*sqrt(3.0) is a compound operation.
I'm going to verify your answer

In the help pages for printf it's written that 
The fprintf command is based on a C standard library command of the same name

Would you be kind enough to explain me the results in the attached file ?
It seems to me that columns 1 and 3 are "in my words correct" and that column 2 is not.
Moreover, I understand the differences between columns 2 and 3 as different "floating point evaluations" for evalf and the C standard library cited above

evaluation.mw
 

restart:

interface(version);

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

(1)

Digits := 10:
for n from 10 to 3 by -1 do
  #blk := cat(seq(" ", k=1..11-n));
  fmt := cat("%-12a  |  %-12a  |  ", "%-3.", n-1, "f\n"):
  printf(fmt,
         evalf[n](sqrt(6.0)),
         evalf[n](sqrt(2.0)*sqrt(3.0)),
         sqrt(2.0)*sqrt(3.0)
  )
end do:

2.449489743   |  2.449489743   |  2.449489743
2.44948974    |  2.44948974    |  2.44948974
2.4494897     |  2.4494898     |  2.4494897
2.449490      |  2.449491      |  2.449490
2.44949       |  2.44948       |  2.44949
2.4495        |  2.4495        |  2.4495
2.449         |  2.449         |  2.449
2.45          |  2.44          |  2.45

 

 

 


Thanks in advance

Download evaluation.mw

 

 

 

What you say is :

  1. X = evalf[3](2) = 1.41
  2. Y = evalf[3](2) = 1.73
  3. X*Y / 10000 = 2.4393
  4. Then : Why evalf[3](6) would give 2.45 ?


This is not a demonstration to the fact that evalf[3](sqrt(2)*sqrt(3)) MUST return 2.44
What you say here is that 2.44 = evalf[3](sqrt(2)) * evalf[3](sqrt(3)), which is obviously not in agreement with the  IEEE 754 standard for arithmetic of computers.

What this standard says is that evalf[3](sqrt(2)*sqrt(3)) should return the value of evalf[3](sqrt(6)), then 2.45
If it is not the case, this is because MAPLE interprets evalf[3](sqrt(2)*sqrt(3))  as evalf[3](sqrt(2)) * evalf[3](sqrt(3)).

The IEEE754 standard says exactly what elementary mathematics says about composition of function : the inner operation (here `*`) must be evaluated before the outer one (evalf[3]).
 

I now the world « bug » is taboo here and I haven’t pronounce it. But if you want to prove me MAPLES does the thing correctly and I'm wrong (a thing  I am willing to accept because nobody, me in particular, is infallible), please give me a more solid explanation.

@Preben Alsholm 


Thanks for your answer Preben.
It lets me unsatisfied but  I should have probably expressed my question differently.

 

Please take a look at that:

for n from 10 to 3 by -1 do

  [evalf[n](sqrt(6.0)),evalf[n](sqrt(2.0)*sqrt(3.0))]

end do;

                   [2.449489743, 2.449489743]

                    [2.44948974, 2.44948974]

                     [2.4494897, 2.4494898]

                      [2.449490, 2.449491]

                       [2.44949, 2.44948]

                        [2.4495, 2.4495]

                         [2.449, 2.449]

                          [2.45, 2.44]

 

The operation evalf[n](sqrt(6.0)) respects the IEEE 754 « rounding to the nearest » rule (aka « RN » rule)

But the result of evalf[n](sqrt(2.0)*sqrt(3.0)) doesn’t seem to respect any rule (not the RN rule, neither the DR « towards -infinity » rule, 

neither the RU « towards +infinity » rule)

 

The only explanation I can give for explaining the differences in the lines above has to be searched in the way Maple does the operationevalf[n](sqrt(2.0)*sqrt(3.0)) 

 

Let E[n] the operation evalf[n], X = sqrt(2.0) and Y = sqrt(3.0)

In evalf[n](sqrt(2.0)*sqrt(3.0))  MAPLE does not realize the compound operation E[n](`*`(X, Y)) but the operation

`*`(E[n](X), E[n](Y)).

Is it correct ?

 

 

Let O some arithmetic operation, O[infinity] its result in infinite precision and O[n] its result in finite « n digits » precision 

The IEEE 754 standard requires that the result of O[n](X,Y) be equal to E[n](O[infinity](X,Y)).

Here 

 

  • IEEE 754              :  O[infinity](`*`(X,Y) = sqrt(6.0)  
  • IEEE 754 RN rule :  O[8](`*`(X,Y)) = 2.4494897  
  • MAPLE                 :  E[8](`*`(X,Y)) i= 2.4494898

 

Unless I'm mistaken MAPLE does not respect the IEEE 754 standard  ???

It seems so odd that I can believe that.

 

Any answer ?

@acer  @vv

Thanks for the this little explanation.
By the way you two made me discover the function "value" and its key role in some circumstances.


From the value help pages I read "During computations, inert functions remain unevaluated (the operations they represent remain unperformed). The value command maps these inert functions, such as Int or %exp, into the corresponding active functions int and exp"
If I'm not mistaken value(r) here has the same role convert(r, sum) had in my code ?

Happy new year to both
 

@vv 

It was not at all obvious !
Thank you for the answer

@Kitonum 
Lack of concentration at the beginning of this year.

Surely this would have been a better answer 
   NewList := proc(L::list, treshold::{integer, float})::list:
      map(u -> if treshold <u then 2*u else u end if, L):
   end proc:


Thanks for you exercised eye.
By the way, I use to use map instead of `if ` (contrary to what is usually done here) : are the twos equivalent or does it exist some subtleties ?

Happy new year
 

Point 1 : I guess (see Kitonum's reply) a[1], a[3], p[3] are arbitrary functions of t ?
In this case you should write your ode this way :
edo := diff(r(t), t)-(diff(a__3(t), t))*r(t)/a__3(t)+a__1(t)*p__3(t)-3*r(t)*r(t) = 0;

Point 2 : every time dsolve doesn't give an answer, use infolevel :
Here 
infolevel[dsolve] := 4:
dsolve(edo, r(t))

What happens ? Maple works hard but doesn't find any solution.
Let us rewrite your ode on the form  edo := diff(r(t), t)+f(t)*r(t)+g(t)+c*r(t)^2 = 0;
This a non linear first order edo of Ricatti type (use infolevel[dsolve] := 4 to verify) and the problem is that, apart from special cases, there is no general technique for obtaining a solution (you can refer to "Advanced Mathematical Methods for Scientists and Engineers", Bender & Orsazg, paragraph 1.6 [the one the red text here comes from] ).
In this reference you will see that "many Ricatti equations can be solved" ... as soon as you are able to "guess just one [particular] solution"
In your case the 3 functions a[1], a[3], p[3] not being in closed form, there exist no mathematical way to obtain a solution !
 

@Les 

just replace 

   DrawGraph(subs(DrawRule(M), M), style=tree, root="World") ;

by 
   DrawGraph(subs(DrawRule(M), M), style=tree, root="World") :
   PXY := GetVertexPositions(M, style=tree, root="World"):
   PYX := map(u -> [-u[2], u[1]], PXY):
   SetVertexPositions(M, PYX):
   print(DrawGraph(subs(DrawRule(M), M), scaling=unconstrained))

A last problem remains : the yellow boxes that contain the labels are not large enough for these labels to appear.
I do not see how to fix this in a simple way ... maybe you will find some way on your side ?
 

 

@Les 

I took your reply as an incentive to do better.
Here is a code that generates the permutation graph for any number of elements to permute.
If this number is large, it would be better to swap the X and Y coordinates before drawing the graph (the graph would extend from left to right, not from top to down)

The code is rather simple and still based on Robert Israel's trick.
I use here a "second trick" to avoid the complex management of the numbered vertices the original trick seems to impose : each vertex (a single letter) is extended by a random integer number.
For instance "a" becomes cat("a", rand()).
It works perfectly well ... while the same extended name is not generated twice !

One solution (as suggested before) is to construct recursively the permutation graph: but I'm not familiar with recursivity and it is always a pain for me to write recursive procedures.
So the code begins by constructing the trivial "order 1" permutation graph "World"--"a", and progressively augment it until all the elements ("b", "c", ...) are included.
I don't know if this  strategy icompares to a recursive one in terms of performances but, I believe it is simpler to read and maintain.
 

PermutationGraphs.mw

@Joe Riel 

Here is a copy paste from the Maple worksheet
X . Vector[column](3, [1, 1, -1]) ; 
X . Vector[column](3, [-1, 2, 0]) ;

By the way, does it exist a way to vizualize extra characters in a worksheet ?

 

 

@tomleslie 

I'm sorry :

two options :

  1. You have read the description of the physical problem I gave to JhonS
    1. I misspoke and you misunderstood the problem: it's normal that you provide me a solution which does not correspond to the one I'm expected for.
      If it is so, my apologies for my poor explanations
    2. My explanations are correct but did not understand the problem and you keep going on your first idea  that I shouldn't have put the problem the way I did.
      In this case, with all due respect, I don't see any good reason to keep talking about it
  2. You did not
    In this case I don't care about the solution you give in eventsProb.mw : it is not the solution of the problem I'm concerned with.

Thank you, however, for the yime you have devoted to me

 

@tomleslie 
Please look to my answer to JohnS, maybe it will enlighten you.


I just have two events.
The first one one [[x(t)=0, v(t)<0], [x(t)=0, v(t)=0]] plays the same role than the event described in the help pages for the bouncing ball test case. Excepted I do not reverse the velocity here, but just force it to be 0 (think not to a boucing ball but o an egg crashing on the floor).

The events are correctly specified: see help pages, example 1 :
dsn := dsolve({diff(y(t),t,t)+y(t)=0,y(0)=0,D(y)(0)=1}, numeric, events=[[[y(t),diff(y(t),t)>0],halt], [[y(t),diff(y(t),t)<0],halt]]);

You write "If I replace your 'events' option with the simpler events=[[x(t)=0, v(t)=0], [x(t)-CMAX, halt]] ..."  ... I will not solve correctly the problem (my reply to JohnS).
Why ? Because the event [x(t)=0, v(t)=0] : 

  1. will preclude any initial movement of the mass as it is trivially verified by the initial condition
  2. ... and even if you use some trick to "pass" this pathological case, the equation dv(t)/dt = f will make the velocity v(t) to evolve and then x(t) to become negative

 


You can also look also to this more complete mw file
ErrorWithDsolve-2.mw

The acceleration is a see-saw function.
The first pattern ranges from 0 to 4 with maximum amplitude 1.
The second pattern ranges from 4 to 8 with maximum amplitude 2

Given c and k, depending on the value of CMAX  :

  1. x(t) = CMAX before t=2: the solution stops because of the event [x(t)-CMAX, halt]
  2. x(t=2) < CMAX : the mass moves backward up to its wall position x=0 and stays still until the acceleration acc(t) (t>4) is large enough for the piston to move forward again
    Then : 
    1. either it reaches CMAX before the acceleration reverses and forces it to return to its wall position
    2. either it reaches CMAX and the solution stops because of the event [x(t)-CMAX, halt]

 

 

@JohnS 

Good remarks.

  1. In the expression of f, there are only two effects : the acceleration, as you obeserved, and a linear spring of length L0 and mounting length Lm < L0. Given the stiffness k the spring develops a resustive force k(Lm-L0) I just wrote c
    So c is not some dampind coefficient.
    Note that in the problem of interest I also have a damping term
     
  2. Why do you expect the condition x = v = 0 to occur ?These equations model the behaviour of a mass (here with value 1) contained in a box submitted toa given acceleration..
    At time 0 the mass is pushed against the bottom of the box by the spring (force of modulus c).
    If the acceleration is high enough to balance this force the mass begins to move to the opposite side of the box.
    This happens at some time Tdec__1 (at this time x(t) is obviously null)
    Two events can happen :
    1. the displacement of the mass equals CMAX and the system locks itself
    2. the acceleration diminishes and the spring forces the mass to return to its initial rest position x(t)=0
       in this case the system is no longer governed by an ODE and its state is trivialy x(t)=v(t)=0.
      This static equilibrium lasts until the acceleration increases again to balance the prestressing c.
      ​​​

So it is necessary to use events to monitor the differents situations the mass can face.

 

@acer 
I just thought it could have been a rather common situation, not case specific, and that the code/data was not necessary.


I'm sorry

First 134 135 136 137 138 139 140 Page 136 of 140