dharr

Dr. David Harrington

1018 Reputation

13 Badges

15 years, 220 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

A.B-1 is rank 1 so there is only one nonzero eigenvalue. If this eigenvalue is found symbolically, then parameter values may be varied in plots as shown.
 

restart

with(LinearAlgebra); A := Matrix(6, 6, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = alpha[1]*S[p], (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = beta, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0}); B := Matrix(6, 6, {(1, 1) = mu, (1, 2) = tau, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = e[o], (2, 1) = 0, (2, 2) = tau+mu, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = beta+eta[1]+sigma+mu, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (1-delta)*alpha[2]-mu, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = alpha[2]*delta, (5, 5) = mu+eta[2], (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = mu+e[o]}); C := A.(1/B); Rank(C); evs := Eigenvalues(C); eig := op(`minus`({entries(evs, nolist)}, {0}))

A := Matrix(6, 6, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = alpha[1]*S[p], (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = beta, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0})

B := Matrix(6, 6, {(1, 1) = mu, (1, 2) = tau, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = e[o], (2, 1) = 0, (2, 2) = tau+mu, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = beta+eta[1]+sigma+mu, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (1-delta)*alpha[2]-mu, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = alpha[2]*delta, (5, 5) = mu+eta[2], (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = mu+e[o]})

C := Matrix(6, 6, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = alpha[1]*S[p]/(beta+eta[1]+sigma+mu), (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = beta/(beta+eta[1]+sigma+mu), (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0})

1

evs := Vector(6, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = alpha[1]*S[p]/(beta+eta[1]+sigma+mu)})

alpha[1]*S[p]/(beta+eta[1]+sigma+mu)

params := {beta = .5, delta = .2115, mu = 0.2041e-1, sigma = .9, tau = .33, S[p] = 34.722406639004, alpha[1] = 0.2e-3, alpha[2] = .2, e[o] = .33, eta[1] = 0.96e-1, eta[2] = 0.2485e-2}

{beta = .5, delta = .2115, mu = 0.2041e-1, sigma = .9, tau = .33, S[p] = 34.722406639004, alpha[1] = 0.2e-3, alpha[2] = .2, e[o] = .33, eta[1] = 0.96e-1, eta[2] = 0.2485e-2}

Plot vs mu, for example

plot(eval(eig, remove(has, params, mu)), mu = 0 .. 1)

 


 

Download eval.mw

For regular 1-D input, a rectangle surrounds the parenthesis that matches the one next to the cursor (vertical line here)

If you don't try to see R0 (terminate with colon), then implicitplot is fine with gridrefine=3 (gridrefine=10 runs into a memory allocation on my machine), and the data is extracted with plottools:-getdata without issue.

implictplot.mw

The ISOSURFACE plot structure (see ?plot,structure) is used internally by implicitplot3d to draw a surface where a function is zero over an evenly spaced grid of values, so it is close to your dataarray, assuming your x,y,z values are evenly spaced. So to get a red surface enclosing -1 values can be done as below. Not sure it makes sense to color the rest of space green since then you couldn't see anything.

restart;with(plots):

Generate dataarray with -1 values within an ellipse, +1 outside

ellipse:=(x,y,z)->if x^2/0.2+y^2/0.8+z^2/0.3 < 1 then -1 else 1 end if;

proc (x, y, z) options operator, arrow; if x^2/.2+y^2/.8+z^2/.3 < 1 then -1 else 1 end if end proc

x:=Vector(1..10,i->(i-5.5)/5):
y:=Vector(1..10,i->(i-5.5)/5):
z:=Vector(1..10,i->(i-5.5)/5):
dataarray:=Array(1..10,1..10,1..10,(i,j,k)->ellipse(x[i],y[j],z[k]),datatype=float[8],order=C_order):

Fill array required for ISOSURFACE structure

isoarray:=Array(1..10,1..10,1..10,1..4,datatype=float[8],order=C_order):

for ix to 10 do
  for iy to 10 do
    for iz to 10 do
      isoarray[ix,iy,iz,1]:=x[ix];
      isoarray[ix,iy,iz,2]:=y[iy];
      isoarray[ix,iy,iz,3]:=z[iz];
      isoarray[ix,iy,iz,4]:=dataarray[ix,iy,iz];
   end do;
 end do;
end do:

display(PLOT3D(ISOSURFACE(isoarray)),view=[-1..1,-1..1,-1..1],color=red);

 


 

Download isosurface2.mw

 

Try arctan(-diff(fy,t),diff(fx,t)), or perhaps arctan(diff(fy,t),diff(fx,t)) I didn't follow the details, but it removes the discontinuity.

It's just very slow, especially for small and large times. You are doing an integral within an integral. If you change the plot range to 1..1000 then you do get a plot in a few seconds.

In addition to @tomleslie 's solutions there is one real one between 0 and 1. I have also made the mininum number of changes to convert from the old linalg package to the LinearAlgebra package. To more consistently find the roots, you could use the NextZero routine in the RootFinding package.

I can't see any typos, so I don't know where the 0.5 root could come from.


 

restart

with(LinearAlgebra)

V0 := 1

1

a := 2

2

b := 1

1

assume(En > 0)

ks := sqrt(2*En)

2^(1/2)*En^(1/2)

qs := sqrt(2*(V0-En))

(2-2*En)^(1/2)

psi11 := proc (x) options operator, arrow; A1*exp(qs*x)+B1*exp(-qs*x) end proc

proc (x) options operator, arrow; A1*exp(qs*x)+B1*exp(-qs*x) end proc

psi21 := proc (x) options operator, arrow; A2*exp(I*ks*x)+B2*exp(-I*ks*x) end proc

proc (x) options operator, arrow; A2*exp(I*ks*x)+B2*exp(-I*ks*x) end proc

lp := a+b

3

psi12 := proc (x) options operator, arrow; exp(I*Kp*lp)*(A1*exp(qs*(x-lp))+B1*exp(-qs*(x-lp))) end proc

proc (x) options operator, arrow; exp(I*Kp*lp)*(A1*exp(qs*(x-lp))+B1*exp(-qs*(x-lp))) end proc

psi22 := proc (x) options operator, arrow; exp(I*Kp*lp)*(A2*exp(I*ks*(x-lp))+B2*exp(-I*ks*(x-lp))) end proc

proc (x) options operator, arrow; exp(I*Kp*lp)*(A2*exp(I*ks*(x-lp))+B2*exp(-I*ks*(x-lp))) end proc

eq1 := psi11(0) = psi21(0)

A1+B1 = A2+B2

eq2 := (D(psi11))(0)-(D(psi21))(0)

A1*(2-2*En)^(1/2)-B1*(2-2*En)^(1/2)-I*A2*2^(1/2)*En^(1/2)+I*B2*2^(1/2)*En^(1/2)

eq3 := expand(psi21(a)-psi12(a))

A2*(exp(I*2^(1/2)*En^(1/2)))^2+B2/(exp(I*2^(1/2)*En^(1/2)))^2-(exp(I*Kp))^3*A1/exp((2-2*En)^(1/2))-(exp(I*Kp))^3*B1*exp((2-2*En)^(1/2))

eq4 := expand((D(psi21))(a)-(D(psi12))(a))

I*A2*2^(1/2)*En^(1/2)*(exp(I*2^(1/2)*En^(1/2)))^2-I*B2*2^(1/2)*En^(1/2)/(exp(I*2^(1/2)*En^(1/2)))^2-(exp(I*Kp))^3*A1*(2-2*En)^(1/2)/exp((2-2*En)^(1/2))+(exp(I*Kp))^3*B1*(2-2*En)^(1/2)*exp((2-2*En)^(1/2))

row1 := Vector[row]([1, 1, -1, -1])

row1 := Vector[row](4, {(1) = 1, (2) = 1, (3) = -1, (4) = -1})

row2 := Vector[row]([coeff(eq2, A1), coeff(eq2, B1), coeff(eq2, A2), coeff(eq2, B2)])

row2 := Vector[row](4, {(1) = (2-2*En)^(1/2), (2) = -(2-2*En)^(1/2), (3) = -I*2^(1/2)*En^(1/2), (4) = I*2^(1/2)*En^(1/2)})

row3 := Vector[row]([coeff(eq3, A1), coeff(eq3, B1), coeff(eq3, A2), coeff(eq3, B2)])

row3 := Vector[row](4, {(1) = -(exp(I*Kp))^3/exp((2-2*En)^(1/2)), (2) = -(exp(I*Kp))^3*exp((2-2*En)^(1/2)), (3) = (exp(I*2^(1/2)*En^(1/2)))^2, (4) = 1/(exp(I*2^(1/2)*En^(1/2)))^2})

row4 := Vector[row]([coeff(eq4, A1), coeff(eq4, B1), coeff(eq4, A2), coeff(eq4, B2)])

row4 := Vector[row](4, {(1) = -(exp(I*Kp))^3*(2-2*En)^(1/2)/exp((2-2*En)^(1/2)), (2) = (exp(I*Kp))^3*(2-2*En)^(1/2)*exp((2-2*En)^(1/2)), (3) = I*2^(1/2)*En^(1/2)*(exp(I*2^(1/2)*En^(1/2)))^2, (4) = -I*2^(1/2)*En^(1/2)/(exp(I*2^(1/2)*En^(1/2)))^2})

C := Matrix(4, 4, 0)

C := Matrix(4, 4, {(1, 1) = 1, (1, 2) = 1, (1, 3) = -1, (1, 4) = -1, (2, 1) = (2-2*En)^(1/2), (2, 2) = -(2-2*En)^(1/2), (2, 3) = -I*2^(1/2)*En^(1/2), (2, 4) = I*2^(1/2)*En^(1/2), (3, 1) = -(exp(I*Kp))^3/exp((2-2*En)^(1/2)), (3, 2) = -(exp(I*Kp))^3*exp((2-2*En)^(1/2)), (3, 3) = (exp(I*2^(1/2)*En^(1/2)))^2, (3, 4) = 1/(exp(I*2^(1/2)*En^(1/2)))^2, (4, 1) = -(exp(I*Kp))^3*(2-2*En)^(1/2)/exp((2-2*En)^(1/2)), (4, 2) = (exp(I*Kp))^3*(2-2*En)^(1/2)*exp((2-2*En)^(1/2)), (4, 3) = I*2^(1/2)*En^(1/2)*(exp(I*2^(1/2)*En^(1/2)))^2, (4, 4) = -I*2^(1/2)*En^(1/2)/(exp(I*2^(1/2)*En^(1/2)))^2})

for i to 4 do C[1, i] := row1[i]; C[2, i] := row2[i]; C[3, i] := row3[i]; C[4, i] := row4[i] end do

C

Matrix([[1, 1, -1, -1], [(2-2*En)^(1/2), -(2-2*En)^(1/2), -I*2^(1/2)*En^(1/2), I*2^(1/2)*En^(1/2)], [-(exp(I*Kp))^3/exp((2-2*En)^(1/2)), -(exp(I*Kp))^3*exp((2-2*En)^(1/2)), (exp(I*2^(1/2)*En^(1/2)))^2, 1/(exp(I*2^(1/2)*En^(1/2)))^2], [-(exp(I*Kp))^3*(2-2*En)^(1/2)/exp((2-2*En)^(1/2)), (exp(I*Kp))^3*(2-2*En)^(1/2)*exp((2-2*En)^(1/2)), I*2^(1/2)*En^(1/2)*(exp(I*2^(1/2)*En^(1/2)))^2, -I*2^(1/2)*En^(1/2)/(exp(I*2^(1/2)*En^(1/2)))^2]])

chareq := Determinant(C)

2*((2*I)*2^(1/2)*En^(1/2)*(2-2*En)^(1/2)*(exp(I*Kp))^6*exp((2-2*En)^(1/2))*(exp(I*2^(1/2)*En^(1/2)))^2-I*2^(1/2)*En^(1/2)*(2-2*En)^(1/2)*(exp(I*Kp))^3*(exp((2-2*En)^(1/2)))^2*(exp(I*2^(1/2)*En^(1/2)))^4-I*2^(1/2)*En^(1/2)*(2-2*En)^(1/2)*(exp(I*Kp))^3*(exp(I*2^(1/2)*En^(1/2)))^4+2*En*(exp(I*Kp))^3*(exp((2-2*En)^(1/2)))^2*(exp(I*2^(1/2)*En^(1/2)))^4-I*2^(1/2)*En^(1/2)*(2-2*En)^(1/2)*(exp(I*Kp))^3*(exp((2-2*En)^(1/2)))^2-(exp(I*Kp))^3*(exp((2-2*En)^(1/2)))^2*(exp(I*2^(1/2)*En^(1/2)))^4-2*En*(exp(I*Kp))^3*(exp(I*2^(1/2)*En^(1/2)))^4-I*2^(1/2)*En^(1/2)*(2-2*En)^(1/2)*(exp(I*Kp))^3+(2*I)*2^(1/2)*En^(1/2)*(2-2*En)^(1/2)*exp((2-2*En)^(1/2))*(exp(I*2^(1/2)*En^(1/2)))^2+(exp(I*Kp))^3*(exp(I*2^(1/2)*En^(1/2)))^4-2*En*(exp(I*Kp))^3*(exp((2-2*En)^(1/2)))^2+(exp(I*Kp))^3*(exp((2-2*En)^(1/2)))^2+2*En*(exp(I*Kp))^3-(exp(I*Kp))^3)/(exp((2-2*En)^(1/2))*(exp(I*2^(1/2)*En^(1/2)))^2)

chareq0 := simplify(subs(Kp = 0, chareq))

(-(2*I)*(2-2*En)^(1/2)*2^(1/2)*En^(1/2)+4*En-2)*exp(-(2-2*En)^(1/2)-(2*I)*2^(1/2)*En^(1/2))+(-(2*I)*(2-2*En)^(1/2)*2^(1/2)*En^(1/2)-4*En+2)*exp(-(2-2*En)^(1/2)+(2*I)*2^(1/2)*En^(1/2))+(-(2*I)*(2-2*En)^(1/2)*2^(1/2)*En^(1/2)-4*En+2)*exp((2-2*En)^(1/2)-(2*I)*2^(1/2)*En^(1/2))+(-(2*I)*(2-2*En)^(1/2)*2^(1/2)*En^(1/2)+4*En-2)*exp((2-2*En)^(1/2)+(2*I)*2^(1/2)*En^(1/2))+(8*I)*(2-2*En)^(1/2)*2^(1/2)*En^(1/2)

Real parts are all zero for En<1

plot(Re(chareq0), En = 0 .. 3, color = red, thickness = 3)

Imaginary parts are all zero for En>1

plot(Im(chareq0), En = 0 .. 3, color = red, thickness = 3)

sol1 := fsolve(chareq0, En = .2 .. .4); sol2 := fsolve(chareq0, En = 2 .. 2.5)

.2669813621+0.*I

2.375133245-0.*I

 


 

Download qm_maple_-_periodic_potentials.mw
 

As @acer suggested, they have different compressions (and formats). Corel PhotoPaint reports the original 125 KB file is 1 bit B&W, 300 dpi, CCITT group 3 subformat (compression), and the 133 KB one written by Maple is 72 dpi, 8 bit grayscale and LZW compression. LZW is a lossless compression, but it is interesting that the output resolution is lower.

Edit: aside from the resolution issue, the above is as indicated on the ImageTools[Formats] help page.

 

There was nothing wrong with  the actual code, but it hadn't been entered after a prompt (in a worksheet) and was somehow two separate regions. If you enter 1-D code as you did, you need to format it yourself. You get some nicer formatting with a code edit region if you are going to enter a lot of code.

For this site, uploading a worksheet as you did is probably best.

proc.mw

Edit: Using a startup code edit region for initialization and procedure definitions can be a nice way to do it.

The help for CayleyTableGroup specifies: The single required argument is the Cayley table (or multiplication table): an n by n
 Matrix or Array, m, describing the group with n elements, where the [i, j] entry contains the positive integer k such that the product of the ith and j h group elements is the kth group element.

So this needs to be a matrix of integers.


 

restart;interface(rtablesize=12);

10

with(GroupTheory):

ct := <<e | p | q | r | s | t | u | v | w | x | y | z>, <p | q | e | y | u | w | z | r | x | t | v | s>, <q | e | p | v | z | x | s | y | t | w | r | u>, <r | z | t | s | e | y | v | x | p | u | q | w>, <s | w | y | e | r | q | x | u | z | v | t | p>, <t | r | z | x | w | u | e | q | y | p | s | v>, <u | x | v | p | y | e | t | z | s | r | w | q>, <v | u | x | z | q | r | y | w | e | s | p | t>, <w | y | s | t | x | z | p | e | v | q | u | r>, <x | v | u | w | t | s | q | p | r | e | z | y>, <y | s | w | u | p | v | r | t | q | z | e | x>, <z | t | r | q | v | p | w | s | u | y | x | e>>;

ct := Matrix(12, 12, {(1, 1) = e, (1, 2) = p, (1, 3) = q, (1, 4) = r, (1, 5) = s, (1, 6) = t, (1, 7) = u, (1, 8) = v, (1, 9) = w, (1, 10) = x, (1, 11) = y, (1, 12) = z, (2, 1) = p, (2, 2) = q, (2, 3) = e, (2, 4) = y, (2, 5) = u, (2, 6) = w, (2, 7) = z, (2, 8) = r, (2, 9) = x, (2, 10) = t, (2, 11) = v, (2, 12) = s, (3, 1) = q, (3, 2) = e, (3, 3) = p, (3, 4) = v, (3, 5) = z, (3, 6) = x, (3, 7) = s, (3, 8) = y, (3, 9) = t, (3, 10) = w, (3, 11) = r, (3, 12) = u, (4, 1) = r, (4, 2) = z, (4, 3) = t, (4, 4) = s, (4, 5) = e, (4, 6) = y, (4, 7) = v, (4, 8) = x, (4, 9) = p, (4, 10) = u, (4, 11) = q, (4, 12) = w, (5, 1) = s, (5, 2) = w, (5, 3) = y, (5, 4) = e, (5, 5) = r, (5, 6) = q, (5, 7) = x, (5, 8) = u, (5, 9) = z, (5, 10) = v, (5, 11) = t, (5, 12) = p, (6, 1) = t, (6, 2) = r, (6, 3) = z, (6, 4) = x, (6, 5) = w, (6, 6) = u, (6, 7) = e, (6, 8) = q, (6, 9) = y, (6, 10) = p, (6, 11) = s, (6, 12) = v, (7, 1) = u, (7, 2) = x, (7, 3) = v, (7, 4) = p, (7, 5) = y, (7, 6) = e, (7, 7) = t, (7, 8) = z, (7, 9) = s, (7, 10) = r, (7, 11) = w, (7, 12) = q, (8, 1) = v, (8, 2) = u, (8, 3) = x, (8, 4) = z, (8, 5) = q, (8, 6) = r, (8, 7) = y, (8, 8) = w, (8, 9) = e, (8, 10) = s, (8, 11) = p, (8, 12) = t, (9, 1) = w, (9, 2) = y, (9, 3) = s, (9, 4) = t, (9, 5) = x, (9, 6) = z, (9, 7) = p, (9, 8) = e, (9, 9) = v, (9, 10) = q, (9, 11) = u, (9, 12) = r, (10, 1) = x, (10, 2) = v, (10, 3) = u, (10, 4) = w, (10, 5) = t, (10, 6) = s, (10, 7) = q, (10, 8) = p, (10, 9) = r, (10, 10) = e, (10, 11) = z, (10, 12) = y, (11, 1) = y, (11, 2) = s, (11, 3) = w, (11, 4) = u, (11, 5) = p, (11, 6) = v, (11, 7) = r, (11, 8) = t, (11, 9) = q, (11, 10) = z, (11, 11) = e, (11, 12) = x, (12, 1) = z, (12, 2) = t, (12, 3) = r, (12, 4) = q, (12, 5) = v, (12, 6) = p, (12, 7) = w, (12, 8) = s, (12, 9) = u, (12, 10) = y, (12, 11) = x, (12, 12) = e})

elems:=convert(ct[1],list);
for i to nops(elems) do f(elems[i]):=i end do: #use remember table

[e, p, q, r, s, t, u, v, w, x, y, z]

M:=map(f,ct);

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

G := Group(M,labels=elems);

_m231969395104

 


 

Download Cayley.mw

I don't get an error message for the plot in Maple 2015, just an empty plot. But eval(S,w=1) gives -1.853089027*10^329399 so your numbers seem to be unreasonably large.


 

restart;

A:=(x^3-13*x^2+55*x-73)/(x-1)=0:
p:=numer(lhs(A));
realroot(p)[];

x^3-13*x^2+55*x-73

[169/64, 339/128]


 

Download realroot.mw

[Edited based on Carl's comment]. At least for real parameters (where the Matrix is Hermitian), there must be four eigenvectors.Not setting phi=Pi circumvents the OPs error message, and 4 eigenvectors are found, but then substituting phi=Pi leads to an error.


 

restart;

with(LinearAlgebra):

assume(k1,real,k2,real,m1,real,m2,real,phi,real);

Leave phi unspecified for now

A := <
          <0, 0, exp(I*k1) + m1, exp(I*k2) + m2>|
          <0, 0, exp(I*phi)*(exp(-I*k2) + m2), exp(-I*k1) + m1>|
          <exp(-I*k1) + m1, exp(-I*phi)*(m2 + exp(I*k2)), 0, 0>|
          <exp(-I*k2) + m2, exp(I*k1) + m1, 0, 0>
     >;

A := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = exp(-I*k1)+m1, (1, 4) = exp(-I*k2)+m2, (2, 1) = 0, (2, 2) = 0, (2, 3) = exp(-I*phi)*(exp(I*k2)+m2), (2, 4) = exp(I*k1)+m1, (3, 1) = exp(I*k1)+m1, (3, 2) = exp(I*phi)*(exp(-I*k2)+m2), (3, 3) = 0, (3, 4) = 0, (4, 1) = exp(I*k2)+m2, (4, 2) = exp(-I*k1)+m1, (4, 3) = 0, (4, 4) = 0})

Under the assumptions that the parameters are real, the Matrix is Hermitian and so must have 4 eigenvectors

simplify(A-A^%H);

Matrix([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])

If we specifiy phi=Pi, we get the OP's error

evals,evecs:=LinearAlgebra:-Eigenvectors(eval(A,phi=Pi)):

Error, (in LinearAlgebra:-Eigenvectors) multiplicity mismatch

If we leave phi unspecified, then we don't get the error message, and get 4 eigenvectors

evals,evecs:=LinearAlgebra:-Eigenvectors(A):
nops(evals),nops(evecs);  #nops always reterns 3, as Carl points out
numelems(evals),numelems(evecs);

3, 3

4, 16

Now specify phi=Pi. Works for the eigenvalues, but not for the eigenvectors. Perhaps some manipulation could improve this, but I'm not sure

evals2:=eval(evals,phi=Pi);
eval(evecs,phi=Pi);

evals2 := Vector(4, {(1) = (2+m1^2+m2^2+2*m1*cos(k1)+2*m2*cos(k2))^(1/2), (2) = -(2+m1^2+m2^2+2*m1*cos(k1)+2*m2*cos(k2))^(1/2), (3) = (2+m1^2+m2^2+2*m1*cos(k1)+2*m2*cos(k2))^(1/2), (4) = -(2+m1^2+m2^2+2*m1*cos(k1)+2*m2*cos(k2))^(1/2)})

Error, numeric exception: division by zero

 


 

Download Hermitian2.mw

When legend is a list there need to be multiple curves, but you had one curve with 3 points. One way to do this is:

display(pointplot([28,.6481496576],color=blue),pointplot([28,.648149657615473],color=orange),pointplot([28,.6512873548],color=red),style=point,labels = ["k", "y(k)"], legend = ["10-digit precision", "15-digit precision", "Floating-point iteration"] ,legendstyle = [font = ["HELVETICA", 9], location = right]);
 

 

Not entirely certain I know what you want, but this may be a start.

[Edit: I changed pi (just a symbol) to Pi (the circle constant thing)]

At the beginning I define i to be sqrt(-1) as you want and D to be a regular variable rather than the differential operator.

Maple gives the answer as a sum over roots of a quartic polynomial. To get a simpler formula, you will need to give values to some variables or parhaps tell Maple their signs (using assuming)

Download Linear_System.mw

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