acer

31819 Reputation

29 Badges

19 years, 180 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Are you going to be looking for floating-point or exact final results? If floating-point, then do you suspect that there might be any numerical difficulties encountered when using 'solve' as an intermediate exact solver here? Do you know anything about X, such as whether it might be symmetric?

I ask all this only because there is a LAPACK routine dtgsyl which might be of interest, depending on the answers to my questions above. There might already be a somewhat nicer interface to using it in Maple 12 (ie. nicer than setting up one's own external call to it, and already handling transformation to Schur form as an initial step).

> restart:
> kernelopts(opaquemodules=false):
>
> eval(DynamicSystems:-SylvesterSolve);
proc(A, trana, schura, B, tranb, schurb, C, isgn, $)
description "solve the Sylvester equation A . X + X . B = C"
     ...
end proc
 
> eval(DynamicSystems:-`SylvesterSolve/HardwareTable`);
table([(complex[8], rectangular) = [ztrsyl_],
       (float[8], rectangular) = [dtrsyl_]
    ])

acer

You might try simplifying the expression under certain assumptions on the variables.

For example,

> expr := sqrt(a^2*p*(1-p)*(N-n)/(n^2*(N-1))):

> simplify(expr) assuming n::real, a::real, p>0, p<1, N>1;
                       1/2        1/2        1/2
                      p    (1 - p)    (N - n)    | a/n |
                      ----------------------------------
                                         1/2
                                  (N - 1)

You might also wish to consider the explanation of the 'symbolic' option of the simplify() routine, as it appears in the help-page ?simplify,details.

acer

Have a look at the help-page ?Optimization,Options and search for the feasibilitytolerance option. It can be used to pass a specific value of the allowed violation in the constraints.

> restart:

> with(Optimization):

> cnsts := {x1+x2+x3+x4+x5 >= 1,
>  x1+x2+x3+x11+x6 >= 1, x7+x9+x10+x4+x5>=1,
>  x7+x9+x10+x4+x5 >= 1, x7+x9+x10+x11+x6 >=1,
>  x8+x9+x10+x4+x5 >=1, x8+x9+x10+x11+x6 >=1}:

> obj := x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11:

> sol := LPSolve(obj,cnsts,assume={nonnegative}):

> eval(convert(cnsts,list),sol[2]);
[1 <= 1.000000000, 1 <= 0.9999999990, 1 <= 1.000000000, 1 <= 0.9999999990,
 
    1 <= 1.000000000, 1 <= 0.9999999990]
 
> sol2 := LPSolve(obj,cnsts,assume={nonnegative},
>                 feasibilitytolerance=0.1e-13):

> eval(convert(cnsts,list),sol2[2]);
[1 <= 1.000000000, 1 <= 1.000000000, 1 <= 1.000000000, 1 <= 1.000000000,
 
    1 <= 1.000000000, 1 <= 1.000000000]

I suppose that it might not always make sense to set that tolerance very small without also increasing the working precision.

You can also set infolevel[Optimization]:=3 and see which value gets used for that tolerance, even when you don't specify it yourself as an optional argument.

acer

Do you still have your original activation code? Would it be possible to try uninstalling and then reinstalling Maple? Just an idea.

acer

I figure that by "resolution" you mean Digits, the setting of the precision for floating-point calculation.

I'm not sure whether the following is useful to you,


> restart:
>
> p := Pi*csc(Pi*alp)
> *((E*x)^m*hypergeom([m],[1-alp,1+m],E*x)
> /(GAMMA(c)*GAMMA(1-alp)*GAMMA(1+m))
> -(E*x)^c*hypergeom([c], [1+alp, 1+c], E*x)
> /(GAMMA(m)*GAMMA(1+alp)*GAMMA(1+c)));
                    /     m
                    |(E x)  hypergeom([m], [1 + m, 1 - alp], E x)
p := Pi csc(Pi alp) |--------------------------------------------
                    \    GAMMA(c) GAMMA(1 - alp) GAMMA(1 + m)

            c                                      \
       (E x)  hypergeom([c], [1 + c, 1 + alp], E x)|
     - --------------------------------------------|
           GAMMA(m) GAMMA(1 + alp) GAMMA(1 + c)    /

>
> thisp := eval(subs(alp=c-m,p),[c=5,m=1/2,E=25]);
            /     1/2  1/2
            |35 25    x    hypergeom([1/2], [-7/2, 3/2], 25 x)
thisp := Pi |-- ----------------------------------------------
            \64                       Pi

                5                                \
       1562500 x  hypergeom([5], [11/2, 6], 25 x)|
     - ------- ----------------------------------|
         567                   Pi                /

>
> evalf[10](eval(thisp,x=40));
                                              23
                              -0.3141592654 10

> evalf[30](eval(thisp,x=40));
                                      0.

> evalf[50](eval(thisp,x=40));
             0.99999999999999999831055348863866510505212220867081

>
> newthisp:=expand(map(convert,convert(thisp,StandardFunctions),expln)):

> newthisp := simplify(newthisp) assuming x>0;

newthisp := -1/192

           (3/2)         1/2         2                              1/2
    (7000 x      + 1395 x    + 5000 x  + 4350 x + 192 - 192 exp(10 x   ))

             1/2
    exp(-10 x   )

>
> evalf[10](eval(newthisp,x=40));
                                 0.9999999999

> evalf[30](eval(newthisp,x=40));
                       0.999999999999999999999982319406

The following form may make it appear more clear,

> collect(expand(newthisp),exp(x^(1/2))^10);
                           3/2        1/2        2
                      875 x      465 x      625 x    725 x
                    - -------- - -------- - ------ - ----- - 1
                         24         64        24      32
                1 + ------------------------------------------
                                        1/2 10
                                   exp(x   )

A loftier goal might be to get a nice form without having to instantiate at specific values of c, m, and E. Maybe it can be done under assumptions, but I do not yet see how.

acer

I get this too, with Maple 12 on an (old, FC2) Fedora Core Linux system. But I don't get it on a newer SuSE 10.1 machine.

acer

> restart:

> Units:-AddUnit(roung,prefix=SI,abbreviation=rg,context=SI,conversion=m/s^2);
> Units:-AddSystem(ChrisSI,Units:-GetSystem(SI),roung);
> Units:-UseSystem(ChrisSI);

> with(Units:-Standard):

> 5*Unit(watt*second)/Unit(m*kg);
                                    5 [rg]
 
> 50*Unit(krg)*Unit(s)/Unit(m);
                                  50000 [1/s]

Do you mean something other than the two units palettes? If you configure your left-bar's palettes to contain either (one for SI, and one for FPS), then you should see a generic "unit" entry which allows entry of any valid unit. That is one way to get the bracketed units in 2D Math input, say.

acer

The Units package is a module. It has both exports (such as UseSystem) and two submodules (Standard and Natural).

The exported routines are for administrating the package (redefining units, switching systems, etc). But it is the submodules which provide environments in which units may be used and recognized by (some) of Maple.

The Standard subpackage is one in which units are entered using the Unit() constructor routine. The Natural subpackage provides an environment in which units are stored internally as Units() calls, but in which units are actually entered or typed in by simply using the common symbols or names alone (eg. m, kg). The problem with Units:-Natural is that there are so very many such symbols that the global namespace is robbed of far too many common variable names.

Try initializing your session with the command with(Units:-Standard).

acer

Yes. But how to do that depends on which OS your are using.

acer

The rhs() routine expects an equation, and will not accept a set (even if its only member is an equation). So you could get your hands on that equation itself, or map the rhs() function over all the elements of the set.

> p1 := x^2+y;
                                        2
                                 p1 := x  + y

> e:=solve({p1}, {x});
                                    1/2             1/2
                      e := {x = (-y)   }, {x = -(-y)   }


> op(e[1]);
                                          1/2
                                  x = (-y)

> x1 := rhs(op(e[1]));
                                           1/2
                                 x1 := (-y)

> x1 := map(rhs,e[1]);
                                           1/2
                                x1 := {(-y)   }

> e[1][1];
                                          1/2
                                  x = (-y)

> x1 := rhs(e[1][1]);
                                           1/2
                                 x1 := (-y)

acer

> eq:=(m+M)=a^3/P^2;
                                              3
                                             a
                              eq := m + M = ----
                                              2
                                             P

> solve(algsubs(P^(-2)=Z,eq),Z);
                                     m + M
                                     -----
                                       3
                                      a

It can be very difficult to control the order of printing of terms in sums, unless one resorts to subverting the mechanism (using the `` operator, for example). I believe that it depends on the order in which the terms appear in some structure (simpl table) internal to the kernel, and that may be based on memory address and thus be session dependent.

> restart:
> eq:=(m+M)=a^3/P^2:
> expand(solve(eq,m));
                                          3
                                         a
                                   -M + ----
                                          2
                                         P

> restart:
> a^3/P^2-M:
> eq:=(m+M)=a^3/P^2:
> expand(solve(eq,m));
                                     3
                                    a
                                   ---- - M
                                     2
                                    P

acer

It's very difficult to give particular advice on what may be the cause without seeing the code. You can upload it in a file, on this site. (Use the green up-arrow FileManager icon on the editor pane's toolbar.)

acer

I'm not sure what you mean by the term "intergration"[sic] in this context..

Are you saying that you want an exact result, instead of a floating-point approximation? If so, then don't apply evalf to the result.

Are you also saying that you are seeking a purely real result, with no nonzero imaginary component? If so then read the help-page ?root and note the item describing the "principal root". Maybe you do not actually intend to use the power 1/5, which will behave according to that description.

> h := t-> (1+t^3*cos(t))^(1/5):
> res := ((D@@5)(h))(Pi):
> simplify(res);
        1/5                              3          2          10           6
- 3 (-1)    (12500 - 60000 Pi + 140000 Pi  - 3125 Pi  - 2856 Pi   - 90000 Pi
 
                4            7           9           8          5         12
     - 468000 Pi  - 244800 Pi  - 63200 Pi  + 20250 Pi  - 5000 Pi  + 700 Pi
 
               11         14    /               3 24/5
     - 13000 Pi   + 875 Pi  )  /  (3125 (-1 + Pi )    )
                              /
 
> evalf(%);
                        -0.1245723166 - 0.09050708583 I
 
> H := t-> surd( (1+t^3*cos(t)), 5 ):
> Res := ((D@@5)(H))(Pi):
> simplify(Res);
                               3          2          10           6
3 (12500 - 60000 Pi + 140000 Pi  - 3125 Pi  - 2856 Pi   - 90000 Pi
 
                4            7           9           8          5         12
     - 468000 Pi  - 244800 Pi  - 63200 Pi  + 20250 Pi  - 5000 Pi  + 700 Pi
 
               11         14    /               3 24/5
     - 13000 Pi   + 875 Pi  )  /  (3125 (-1 + Pi )    )
                              /
 
> evalf(%);
                                 0.1539798515

> evalf(-(-1)^(1/5));
                        -0.8090169944 - 0.5877852523 I
 
> % * %%;
                       -0.1245723167 - 0.09050708586 I

So, the second result, Res, is purely real, and exact.

Notice that the exact results res and Res differ by a factor of -(-1)^(1/5), and that Maple takes that using the principal root concept mentioned in that help-page.

Do you truly intent to convey that the answer should be an exactly equal to some rational number, and if so, why? Or did you simply intead "real", instead of "rational"?

acer

The T that you used in your decAdd operator is not the same as the T being used by G.

Try using G:-variable instead.

Note also that what you are passing to eval() may look (get printed) as a polynomial but underneath is a modp1(ConvertIn(..)) object. So it won't work properly, using eval that way, even with G:-variable. That's why Alec's suggestion to use G:-output is better.

acer

You might consider using option overload. See ?overload

acer

First 307 308 309 310 311 312 313 Last Page 309 of 330