Here are two approaches, the second of which gets that precise form.

(There are more ways to do this particular example, but from just this one stated example and no more description it is not possible to completely characterize some general task that you may or may not have in mind.)

Download denom.mw

acer

It's not clear what you want, exactly. But you may be able to adjust this to suit.

restart:

M:=Matrix(10, 10,
 {(1, 1) = 1, (1, 2) = 1, (1, 3) = 1, (1, 4) = 1, (1, 5) = 1, (1, 6) = 1,
 (1, 7) = 1, (1, 8) = 1, (1, 9) = 1, (1, 10) = 1, (2, 1) = 1, (2, 2) = 1,
 (2, 3) = 2, (2, 4) = 2, (2, 5) = 2, (2, 6) = 2, (2, 7) = 1, (2, 8) = 1,
 (2, 9) = 1, (2, 10) = 1, (3, 1) = 1, (3, 2) = 2, (3, 3) = 2, (3, 4) = 2,
 (3, 5) = 2, (3, 6) = 2, (3, 7) = 2, (3, 8) = 2, (3, 9) = 1, (3, 10) = 1,
 (4, 1) = 1, (4, 2) = 2, (4, 3) = 2, (4, 4) = 2, (4, 5) = 2, (4, 6) = 2,
 (4, 7) = 2, (4, 8) = 2, (4, 9) = 1, (4, 10) = 1, (5, 1) = 2, (5, 2) = 2,
 (5, 3) = 2, (5, 4) = 2, (5, 5) = 2, (5, 6) = 2, (5, 7) = 2, (5, 8) = 2,
 (5, 9) = 2, (5, 10) = 1, (6, 1) = 2, (6, 2) = 2, (6, 3) = 2, (6, 4) = 2,
 (6, 5) = 2, (6, 6) = 2, (6, 7) = 2, (6, 8) = 2, (6, 9) = 2, (6, 10) = 1,
 (7, 1) = 2, (7, 2) = 2, (7, 3) = 2, (7, 4) = 2, (7, 5) = 2, (7, 6) = 2,
 (7, 7) = 2, (7, 8) = 2, (7, 9) = 2, (7, 10) = 1, (8, 1) = 2, (8, 2) = 2,
 (8, 3) = 2, (8, 4) = 2, (8, 5) = 2, (8, 6) = 2, (8, 7) = 2, (8, 8) = 2,
 (8, 9) = 2, (8, 10) = 1, (9, 1) = 2, (9, 2) = 2, (9, 3) = 2, (9, 4) = 2,
 (9, 5) = 2, (9, 6) = 2, (9, 7) = 2, (9, 8) = 2, (9, 9) = 2, (9, 10) = 1,
 (10, 1) = 2, (10, 2) = 2, (10, 3) = 2, (10, 4) = 2, (10, 5) = 2, (10, 6) = 2,
 (10, 7) = 2, (10, 8) = 2, (10, 9) = 2, (10, 10) = 1}):

m,n:=LinearAlgebra:-Dimensions(M):

Plines:=plots:-display(
  seq(seq(`if`(M[m+1-j,i-1]<>M[m+1-j,i],
               plottools:-line([i-0.5,j-1+0.5],[i-0.5,j+0.5]),NULL),
          j=1..m),i=2..n),
  seq(seq(`if`(M[m+1-j,i]<>M[m+2-j,i],
               plottools:-line([i-0.5,j-1+0.5],[i-0.5+1,j-1+0.5]),NULL),
          j=2..m),i=1..n),
              view=[1..m,1..n]):

Ppoints:=plots:-display(
  plots:-sparsematrixplot(Matrix(m,n,(i,j)->M[m+1-i,j]-1),color=red),
  plots:-sparsematrixplot(Matrix(m,n,(i,j)->M[m+1-i,j]-2),color=green)):

plots:-display( Plines, Ppoints );

acer

In your example, `c` has been assigned the value that you see. The fact that the value came from evaluating a+b is lost (except mathematically, but you don't seem to be wanting that), by which I mean that in your example the literal 'a+b' is not a part of the value assigned to name `c`.

If you reverse the order of operations, then yes you can evaluate `c` to one level, and get the unevaluated 'a+b'.

> restart:

> c:=a+b;

                                  c := a + b

> a:=5*x^2-3*y^3:
 
> b:=3*x^2+10*y^3:

> c;

                                     2      3
                                  8 x  + 7 y

> eval(c,1);

                                     a + b

In my example above it is full evaluation of `c` that produces and displays that expression in x and y.

Another variation is,

> restart:

> a:=5*x^2-3*y^3:      
> b:=3*x^2+10*y^3:     

> c:='a+b';

                                  c := a + b

> c;

                                     2      3
                                  8 x  + 7 y

> eval(c,1);

                                     a + b

You could also study the evaluation rules within procedures, which differ from those for the "top level." The manuals should help here.

You might also consider the somewhat related topic of command history.

acer

One can programmatically access the body of a procedure using the `eval` command. Eg. to subs into it, or get its inert form, or assign a copy of it to something else, etc.

By default one can also use just the `print` command to display the body of a procedure. This differs from `eval` in that it only prints. The return value from this operation is NULL, so it doesn't give programmatic access. (See also ?verboseproc)

To display line numbers, one can use the `showstat` command.

> restart:

> u:=x->x^2;
                      
                            2
                      x -> x 

> f:=eval(u);

                            2
                      x -> x 

> f(3);

                         9

> print(u);

                            2
                      x -> x

> showstat(u);

u := proc(x)
   1   x^2
end proc

Note that all of these will display without showing any comments present in the original source.

acer

I don't really understand why you wouldn't post your full example in the original question. All these followups with additional examples just makes this thread hard to read (or organize).

This is as far as I got with your 3rd example. You haven't stated what assumptions you might be prepared to accept on the variables.

Trying to combine or simplify under simultaneous pair of assumptions that numerator and denominator of the final results radicals were both positive just seemed to consume time and memory, and I interrupted it.

restart:
ex:=-(1/6)*(4*e+a^2*g*2^(1/2))^(1/2)*2^(1/4)*(-64*EllipticF((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2),

 I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))*e^4

+2*EllipticF((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2))*g^3*a^6*2^(1/2)*e+g^4*a^8*EllipticE((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2),

 I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))

+8*EllipticF((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2))*e^2*g^2*a^4+2*g^3*a^6*2^(1/2)*e*EllipticE((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2)

, I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))

-16*g*a^2*2^(1/2)*e^3*EllipticE((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2

+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4

+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))-16*EllipticF((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2),

 I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))*g

*a^2*2^(1/2)*e^3-8*g^2*a^4*EllipticE((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2

+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4+4*g

*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))*e^2)/((g*a^2+2*2^(1/2)

*e)^2*g^(1/2)*(8*e^2-g^2*a^4)^(1/2)):

Q:=combine(ex,radical) assuming -g*a^2+2*2^(1/2)*e > 0:
simplify(Q,size) assuming 8*e^2-g^2*a^2 > 0:
simplify(%) assuming 8*e^2-g^2*a^2 > 0:
simplify(%,size);

                              /       /                         
                              |       |                         
              1               | (1/4) |/   2      (1/2)  \    2 
----------------------------- |2      |\g a  + 2 2      e/ g a  
                            2 |       |                         
   (1/2) /   2      (1/2)  \  |       |                         
6 g      \g a  + 2 2      e/  \       \                         

           /                                  (1/2)\     
           |  /   2    2  4        2  (1/2)  \     |     
           |I \8 e  + g  a  + 4 g a  2      e/     |     
  EllipticE|---------------------------------------| + 2 
           |                        (1/2)          |     
           |          /   2    2  4\               |     
           \          \8 e  - g  a /               /     

           /                                  (1/2)\       
           |  /   2    2  4        2  (1/2)  \     |       
           |I \8 e  + g  a  + 4 g a  2      e/     |   /   
  EllipticK|---------------------------------------| e \4 e
           |                        (1/2)          |       
           |          /   2    2  4\               |       
           \          \8 e  - g  a /               /       

                 \                                             \
                 |                    (1/2)               (1/2)|
      2    (1/2)\| /       2    (1/2)\      /   2    2  4\     |
   + a  g 2     /| \4 e + a  g 2     /      \8 e  - g  a /     |
                 |                                             |
                 |                                             |
                 /                                             /

acer

restart:

ex:= (g*a^2+2*2^(1/2)*E^(1/2))^(1/2)
     *(-g*a^2+2*2^(1/2)*E^(1/2))^(1/2)/(8*E-g^2*a^4)^(1/2):

simplify(combine(ex,radical)) assuming 2*2^(1/2)*E^(1/2) > g*a^2;

                               1

simplify(combine(ex,radical)) assuming 2*2^(1/2)*E^(1/2) < g*a^2;

                               1

simplify(combine(ex,radical)) assuming 2*2^(1/2)*E^(1/2)=g*a^2;

                               0

[edit] What happens to the denominator of `ex` under that equality assumption may be worth considering.

acer

If you are going to use the name `x` for the piecewise as well as for the name in subsequent evaluation and plotting then you could pass that name into the procedure as an argument.

As you originally had it (as commented code, apparently with `posi4` declaring `x` as a local) the `x` in the piecewise would be an escaped local and would not be the same name as the `x` in your later `plot` call.

Also, the business of using `assign` on the result from `solve` can get awkward. In the attached sheet I use 2-argument `eval` instead.

posi1_C_modif.mw

This revision to the methodology works both with and without a wrapping procedure. (Here it is without a proc.)

acer

Quite often a good way to do this is to set it up with dsolve/numeric's events functionality. The solver itself is in charge of tracking its own tolerances internally, and adjusting stepping dynamically, and letting it do that with the goal of noticing a particular event (such as hitting a zero) makes a lot of sense. Slightly less good (sometimes, b/c the accuracy of the procs emitted by dsolve/numeric are only as good as the piecewise polynomial interpolants that get set up) is to do rootfinding (fsolve, DirectSearch, etc) on option listprocedure output. And interpolating an array of data points (output from dsolve/numeric at steps chosen by the programmer) can be worse, since that data can lack any finer grained stepping that dsolve/numerics solver might have had to do near a root.

By the way, you should be able to find several older posts that relate; just use the search field at the top of this site's pages. (eg. searching for events gets here, here, etc.)

acer

If you're willing to write out the `legend` entries by hand (or do a little coding to manipulate given values of the parameter), you might try it like,

  legend=[typeset(Beta[0]=10^(-``*9))]

That should give you 9 in an upright font. Using 10^(-`9`) would show the 9 in italics, a mismatch with the upright 10. And 10^``(-9) would display with surrounding brackets.

acer

If I understand the issue rightly, then you might be able to use fsolve's `avoid` option to find (potentially) up to all three roots. And then you could just pick whichever of those had the least `t` value.

I'm not sure that I understand why Digits=180 is necessary.

The method I've outlined will have to be allowed to fail out, by which I mean that at every iteration if will have to be allowed to search until it hits its own internal number-of-attempts-limit. This means that it will run slower as a whole.

When using a type-check to test the return value from fsolve computed at the top-level, you should use eval(soln,1), If you don't then in the case that fsolve failed the returned unevaluated `fsolve` call will get re-evaluated by virtue of its being an argument in the call to `type`.

billiard_simulatio_m.mw

Let me know if I've misunderstood.

acer

Here's one way, without changing your defns for `dist` and `v`. (This is 1D input, but you should be able to do it in 2D input as well.)

restart:

v:=3.9*Unit(km/h):
dist:=t->v*t;

plot( combine(dist(t*Unit(s)),units), t = 0*Unit(s)..3*Unit(s));

acer

You've forgotten the colon in := when trying to assign to `Rounding`.

restart;

ans:=fsolve(sin(x)=0, x=2..4);

                          3.141592654

evalf[5](ans);

                             3.1416

Rounding:=-infinity:

evalf[5](ans);

                             3.1415

acer

Your subscripted thetac is a table reference. You sheet is doing the 2D Math input equivalent of this 1D input,

`θc`:=1:

`θc`;

                               1

`θc`[deg]:=4:

`θc`;
                            θc

By assigning to the table entry you clobber the previously assigned value of the base name. You can get around this by making the 2D Math input of the subscripted thetac[deg] instead be a unique name (so that it does not collide with base name thetac). To do this, highlight all of your subscripted thetac[deg] input, on the left hand side of that assignment statement, and use the right-click context-menu action 2-D Math -> Convert To -> Atomic Identifier.

If you do this then you'll have to make sure that you also use the same atomic identifier whenever you use it or refer to it later. You could repeat the context-menu action each time, or cut & paste it from earlier.

acer

This is a cross-platform regression in the Standard GUI's 2D plot renderer, introduced it seems in 16.02.

You could try reinstalling 16.01, if you are able. Or you could try Preben's CURVE splitting fix-up procedure.

It's been reported several times, in various plotting commands, since December. (eg, 1, 2, 3, 4, 5)

acer

It's not clear, right now, that you know for sure that there is in fact a root with a1:=1.514 and a2:=5.049. Perhaps you are willing to allow a small violation of the equations, in which case those could be constraints of an optimization problem for which you'd likely need a global solver.

It might be interesting to see how the other variables change, w.r.t. a1 and a2, near the previous full precision root you cited.

In your earlier thread, Carl Love mentioned controlling accuracy and precision separately. I too had been working along those lines, to investigate. I did not find a root, but it's not clear whether that is because the ranges are so great, or the functions are so nonlinear, or perhaps there is no root for those chopped a1 and a2 values.

Anyway, here is one way to have the evaluations done at higher precision while still allowing fsolve to do its usual thing and compute an accuracy tolerance based on the inbound Digits value.

restart:

read "Cvector_with_a.m":

Cnew:=eval(C,[a1=1.514,a2=5.049]):

for i from 1 to 9 do
   Cproc[i]:=codegen[optimize](
                unapply(Cnew[i],[indets(Cnew[5],name)[]])):
end do:

save Cproc, "Cproc.m":  # you only ever want to do this once

And then,

restart:

read "Cproc.m":

for i from 1 to 9 do
   f[i]:=subs(counter=i,proc(cc,C1,C2,C3,C4,C5,C6,C7,HB)
      Digits:=100:
      Cproc[counter](cc,C1,C2,C3,C4,C5,C6,C7,HB);
   end proc);
end do:

Cproc[5](1.0,
         246.851, 9.597, 7.153, 56.14, 18.91, 2.656, 7.73,
         0.012);

                                              11
                              -0.9931016321 10

f[5](1.0,
     246.851, 9.597, 7.153, 56.14, 18.91, 2.656, 7.73,
     0.012);

-0.992357302391520534067512587085549339330928664664135982173091890590934784\

                                   11
    6556567851795358311687275614 10

and so on. Of course not all those decimal digits in the result from f[5](...) above are correct. But hopefully there are enough correct (more than from Cproc[5] called at default working precision, anyway) to allow a numerical solver to get a grip and converge if it can.

The procedures f[i] use the codegen[optimized] procedures's formulas -- by virtue of calling Cproc[i] -- but they each temporarily raise Digits to 100 before making those calls.

You can then try feeding those to a rootfinder (or global optimization scheme, if that's what you want). And you could replace the hard-coded 100 by some name, so that the higher working precision was under your control.

acer