tomleslie

12497 Reputation

19 Badges

12 years, 231 days

MaplePrimes Activity


These are answers submitted by tomleslie

because

  1. implicitplot() finds three solutions for Eq1
  2. implicitplot() finds no solutions for for Eq2 - not sure why this failed, because
  3. working out solutions for Eq2 over a fixed number of points using fsolve(), one can obtain a solution - which appears to correspond with one of the solutions for Eq1(?)

See the attached (NB plots "look better" in the worksheet than they do on this site - honest!)

restart;  
Digits:=16:
Eqmin:=-deltac*((kc-1/2)/(kc-3/2))+mu*deltab/((M-u0b)**2)+1/(M**2-3*sigmai):

Eqmax:=deltac*(1-(1-(1/2)*(M-sqrt(3)*sqrt(sigmai))^2/(kc-3/2))^(-kc+3/2))+deltab*(M-u0b)^2*(1-sqrt(1+mu*(M-sqrt(3)*sqrt(sigmai))^2/(M-u0b)^2))/mu-(1/18)*sqrt(3)*(((M+sqrt(3)*sqrt(sigmai))^2-(M-sqrt(3)*sqrt(sigmai))^2)^(3/2)-(M+sqrt(3)*sqrt(sigmai))^3+(M-sqrt(3)*sqrt(sigmai))^3)/sqrt(sigmai)*(sqrt(3)*(((M+sqrt(3)*sqrt(sigmai))^2-(M-sqrt(3)*sqrt(sigmai))^2)^(3/2)-(M+sqrt(3)*sqrt(sigmai))^3+(M-sqrt(3)*sqrt(sigmai))^3)):

params:=[kc=5,u0b=0.05,deltab=0.00025,deltac=0.9975,mu=1836]:

R:=sigmai=0.01..0.7;

Eq1:=eval(Eqmin, params):
Eq2:=eval(Eqmax, params):
#
# Hmmmm three solutions for Eq1
#
  p1:=plots:-implicitplot(Eq1, R, M=0.1..3, color=red, thickness=4, gridrefine=4);
#
# Hmmmm no solutions for Eq1
#
  p2:=plots:-implicitplot(Eq2, R, M=0.1..3, color=red, gridrefine=4);
#
# But if one solves Eq2 "the hard way", then it has a solution
# which appears(?) to coincide with one of the solutions for Eq1
#
  p3:=plot([seq( [i, fsolve(eval(Eq2, sigmai=i),M=0.1..3)], i=0.01..0.7,0.01)], color=black);
  plots:-display( [ p1, p3],
                  legend=["Eq1", "Eq2"]
                );
 

sigmai = 0.1e-1 .. .7

 

 

 

 

 

 

``

NULL

NULL

NULL

NULL

Download oddPlot.mw

you may have stumbled across a bug - namely that the reflection() command appears to "error" when reflecting a circle through a point. I'm still investigating this, and if I decide that it really really is a bug then I will report it to MapleSoft.

Meanwhile, in the attached I think I have come up with a "workaround" for your specific example - shown in the attached

  restart;
  with(plots):
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:

  R := 7:
  point(A, [0, R]):
  line(L1, y = sqrt(3)*x + R):
  line(L2, y = -sqrt(3)*x + R):
  line(L3, y = R/3):
  intersection(B, L1, L3):
  intersection(C, L2, L3):
  triangle(ABC, [A, B, C]):
  circle(C1, [point(P1, [0, 0]), R]):
  circle(C2, [point(P2, [0, R/3 + (2*R)/9]), (2*R)/9]):

  circle(C3, [ reflection(P3, center(C2), C), radius(C2)]):
  circle(C4, [ reflection(P4, center(C2), B), radius(C2)]):

  draw( [ L1(color = blue),
          ABC(color = red, transparency = 0.5, filled = true),
          L2(color = blue),
          L3(color = blue),
          C1(color = blue, thickness = 3),
          C1(color = yellow, transparency = 0.8, filled = true),
          C2(color = green, filled = true),
          C3(color = green, filled = true),
          C4(color = green, filled = true)
        ],
        axes = normal,
        view = [-R .. R, -R .. R],
        scaling = constrained
      );

 

 

 

Download reflCir.mw

 

code returns a reference to a table which is generated by your procedure, rather than the table itself. If you want to return the tables contents, the last line in the procedure should be eval(V).

In fact you can make either method "work", just that the first of these involves one more "look-up" step when accessing values of the table, so is generally less efficient.

See the attached, where both approaches are shown.

  restart;

  A := Matrix([[2*t, 2, 3], [4, 5*t, 6], [7, 8*t, 9]]);
  B := Matrix([t, 2*t, 3*t])^%T;
  test1:=proc(n)
              local s,t,M,V;
              s:=1:
              for t from 1 to n do
                  M:=A+A^(-1);
                  V[s]:=(M.B):
                  s:=s+1:
              end do;
              V;
         end proc:

test2:=proc(n)
              local s,t,M,V;
              s:=1:
              for t from 1 to n do
                  M:=A+A^(-1);
                  V[s]:=(M.B):
                  s:=s+1:
              end do;
              eval(V);
         end proc:

Matrix(3, 3, {(1, 1) = 2*t, (1, 2) = 2, (1, 3) = 3, (2, 1) = 4, (2, 2) = 5*t, (2, 3) = 6, (3, 1) = 7, (3, 2) = 8*t, (3, 3) = 9})

 

Matrix(%id = 36893488148087005300)

(1)

  R:=test1(4);
  R[1];

R := V

 

Matrix(%id = 36893488148087007940)

(2)

  R:=test2(4);

table(%id = 36893488148261674268)

 

Matrix(%id = 36893488148086988300)

(3)

 

Download procref.mw

the rotation() command applies to individual geometric objects, so for example if 'L1' is a line, you can rotate this line to produce a new line.

However it appears that you may want to rotate a plot(?) - for which  you should use the plottools:-rotate() command

that the value of your expression does not depend very much on the value of the parameter 'k' whihc you supply (at least in the range 0.1..1, with R=10 and B=0.5).

Consequently, the roots don't vary noticeably with 'k' - these are at x=-1 and x=1, for all values of k in the range 0.1..1.

See further explanation in comments in the attached


 

restart

P := -R*(x+1)*((((-(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-3*k*(1/2))*x^4+(-(2*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-6*k-9/2)*x^2-9+(5*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-6*k)*R^2-(2*((-(1/4)*k-1/4)*x^4+((2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))+(1/2)*k-5/2)*x^2+(5*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))+17*k*(1/4)-7/4))*(-1+B)*R+(x^2+5)*(-1+B)^2*(k+1))*exp(-R*(x^2+2)/(-1+B))+((((2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-3/2)*x^4+((2*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-21/2)*x^2-15+(13*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1)))*R^2+(2*((-(1/4)*k-1/4)*x^4+((2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-k-4)*x^2+(5*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-13*k*(1/4)-37/4))*(-1+B)*R-(x^2+5)*(-1+B)^2*(k+1))*exp(-R*(x^2+5)/(-1+B))+((((2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-3/2)*x^4+((2*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-9*k*(1/2)-6)*x^2-6-(5*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-9*k)*R^2-(2*(((1/4)*k+1/4)*x^4+((2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))+5*k*(1/2)-1/2)*x^2+(5*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))+7*k*(1/4)-17/4))*(-1+B)*R+(x^2+5)*(-1+B)^2*(k+1))*exp(-3*R/(-1+B))+((-(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-3*k*(1/2))*x^4+(-(2*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-21*k*(1/2))*x^2-(13*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))-15*k)*R^2+(2*(((1/4)*k+1/4)*x^4+((2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1)*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))+4*k+1)*x^2+(5*(2*R*(1+(1/2)*x^2)*exp(2*R*(1+(1/2)*x^2)/(1-B))/(1-B)+2*R*k*(1+(1/2)*x^2)/(1-B)-exp(2*R*(1+(1/2)*x^2)/(1-B))*k-exp(2*R*(1+(1/2)*x^2)/(1-B))+k+1))*(1-B)/(2*R*(exp(2*R*(1+(1/2)*x^2)/(1-B))-1))+37*k*(1/4)+13/4))*(-1+B)*R-(x^2+5)*(-1+B)^2*(k+1))/((6*(((1/2)*R*x^2+B+R-1)*exp(-R*(x^2+2)/(-1+B))+(1/2)*R*x^2-B+R+1))*((B+3*R*(1/2)-1)*exp(-3*R/(-1+B))-B+3*R*(1/2)+1)*(x^2+2))

  P:=simplify(P);
#
# Plot P for k=0.1..1 in steps of 0.1. Can't really see
# much difference
#
  plots:-display
         ( [ seq
             ( plot
               ( eval
                 ( P,
                   [R = 10, B = 0.5, k = j]
                 ),
                 x = -2 .. 2
               ),
               j = 0.1 .. 1, 0.1
             )
          ]
        );
#
# Above plots indicate potential roots close to -1 and 1
# for all values of 'k'. Check this
#
  r0 := -2:
  for j from 0.1 by 0.1 to 1 do
      while true do
            r0 := RootFinding:-NextZero
                               ( unapply
                                 ( eval
                                   ( P,
                                     [R = 10, B = 0.5, k = j]
                                   ),
                                   x
                                 ),
                                 r0
                               );
            if   r0 = FAIL
            then break;
            else printf(" k=%3.1f, root=%10.8f\n", j, r0);
            end if;
      end do;
      r0 := -2;
      printf( "\n" );
  end do:

(1/6)*(x+1)*R^2*(-((1/2)*R*x^2+B+R-1)*(x+1)*(x-1)*exp((-2*x^2-7)*R/(-1+B))+((1/2)*R*(k+1)*x^4+((2*k+2)*R+(k-1)*(-1+B))*x^2+(2*k+2)*R+(-4*B+4)*k-2*B+2)*exp(-R*(x^2+2)/(-1+B))+(-(1/2)*R*(k+1)*x^4+((-2*k-2)*R-(k-1)*(-1+B))*x^2+(-2*k-2)*R-2*(k+2)*(-1+B))*exp(-R*(x^2+5)/(-1+B))+3*((1/2)*R*x^2+B+R-1)*((1/3)*x^2+k+2/3)*exp(-2*R*(x^2+2)/(-1+B))+(-(1/2)*R*x^2+B-R-1)*((k*x^2+2*k+3)*exp(-3*R/(-1+B))-k*x^2+k))/(((B+(3/2)*R-1)*exp(-3*R/(-1+B))-B+(3/2)*R+1)*(((1/2)*R*x^2+B+R-1)*exp(-R*(x^2+2)/(-1+B))+(1/2)*R*x^2-B+R+1)*(exp(-R*(x^2+2)/(-1+B))-1))

 

 

 k=0.1, root=-1.00000000
 k=0.1, root=1.00000000

 k=0.2, root=-1.00000000
 k=0.2, root=1.00000000

 k=0.3, root=-1.00000000
 k=0.3, root=1.00000000

 k=0.4, root=-1.00000000
 k=0.4, root=1.00000000

 k=0.5, root=-1.00000000
 k=0.5, root=1.00000000

 k=0.6, root=-1.00000000
 k=0.6, root=1.00000000

 k=0.7, root=-1.00000000
 k=0.7, root=1.00000000

 k=0.8, root=-1.00000000
 k=0.8, root=1.00000000

 k=0.9, root=-1.00000000
 k=0.9, root=1.00000000

 k=1.0, root=-1.00000000
 k=1.0, root=1.00000000
 

 

#
# Above is interesting, it appears that th value of the
# root (and the expression!) is pretty much independent
# of the value of 'k'.
#
# As a crude check, evaluate the expression at various
# x-values across -2..2. It would seem that the contribution
# of teh factor containing 'k' is vanishingly small
# compared with other terms!
#
  collect( eval( P, [R = 10, B = 0.5, x=-2]), k);
  collect( eval( P, [R = 10, B = 0.5, x=-1]), k);
  collect( eval( P, [R = 10, B = 0.5, x=0]), k);
  collect( eval( P, [R = 10, B = 0.5, x=1]), k);
  collect( eval( P, [R = 10, B = 0.5, x=2]), k);

3.448275863-0.3019486472e-25*k

 

-0.

 

1.149425287-0.9766331587e-17*k

 

-0.5652725092e-51

 

-10.34482759+0.9058459412e-25*k

(1)

 

Download kAndRoots.mw

from the help page:

timelimit(t, x)

The timelimit function evaluates the expression x, but gives up if the evaluation takes longer than the number of seconds specified by t.

Thiu will generate an error if the timelimit is exceeded.You will have to "trap" this using and appropriate try..catch construct to ensure that your loop will keep executing afetr the error.

DPI would "normally" be considered a property of the monitor on which you are displaying the plot.

Perhaps you are referring to the number of points at which the function [x, f(x)] is evaluate? If so, check the plot option numpoints. From the help page for the plot,options command (available by entering ?plot,options at the command prompt), one can read

numpoints=n
Specifies the minimum number of points to be generated.  The default is "200".
Note: plot employs an adaptive plotting scheme which automatically does more work where the function values do not lie close to a straight line.  Hence, plot often generates more than the minimum number of points.

 

In the attached I have provided code for each of the product graphs listed on

https://en.wikipedia.org/wiki/Graph_product

where the logical conditions are given in the "overview Table" of the above article.

I have written a procedure for each of the product graphs separately. These procedures only differ in the logical condition specified in the 'if' clause. The logical conditions used are lifted directly from the second column of the "Overview Table" in the above Wikipedia page.

See further comments in the attached.

  restart:
  with(GraphTheory):
  with(RandomGraphs):
  with(StringTools):
#
# A couple of random graphs for test purposes
#
  G1:= RandomGraph(6,10);
  G2:= RandomGraph(7,12);
#
# Generate all the vertices for the product graph
#
  vertexList:=[ seq
                ( seq
                  ( Join
                    ( [i, ":", j]
                    ),
                    i in convert~(Vertices(G1), string)
                  ),
                  j in convert~(Vertices(G2), string)
                )
              ]:
  coNorm:= proc( verts::list)
                 uses StringTools, GraphTheory:
                 local e1:= parse~( Split(verts[1], ":") ),
                       a1:= e1[1],
                       a2:= e1[2],
                       e2:= parse~( Split(verts[2], ":") ),
                       b1:= e2[1],
                       b2:= e2[2]:
                 if  `or`( HasEdge(G1, {a1, b1}),  HasEdge(G2, {a2, b2}))
                 then AddEdge(G3, {verts[1], verts[2]}):
                 fi;
           end proc:

  cartProd:= proc( verts::list)
                   uses StringTools, GraphTheory:
                   local e1:= parse~( Split(verts[1], ":") ),
                         a1:= e1[1],
                         a2:= e1[2],
                         e2:= parse~( Split(verts[2], ":") ),
                         b1:= e2[1],
                         b2:= e2[2]:
                   if  `or`( `and`( a1=b1, HasEdge(G2, {a2, b2})),
                             `and`( HasEdge(G1, {a1, b1}), a2=b2)
                           )
                   then AddEdge(G3, {verts[1], verts[2]}):
                   fi;
             end proc:
  tensProd:= proc( verts::list)
                   uses StringTools, GraphTheory:
                   local e1:= parse~( Split(verts[1], ":") ),
                         a1:= e1[1],
                         a2:= e1[2],
                         e2:= parse~( Split(verts[2], ":") ),
                         b1:= e2[1],
                         b2:= e2[2]:
                   if  `and`( HasEdge(G1, {a1, b1}), HasEdge(G2, {a2, b2}))
                   then AddEdge(G3, {verts[1], verts[2]}):
                   fi;
             end proc:
  lexiProd:= proc( verts::list)
                   uses StringTools, GraphTheory:
                   local e1:= parse~( Split(verts[1], ":") ),
                         a1:= e1[1],
                         a2:= e1[2],
                         e2:= parse~( Split(verts[2], ":") ),
                         b1:= e2[1],
                         b2:= e2[2]:
                   if  `or`( HasEdge(G1, {a1, b1}),
                             `and`( a1=b1, HasEdge(G2, {a2, b2}))
                           )
                   then AddEdge(G3, {verts[1], verts[2]}):
                   fi;
             end proc:
  strongProd:= proc( verts::list)
                     uses StringTools, GraphTheory:
                     local e1:= parse~( Split(verts[1], ":") ),
                           a1:= e1[1],
                           a2:= e1[2],
                           e2:= parse~( Split(verts[2], ":") ),
                           b1:= e2[1],
                           b2:= e2[2]:
                     if  `or`( `and`( a1=b1, HasEdge(G2, {a2, b2})),
                               `and`( HasEdge(G1, {a1, b1}), a2=b2),
                               `and`(HasEdge(G1, {a1, b1}), HasEdge(G2, {a2, b2}))
                             )
                     then AddEdge(G3, {verts[1], verts[2]}):
                     fi;
             end proc:
  moduProd:= proc( verts::list)
                   uses StringTools, GraphTheory:
                   local e1:= parse~( Split(verts[1], ":") ),
                         a1:= e1[1],
                         a2:= e1[2],
                         e2:= parse~( Split(verts[2], ":") ),
                         b1:= e2[1],
                         b2:= e2[2]:
                   if  `or`( `and`( HasEdge(G1, {a1, b1}), HasEdge(G2, {a2, b2})),
                             `and`( not HasEdge(G1, {a1, b1}), HasEdge(G2, {a2, b2}))
                           )
                   then AddEdge(G3, {verts[1], verts[2]}):
                   fi:
             end proc:
  homoProd:= proc( verts::list)
                   uses StringTools, GraphTheory:
                   local e1:= parse~( Split(verts[1], ":") ),
                         a1:= e1[1],
                         a2:= e1[2],
                         e2:= parse~( Split(verts[2], ":") ),
                         b1:= e2[1],
                         b2:= e2[2]:
                   if  `or`(  a1=b1,
                             `and`( HasEdge(G1, {a1, b1}), not HasEdge(G2, {a2, b2}))
                           )
                   then AddEdge(G3, {verts[1], verts[2]}):
                   fi:
             end proc:

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6], Array(1..6, {(1) = {2, 3, 4, 5}, (2) = {1, 3, 6}, (3) = {1, 2, 4}, (4) = {1, 3, 5, 6}, (5) = {1, 4, 6}, (6) = {2, 4, 5}}), `GRAPHLN/table/2`, 0)

 

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7], Array(1..7, {(1) = {3, 5, 6}, (2) = {3, 4, 7}, (3) = {1, 2, 6, 7}, (4) = {2, 6, 7}, (5) = {1, 6}, (6) = {1, 3, 4, 5, 7}, (7) = {2, 3, 4, 6}}), `GRAPHLN/table/4`, 0)

(1)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the conormal product graph G3
#
  coNorm~(combinat:-choose(Vertices(G3), 2) ):
#
# Check that the edge count in the product
# graph corresponds with that given in
# https://en.wikipedia.org/wiki/Graph_product
#
  is( numelems(Edges(G3)) = numelems( Vertices(G1))^2*numelems(Edges(G2))
                            +
                            numelems( Vertices(G2))^2*numelems(Edges(G1))
                            -
                            2*numelems(Edges(G2))*numelems(Edges(G1))
    );

true

(2)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the cartesian product graph G3. NB MAple already
# has a command for this, so this execution group is only
# here for illustration
#
  cartProd~(combinat:-choose(Vertices(G3), 2) ):
#
# Check that the edge count in the product
# graph corresponds with that given in
# https://en.wikipedia.org/wiki/Graph_product
#
  is( numelems(Edges(G3))= numelems( Edges( CartesianProduct(G1, G2)))
    );

true

(3)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the tensor product graph G3.
#
  tensProd~(combinat:-choose(Vertices(G3), 2) ):
#
# Check that the edge count in the product
# graph corresponds with that given in
# https://en.wikipedia.org/wiki/Graph_product
#
  is( numelems(Edges(G3))= 2*numelems(Edges(G1))*numelems(Edges(G2))
    );

true

(4)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the lexicographic product graph G3.
#
  lexiProd~(combinat:-choose(Vertices(G3), 2) ):
#
# Check that the edge count in the product
# graph corresponds with that given in
# https://en.wikipedia.org/wiki/Graph_product
#
  is( numelems(Edges(G3))= numelems(Vertices(G1))*numelems(Edges(G2))
                           +
                           numelems(Edges(G1))*numelems(Vertices(G2))^2
    );

true

(5)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the strong product graph G3.
#
  strongProd~(combinat:-choose(Vertices(G3), 2) ):
#
# Check that the edge count in the product
# graph corresponds with that given in
# https://en.wikipedia.org/wiki/Graph_product
#
  is( numelems(Edges(G3))= numelems(Vertices(G1))*numelems(Edges(G2))
                           +
                           numelems(Vertices(G2))*numelems(Edges(G1))
                           +
                           2*numelems(Edges(G1))*numelems(Edges(G2))
    );

true

(6)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the modular product graph G3.
#
  moduProd~(combinat:-choose(Vertices(G3), 2) ):
#
# Don't know how many edges this should have" No
# formula given in
# https://en.wikipedia.org/wiki/Graph_product
#
  G3;

GRAPHLN(undirected, unweighted, ["1 : 1", "2 : 1", "3 : 1", "4 : 1", "5 : 1", "6 : 1", "1 : 2", "2 : 2", "3 : 2", "4 : 2", "5 : 2", "6 : 2", "1 : 3", "2 : 3", "3 : 3", "4 : 3", "5 : 3", "6 : 3", "1 : 4", "2 : 4", "3 : 4", "4 : 4", "5 : 4", "6 : 4", "1 : 5", "2 : 5", "3 : 5", "4 : 5", "5 : 5", "6 : 5", "1 : 6", "2 : 6", "3 : 6", "4 : 6", "5 : 6", "6 : 6", "1 : 7", "2 : 7", "3 : 7", "4 : 7", "5 : 7", "6 : 7"], Array(1..42, {(1) = {13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}, (2) = {13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}, (3) = {13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}, (4) = {13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}, (5) = {13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}, (6) = {13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}, (7) = {13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 37, 38, 39, 40, 41, 42}, (8) = {13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 37, 38, 39, 40, 41, 42}, (9) = {13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 37, 38, 39, 40, 41, 42}, (10) = {13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 37, 38, 39, 40, 41, 42}, (11) = {13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 37, 38, 39, 40, 41, 42}, (12) = {13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 37, 38, 39, 40, 41, 42}, (13) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (14) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (15) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (16) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (17) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (18) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (19) = {7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (20) = {7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (21) = {7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (22) = {7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (23) = {7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (24) = {7, 8, 9, 10, 11, 12, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42}, (25) = {1, 2, 3, 4, 5, 6, 31, 32, 33, 34, 35, 36}, (26) = {1, 2, 3, 4, 5, 6, 31, 32, 33, 34, 35, 36}, (27) = {1, 2, 3, 4, 5, 6, 31, 32, 33, 34, 35, 36}, (28) = {1, 2, 3, 4, 5, 6, 31, 32, 33, 34, 35, 36}, (29) = {1, 2, 3, 4, 5, 6, 31, 32, 33, 34, 35, 36}, (30) = {1, 2, 3, 4, 5, 6, 31, 32, 33, 34, 35, 36}, (31) = {1, 2, 3, 4, 5, 6, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 37, 38, 39, 40, 41, 42}, (32) = {1, 2, 3, 4, 5, 6, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 37, 38, 39, 40, 41, 42}, (33) = {1, 2, 3, 4, 5, 6, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 37, 38, 39, 40, 41, 42}, (34) = {1, 2, 3, 4, 5, 6, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 37, 38, 39, 40, 41, 42}, (35) = {1, 2, 3, 4, 5, 6, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 37, 38, 39, 40, 41, 42}, (36) = {1, 2, 3, 4, 5, 6, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 37, 38, 39, 40, 41, 42}, (37) = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 31, 32, 33, 34, 35, 36}, (38) = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 31, 32, 33, 34, 35, 36}, (39) = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 31, 32, 33, 34, 35, 36}, (40) = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 31, 32, 33, 34, 35, 36}, (41) = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 31, 32, 33, 34, 35, 36}, (42) = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 31, 32, 33, 34, 35, 36}}), `GRAPHLN/table/11`, 0)

(7)

#
# Initialise the product graph
#
  G3:=Graph(vertexList):
#
# Compute the homomorphic product graph G3.
#
  homoProd~(combinat:-choose(Vertices(G3), 2) ):
#
# Don't know how many edges this should have" No
# formula given in
# https://en.wikipedia.org/wiki/Graph_product
#
  G3;

GRAPHLN(undirected, unweighted, ["1 : 1", "2 : 1", "3 : 1", "4 : 1", "5 : 1", "6 : 1", "1 : 2", "2 : 2", "3 : 2", "4 : 2", "5 : 2", "6 : 2", "1 : 3", "2 : 3", "3 : 3", "4 : 3", "5 : 3", "6 : 3", "1 : 4", "2 : 4", "3 : 4", "4 : 4", "5 : 4", "6 : 4", "1 : 5", "2 : 5", "3 : 5", "4 : 5", "5 : 5", "6 : 5", "1 : 6", "2 : 6", "3 : 6", "4 : 6", "5 : 6", "6 : 6", "1 : 7", "2 : 7", "3 : 7", "4 : 7", "5 : 7", "6 : 7"], Array(1..42, {(1) = {2, 3, 4, 5, 7, 8, 9, 10, 11, 13, 19, 20, 21, 22, 23, 25, 31, 37, 38, 39, 40, 41}, (2) = {1, 3, 6, 7, 8, 9, 12, 14, 19, 20, 21, 24, 26, 32, 37, 38, 39, 42}, (3) = {1, 2, 4, 7, 8, 9, 10, 15, 19, 20, 21, 22, 27, 33, 37, 38, 39, 40}, (4) = {1, 3, 5, 6, 7, 9, 10, 11, 12, 16, 19, 21, 22, 23, 24, 28, 34, 37, 39, 40, 41, 42}, (5) = {1, 4, 6, 7, 10, 11, 12, 17, 19, 22, 23, 24, 29, 35, 37, 40, 41, 42}, (6) = {2, 4, 5, 8, 10, 11, 12, 18, 20, 22, 23, 24, 30, 36, 38, 40, 41, 42}, (7) = {1, 2, 3, 4, 5, 8, 9, 10, 11, 13, 19, 25, 26, 27, 28, 29, 31, 32, 33, 34, 35, 37}, (8) = {1, 2, 3, 6, 7, 9, 12, 14, 20, 25, 26, 27, 30, 31, 32, 33, 36, 38}, (9) = {1, 2, 3, 4, 7, 8, 10, 15, 21, 25, 26, 27, 28, 31, 32, 33, 34, 39}, (10) = {1, 3, 4, 5, 6, 7, 9, 11, 12, 16, 22, 25, 27, 28, 29, 30, 31, 33, 34, 35, 36, 40}, (11) = {1, 4, 5, 6, 7, 10, 12, 17, 23, 25, 28, 29, 30, 31, 34, 35, 36, 41}, (12) = {2, 4, 5, 6, 8, 10, 11, 18, 24, 26, 28, 29, 30, 32, 34, 35, 36, 42}, (13) = {1, 7, 14, 15, 16, 17, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 31, 37}, (14) = {2, 8, 13, 15, 18, 19, 20, 21, 24, 25, 26, 27, 30, 32, 38}, (15) = {3, 9, 13, 14, 16, 19, 20, 21, 22, 25, 26, 27, 28, 33, 39}, (16) = {4, 10, 13, 15, 17, 18, 19, 21, 22, 23, 24, 25, 27, 28, 29, 30, 34, 40}, (17) = {5, 11, 13, 16, 18, 19, 22, 23, 24, 25, 28, 29, 30, 35, 41}, (18) = {6, 12, 14, 16, 17, 20, 22, 23, 24, 26, 28, 29, 30, 36, 42}, (19) = {1, 2, 3, 4, 5, 7, 13, 14, 15, 16, 17, 20, 21, 22, 23, 25, 26, 27, 28, 29, 31, 37}, (20) = {1, 2, 3, 6, 8, 13, 14, 15, 18, 19, 21, 24, 25, 26, 27, 30, 32, 38}, (21) = {1, 2, 3, 4, 9, 13, 14, 15, 16, 19, 20, 22, 25, 26, 27, 28, 33, 39}, (22) = {1, 3, 4, 5, 6, 10, 13, 15, 16, 17, 18, 19, 21, 23, 24, 25, 27, 28, 29, 30, 34, 40}, (23) = {1, 4, 5, 6, 11, 13, 16, 17, 18, 19, 22, 24, 25, 28, 29, 30, 35, 41}, (24) = {2, 4, 5, 6, 12, 14, 16, 17, 18, 20, 22, 23, 26, 28, 29, 30, 36, 42}, (25) = {1, 7, 8, 9, 10, 11, 13, 14, 15, 16, 17, 19, 20, 21, 22, 23, 26, 27, 28, 29, 31, 37, 38, 39, 40, 41}, (26) = {2, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 24, 25, 27, 30, 32, 37, 38, 39, 42}, (27) = {3, 7, 8, 9, 10, 13, 14, 15, 16, 19, 20, 21, 22, 25, 26, 28, 33, 37, 38, 39, 40}, (28) = {4, 7, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 27, 29, 30, 34, 37, 39, 40, 41, 42}, (29) = {5, 7, 10, 11, 12, 13, 16, 17, 18, 19, 22, 23, 24, 25, 28, 30, 35, 37, 40, 41, 42}, (30) = {6, 8, 10, 11, 12, 14, 16, 17, 18, 20, 22, 23, 24, 26, 28, 29, 36, 38, 40, 41, 42}, (31) = {1, 7, 8, 9, 10, 11, 13, 19, 25, 32, 33, 34, 35, 37}, (32) = {2, 7, 8, 9, 12, 14, 20, 26, 31, 33, 36, 38}, (33) = {3, 7, 8, 9, 10, 15, 21, 27, 31, 32, 34, 39}, (34) = {4, 7, 9, 10, 11, 12, 16, 22, 28, 31, 33, 35, 36, 40}, (35) = {5, 7, 10, 11, 12, 17, 23, 29, 31, 34, 36, 41}, (36) = {6, 8, 10, 11, 12, 18, 24, 30, 32, 34, 35, 42}, (37) = {1, 2, 3, 4, 5, 7, 13, 19, 25, 26, 27, 28, 29, 31, 38, 39, 40, 41}, (38) = {1, 2, 3, 6, 8, 14, 20, 25, 26, 27, 30, 32, 37, 39, 42}, (39) = {1, 2, 3, 4, 9, 15, 21, 25, 26, 27, 28, 33, 37, 38, 40}, (40) = {1, 3, 4, 5, 6, 10, 16, 22, 25, 27, 28, 29, 30, 34, 37, 39, 41, 42}, (41) = {1, 4, 5, 6, 11, 17, 23, 25, 28, 29, 30, 35, 37, 40, 42}, (42) = {2, 4, 5, 6, 12, 18, 24, 26, 28, 29, 30, 36, 38, 40, 41}}), `GRAPHLN/table/12`, 0)

(8)

 

Download productGraphs.mw

One of the easiest is (probably) to use the plotttools:-transform() command to interchange 'x' and 'y', as shown in the attached.

restart

NULL

y := 3*x^2-4*x+5

3*x^2-4*x+5

(1)

p1 := plot(y, x, axes = box, labels = ["x", "y"])

 

t1 := plottools:-transform(proc (x, y) options operator, arrow; [y, x] end proc); plots:-display(t1(p1))

 

 

Download transplot.mw

?Deep Learning.

I have never attempted to use this package, and have no knowledge of Google's "Tensorflow" software, so I have no idea how useful this might be. So far as I can tell, the DeepLearning() package is essentially a Maple "interface" to Google TensorFlow software, so I suspect that one would haave to be reasonably knowledgeable about the latter, before using the former

that by

  1. "long integer representing the ticks in JS", you mean an integer which represents time in units of 10-7 seconds, and by
  2. "convert them to a date in maple", you mean with respect to the "start time" of 00.00am on January 1, 1970 (UTC)

then the attached should fulfil the requirement.

Depending on your local time zone and default date format, the resulting date output may be different.

  restart;
  JStickToDate:= proc( nT::integer, startDate:=Date( 1970, 1, 1, 0, 0, timezone="UTC" ) )
                       uses Calendar:
                       local nD:= trunc( nT/24/60/60/10^7 ),
                             nH:= trunc( (nT-nD*24*60*60*10^7 )/60/60/10^7 ),
                             nM:= trunc( (nT-(nD*24+nH)*60*60*10^7)/60/10^7 ),
                             nS:= trunc( (nT-((nD*24+nH)*60+nM)*60*10^7)/10^7 ):
                       return AdjustDateField
                              ( AdjustDateField
                                ( AdjustDateField
                                  ( AdjustDateField
                                    ( startDate,
                                      "day_of_month",
                                      nD
                                    ),
                                    "hour",
                                    nH
                                  ),
                                  "minute",
                                  nM
                                ),
                                "second",
                                nS
                              );
                 end proc:
#
# Usage: first argument is number of ticks;
#        (optional) second argument is start date
#
  JStickToDate(12345678901234567);
  JStickToDate(13345678901234567);
  JStickToDate(14345678901234567);

_m615572320

 

_m725718304

 

_m725342496

(1)

 

Download TickToDate.mw

 

 

the attached is a bit simpler and neater

  restart:
  with(plots):
  with(plottools):
  with(NumberTheory):
  randomize():
  rnge:=100:
  n:=rand(1..rnge)():
  m:=rand(n+1..rnge+n+1)():
#
# Check whether "triple" is primitive or not
# and set plot options appropriately
#
  if  `and`
       ( AreCoprime(m,n),
         `or`
          ( type(m, even),
            type(n, even)
          )
       )
  then cl:= red:
       tl:= "Primitive Pythagorean Triple":
  else cl:= blue:
       tl:= "Pythagorean Triple":
  fi:
#
# Get sidelengths
#
  a:=m^2-n^2:
  b:=2*m*n:
  c:=m^2+n^2:

  display( [ polygon( [[0, 0], [0, a], [b, 0]],
                      color=cl
                    ),
             textplot( [ [ 0,   a/2, "a", 'align'={'left'}],
                         [ b/2, 0,   "b", 'align'={'below'}],
                         [ b/2, a/2, "c", 'align'={'above', 'right'}]
                      ],
                      font= [times, bold, 16]
                    )
           ],
           title=tl,
           titlefont=[times, bold, 20],
           axes=none
         );
  

 

 


 

Download pyTrip.mw

cat() does not evaluate what it concatenates. From the help page for cat() - emphasis added

The cat function is commonly used to concatenate strings and names together. The arguments are concatenated to form either a string, a name, or an object of type `||` .

So you need an "extra" eval() command, as in the attached

Example: Assume expressions assigned to names

xa, xb, xc := y = x, y = x^2, x = x^3

y = x, y = x^2, x = x^3

(1)

where you want to apply a command (rhs in this example) to all names fitting to a certain scheme (here a and c)

seq(cat(x, i), i = [a, c])

xa, xc

(2)

`~`[rhs]([xa, xc])[]

x, x^3

(3)

Combining (2) and (3) in one line

seq(rhs(eval(cat(x, i))), i = [a, c])

x, x^3

(4)

 

NULL

Download catEval.mw

You can use

simplify(expr, zero)

as shiown in the final execution group of the attached - but I don't see anywhere in yur code whihc produces

I want to  return {x = 0., y = 1.158748796 + 0. I}  as  {x = 0., y = 1.158748796 }. 

When I check the output of

soln3:= fsolve({b1, b2}, {x = 0 .. infinity, y = 0 .. infinity});
soln4:= fsolve({b1, b2}, {x = -infinity ..0, y = -infinity .. 0}); 

both are producing {x = 0., y = 0.} ??????

restart:

doCalc:= proc( xi )

                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.0001,
                       eta__1:= 0.240e-1,
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
                       theta_init:= theta0*sin(Pi*z/d),
                       n:= 30,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar,
                       soln3, soln4, Imagroot1, Imagroot2;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

  g:= q-(1-alpha)*tan(q)+ c*tan(q):
  f:= subs(q = x+I*y, g):
  b1:= evalc(Re(f)) = 0:
  b2:= y-(1-alpha)*tanh(y) -(alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 0:
  qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
  OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

  ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
  qstarTemporary:= min(ModifiedOddAsym);
  indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
  qstar2:= OddAsymptotes(indexOfqstar2);

# Compute complex roots

  AreThereComplexRoots:= type(true, 'truefalse');
  try
   soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
   soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
   qcomplex1:= subs(soln1, x+I*y);
   qcomplex2:= subs(soln2, x+I*y);
   catch:
   AreThereComplexRoots:= type(FAIL, 'truefalse');
  end try;

# Compute the rest of the Roots
  soln3:= fsolve({b1, b2}, {x = 0 .. infinity, y = 0 .. 10});
  printf("%a\n", soln3);
  soln4:= fsolve({b1, b2}, {x = -infinity ..0, y = -infinity .. 0});
  printf("%a\n", soln4);
  Imagroot1:=subs(soln3, I*y);
  Imagroot2:= subs(soln4, I*y);
  OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));

  if AreThereComplexRoots
  then gg:= [qcomplex1, qcomplex2,op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  elif not AreThereComplexRoots
  then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
  end if:

# Remove the repeated roots if any & Redefine n

  qq:= MakeUnique(gg):
  m:= numelems(qq):

## Return all the plots
            return qq, Imagroot1,Imagroot2, p, soln3, soln4;
  end proc:

ans:=[doCalc(0.06)];

{x = 0., y = 0.}
{x = 0., y = 0.}

 

[[4.675134212, 7.845062056, 10.99219719, 14.13554284, 17.27785480, 20.41979538, 23.56157656, 26.70328023, 29.84494257, 32.98658122, 36.12820550, 39.26982064, 42.41142975, 45.55303474, 48.69463685, 51.83623687, 54.97783536, 58.11943271, 61.26102919, 64.40262499, 67.54422027, 70.68581513, 73.82740965, 76.96900390, 80.11059793, 83.25219178, 86.39378547, 89.53537904, 92.67697249, 95.81856586], 0.*I, 0.*I, p, {x = 0., y = 0.}, {x = 0., y = 0.}]

(1)

a:={ 2+0.*I, 3 +0.*I};
simplify~(a, zero);

{2.+0.*I, 3.+0.*I}

 

{2., 3.}

(2)

 

Download oddProc3.mw

now uo are creating a table called 'w' and there are "special" considerations when return (potentially large) data structures - I' probably be using

return eval(w);

as in the attached

restart:

doCalc:= proc( xi )

                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.0001,
                       eta__1:= 0.240e-1,
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
                       theta_init:= theta0*sin(Pi*z/d),
                       n:= 30,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, w, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar,
                       soln3, soln4, Imagroot1, Imagroot2;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

                       qq := [2.106333379+.6286420119*I, 2.106333379-.6286420119*I, 4.654463885, 7.843624703, 10.99193295,14.13546782, 17.27782732, 20.41978346, 23.56157073, 26.70327712, 29.84494078, 32.98658013,         36.12820481, 39.26982019, 42.41142944, 45.55303453, 48.69463669, 51.83623675, 54.97783528,                     58.11943264, 61.26102914, 64.40262495, 67.54422024, 70.68581510, 73.82740963,                                 76.96900389, 80.11059792, 83.25219177, 86.39378546, 89.53537903, 92.67697249];
                        m:= numelems(qq):

                        for i from 1 to m do
                        w[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1);
                        end do;

## Return all the plots
  return eval(w);
  end proc:

ans:=doCalc(0.06);
ans[1];
ans[3];

table( [( 1 ) = 19.26175258-5.902958223*I, ( 2 ) = 19.26175258+5.902958223*I, ( 3 ) = 6.932751694, ( 4 ) = 2.753709683, ( 5 ) = 1.451748259, ( 6 ) = .8907623630, ( 7 ) = .6006863100, ( 9 ) = .3252930815, ( 8 ) = .4319166022, ( 11 ) = .2033707803, ( 10 ) = .2537145109, ( 13 ) = .1390101296, ( 12 ) = .1666327589, ( 15 ) = .1009699238, ( 14 ) = .1177214819, ( 18 ) = 0.6764836348e-1, ( 19 ) = 0.6014940869e-1, ( 16 ) = 0.8755294957e-1, ( 17 ) = 0.7664134466e-1, ( 22 ) = 0.4385090431e-1, ( 23 ) = 0.3987062318e-1, ( 20 ) = 0.5383110866e-1, ( 21 ) = 0.4845811396e-1, ( 27 ) = 0.2835149444e-1, ( 26 ) = 0.3071130419e-1, ( 25 ) = 0.3337839456e-1, ( 24 ) = 0.3640850864e-1, ( 31 ) = 0.2118804423e-1, ( 30 ) = 0.2270014307e-1, ( 29 ) = 0.2438004471e-1, ( 28 ) = 0.2625352267e-1 ] )

 

19.26175258-5.902958223*I

 

6.932751694

(1)

 

Download oddProc2.mw

 

1 2 3 4 5 6 7 Last Page 3 of 191