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

You could try,

seq([eval(t,1),length(t)], t in anames(':-user'));

or make your own more complicated version. See ?length

acer

You can apply both those rules, together or separately. Below they are applied separately. Is that what you mean?

restart:

ee:=cos(3*x)*cos(x) + sin(3*x)*sin(x):

applyrule(cos(a::integer*x)*cos(b::integer*x)
          =1/2*(cos((a-b)*x)+cos((a+b)*x)),
          ee);

           1            1                           
           - cos(2 x) + - cos(4 x) + sin(3 x) sin(x)
           2            2                           

applyrule(sin(a::integer*x)*sin(b::integer*x)
          =1/2*(cos((a-b)*x)-cos((a+b)*x)),
          ee);

                             1            1         
           cos(3 x) cos(x) + - cos(2 x) - - cos(4 x)
                             2            2         

If both are applied, then the terms collapse, and produce just cos(2*x).

acer

A := Matrix([[a, b], [b, c]], shape = symmetric, attributes = [positive_definite]):

B := LinearAlgebra:-LUDecomposition(A, method = Cholesky);
Bhat := simplify(B) assuming real, positive;
map(normal, Bhat.Bhat^%T - A );

F := LinearAlgebra:-LUDecomposition(A, method = Cholesky, conjugate=false);
map(normal, F.F^%T - A );

If you go with conjugate=false then you might still want to simplify under the assumption of being positive, depending on example.

acer

In 16.01 on 64bit Linux on an Intel i5, after about 5 minutes,

SolveTools:-SemiAlgebraic({f1+f2 = 5, f3+f4 = 11/2,
 0 <f1, 0 < f2, 0 < f3, 0 < f4,
 6 < f2+f3, 6 < f2+f4,
 7 < f1+f4, 8 < f1+f3},
 [f1,f2,f3,f4]);

 []

which, with [] as output, means no solutions.

More quickly,

SolveTools:-SemiAlgebraic({ a + b = 1, a > 0, b > 1, c > 0 }, {a,b,c});

 []

However, the computation time grows quickly with problem size, for this approach.

The `solve` command is using an older, faster, and buggier algorithm here.

acer

Instead of using a table (or list) and complexplot you could instead use a hardware datatype Matrix and call plots:-pointplot.

First, create a Matrix M like so,

   M := Matrix( 10^6, 2, datatype=float[8] ):

Then populate the first column of M with the real values and the second column with the imaginary values. If possible, alter your original code to do this directly instead of producing a table (or list) of complex values and copying values over into M by using `Re` and `Im`.

Then you could plot it with something like,

   plots:-pointplot( M );

If done right then this should be able to display the graphic of the plot in a matter of seconds rather than hours.

Your statement that you don't want an approximation is unclear. All plots turn any exact data into floating-point numbers at some point before rendering.

acer