You could, probably, use the surfdata command instead of pointplot3d. It produces a smooth surface from a rather small number of points.

Alec

The second equation can be solved for q1, and the result can be substituted in eq1, which gives one equation with one variable q2.

solve(eq2,q1);
                             5.306666667
                       -----------------------
                       sinh(1.740000000 q2) q2

eval(eq1,q1=%);

                2.653333334 cosh(1.740000000 q2)
  0.760044444 + --------------------------------
                    sinh(1.740000000 q2) q2

                  1.760044444
         - ------------------------- = 0
                               2   2
           sinh(1.740000000 q2)  q2

That can be used to find some roots,

numer(combine(lhs(%)));
                                     2
  1.140066666 sinh(1.740000000 q2) q2

                          2
         - 0.3800222220 q2  sinh(5.220000000 q2)

         + 1.326666667 cosh(1.740000000 q2) q2

         - 1.326666667 q2 cosh(5.220000000 q2)

         + 3.520088888 sinh(1.740000000 q2)

RootFinding:-Analytic(%,q2=-10-10*I..10+10*I);

  0. + 0. I, -5.3549598782850 I, -4.0974193995796 I,

        -5.4165390579135 I, -7.5157023807390 I, -3.6110260386090 I,

        -7.2220520772180 I, -3.5164513023701 I, -5.7855081270120 I,

        -2.4881678566234 I, -1.8055130193046 I, -7.1762291027045 I,

        -1.5797888398084 I, -8.9910373906423 I, -1.1037678541635 I,

        -9.0275650965225 I, -0.47765517662736, -9.2699820114995 I,

        0.47765517662736, 5.3549598782855 I, 5.4165390579130 I,

        1.8055130193045 I, 5.7855081270120 I, 4.0974193995796 I,

        7.1762291027045 I, 3.6110260386090 I, 7.2220520772180 I,

        3.5164513023701 I, 7.5157023807390 I, 2.4881678566234 I,

        8.9910373906467 I, 1.5797888398084 I, 9.0275650965225 I,

        1.1037678541635 I, 9.2699820114990 I

Zero should be excluded (it contradicts to eq2), the 2 real solutions are the same that we already found using fsolve, and one can see that there are many imaginary solutions.

Alec

You could try fsolve,

fsolve({eq1, eq2});

                {q1 = 11.94394795, q2 = 0.4776551767}

That doesn't give all solutions though. In particular, the same q1 and -q2 instead of q2 is also a solution.

fsolve({eq1,eq2},{q1=11,q2=-0.5});

                {q1 = 11.94394795, q2 = -0.4776551767}

Alec

Add a decimal dot somewhere - for example, 22. instead of 22

e1:=22. - 44*t + 216*t^2 - 66*u - 450*u^2 + 362*u^3;
e2:=150 - 300*t + 185*t^2 - 450*u + 795*u^2 - 495*u^3;

solve({e1,e2});

  {t = 0.3553932236, u = 0.2256393210},

        {t = -0.3927202492, u = 1.261618450}, {

        t = 1.089017838 + 0.1538967327 I,

        u = 0.7930027610 + 0.4135809565 I}, {

        t = 0.05432219404 - 0.6134738525 I,

        u = -0.07035986551 + 0.3340858801 I}, {

        t = 0.05432219404 + 0.6134738525 I,

        u = -0.07035986551 - 0.3340858801 I}, {

        t = 1.089017838 - 0.1538967327 I,

        u = 0.7930027610 - 0.4135809565 I}

Alec

An obvious suggestion is to plot less number of points - few thousands, maybe, instead of 1 million. You won't be able to see them anyway, so what is the point in plotting all of them?

Alec

Just don't assign it. Type something like

plot(x^2,x=-1..1);

Hit enter and it will be displayed.

Alec

Well, at least for ww there is no any difference, and the piecewise is not needed, because matrices A.B and B.A have the same trace for any A and B, and ww can be written as

ww:=unapply(LinearAlgebra:-Trace(M1(g).M2(z)),[z,g]);

                   ww := (z, g) -> 2 + (-2 + 2 g) z

which is a suitable form for the Compiler:-Compile.

Two expressions in the piecewise for ss have different signs, with a discontinuity at z=0.5.

Also, I wouldn't call the current form of w "readable".

Alec

One of the best ways is to correct the code so that Compiler:-Compile can be used. Some simplifications like the one above, also can be tried.

You use the piecewise in a wrong format. The correct way is something like

piecewise(x<0.195, an expression, x<0.205, another expression, x<0.5, other expression, final expression)

Perhaps, with <= instead of <, if that works better. What would be definitely better, is to get rid of piecewise and use functions w_1 on one interval, w_2 on another one, w_3 on the third one etc.

Matrices in the form similar to given originally in w, should be avoided.

Alec

It seems as if all those piecewise and Matrices are not necessary,

simplify(w(z,g)) assuming z<.195,g<1,g>0;
                                       2                2    1/2
     0.001038783189 | -1. + 2. z - 2. z  - 2. z g + 2. z  g |
     -----------------------------------------------------------
                                     2 1/4
                           (1. - 1. g )

simplify(w(z,g)) assuming z>.195,z<.205,g<1,g>0;
                                       2                2    1/2
     0.001038783189 | -1. + 2. z - 2. z  - 2. z g + 2. z  g |
     -----------------------------------------------------------
                                     2 1/4
                           (1. - 1. g )

simplify(w(z,g)) assuming z>.205,z<0.5,g<1,g>0;
                                       2                2    1/2
     0.001038783189 | -1. + 2. z - 2. z  - 2. z g + 2. z  g |
     -----------------------------------------------------------
                                     2 1/4
                           (1. - 1. g )

simplify(w(z,g)) assuming z>.5,z<1,g<1,g>0;
                                       2                2    1/2
     0.001038783189 | -1. + 2. z - 2. z  - 2. z g + 2. z  g |
     -----------------------------------------------------------
                                     2 1/4
                           (1. - 1. g )

All answers are the same, so w can be written as

W:=unapply(%,[z,g]);

  W := (z, g) ->

                                          2                2    1/2
        0.001038783189 | -1. + 2. z - 2. z  - 2. z g + 2. z  g |
        -----------------------------------------------------------
                                        2 1/4
                              (1. - 1. g )

Alec

A procedure for calculating the date after adding n months can be also written as

`&+`:=proc(s::string,n::posint)
    local m,d,y,ys,y_new;
    m,d,y:=map(parse,StringTools:-Split(s,"/"))[];
    m:=irem(m+n, 12, ys);
    y_new:=y+ys;
    if y<100 and y_new>=100 then y:=y_new+1900 
        else y:=y_new fi;
    sprintf("%02d/%02d/%02d",m,d,y)
end:

For example,

"09/03/73" &+ 6;

                              "03/03/74"

"08/25/99" &+ 6;

                             "02/25/2000"

The number of days between 2 dates can be calculated by dividing the difference between their epochs by 60*60*24.

The epoch can be found, for example, using a procedure mktime from my blog.

For example, in Linux, using the recent version with WRAPPER,

`&-` := proc (a::string, b::string) 
    local a_m, a_d, a_y, a_T, b_m, b_d, b_y, b_T; 
    global mktime; 
    uses StringTools;
    a_m, a_d, a_y := map(parse, Split(a, "/"))[]; 
    b_m, b_d, b_y := map(parse, Split(b, "/"))[]; 
    a_T := table([tm_mon = a_m-1, tm_mday = a_d, 
        tm_year = `if`(a_y < 100, a_y, a_y-1900), 
        tm_sec = 1, tm_min = 0, tm_hour = 0, tm_wday = 0, tm_yday = 0, 
        tm_isdst = 0, tm_gmtoff = 0, tm_zone = ""]); 
    b_T := table([tm_mon = b_m-1, tm_mday = b_d, 
        tm_year = `if`(b_y < 100, b_y, b_y-1900), 
        tm_sec = 1, tm_min = 0, tm_hour = 0, tm_wday = 0, tm_yday = 0, 
        tm_isdst = 0, tm_gmtoff = 0, tm_zone = ""]); 
    (mktime(a_T)-mktime(b_T))/86400 
end: 

For example,

"03/03/74" &- "09/03/73";

                                     181

"02/25/2000" &- "08/25/99";

                                     184

Alec

It looks as if you also could use hasmset instead of lists, from the very beginning. Something like the following,

use hashmset in 
    a:=new();
    b:=new();
    insert(2,2,a),insert(3,3,a),insert(5,1,a);
    insert(2,3,b),insert(3,1,b),insert(7,1,b);
    entries(a intersect b);
    map(`$`@op,%);
end;
                        a := table(sparse, [])
                        b := table(sparse, [])
                               2, 3, 1
                               3, 1, 1
                           [[2, 2], [3, 1]]
                              [2, 2, 3]

Alec

You could use evalf with more digits,

evalf[20](Int(w(z,0.3),z=0..0.2));

                   0.00019925639830378505531 + 0. I

evalf[20](Int(w(z,0.3),z=0..1));

                   0.00092916665349652607360 + 0. I

It can be plotted as

plot('Re(evalf[20](Int(w(z,g),z=0..1)))',g=0..1);

Alec

First, you don't have to use that mapping. It is much easier to use the Groebner:-Basis without it, defining in your example x as RootOf(z^2+z+1).

Second, any other mapping can be used (preferably with 0 being mapped to 0). This one has an advantage of being easy to calculate - substitute x=2 going in one direction, and using base 2 representation of a number, going in another direction.

Alec

assign(%);

Alec

Mario,

You are welcome. You have a gift for starting interesting threads.

Alec