dharr

Dr. David Harrington

7044 Reputation

21 Badges

20 years, 142 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 retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are replies submitted by dharr

@WA573 

restart

with(LinearAlgebra)

L := Matrix([[(1+I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp(-(1/2)*k), I*rho*(exp((1/2)*k)-exp(-(1/2)*k))/lambda], [I*rho*(exp((-k)*(1/2))-exp((1/2)*k))/lambda, (1-I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp(-(1/2)*k)]]); tL := simplify(Trace(L))

Matrix(%id = 36893490171766514140)

2*exp(-(1/2)*k)

X := Matrix([[-I*lambda*(1/2)-(1/2)*w, -rho], [rho, I*lambda*(1/2)+(1/2)*w]]); Trace(X)

Matrix(%id = 36893490171766521972)

0

X has a pair of eigenvalues = +/-mu; their average is at the origin. To shift the eigenvalues of L so their average is at the origin, we add a matrix S

S := -(1/2)*tL*IdentityMatrix(2); Q := S+L; simplify(Trace(Q))

Matrix(2, 2, {(1, 1) = -exp(-(1/2)*k), (1, 2) = 0, (2, 1) = 0, (2, 2) = -exp(-(1/2)*k)})

Matrix(%id = 36893490171884120892)

0

Now a rotation/contraction in the complex plane by multiplying by the ratio of the eigenvalues b will make bQ and X similar

b := Eigenvalues(X)[1]/Eigenvalues(Q)[1]; bQ := b*Q

(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2)*(1+cos(k))*exp(-(1/2)*k)*lambda/((1+cos(k))*((exp(-(1/2)*k))^4*cos(k)*rho^2+(exp(-(1/2)*k))^4*w^2*cos(k)+(exp(-(1/2)*k))^4*rho^2-(exp(-(1/2)*k))^4*w^2-2*(exp(-(1/2)*k))^2*cos(k)*rho^2-2*(exp(-(1/2)*k))^2*rho^2+cos(k)*rho^2+rho^2))^(1/2)

Eigenvalues(bQ, output = list); Eigenvalues(X, output = list)

[(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2), -(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2)]

[(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2), -(1/2)*((2*I)*lambda*w-lambda^2-4*rho^2+w^2)^(1/2)]

Now bQ and X are similar. So what Matrix relates them?

evalsbQ, evecsbQ := Eigenvectors(bQ)

evalsX, evecsX := Eigenvectors(X)

So we have 1/evecsX.X.evecsX = 1/evecsbQ.bQ.evecsbQ or X = evecsX.(1/evexcbQ).bQ.evecsbQ.(1/evecsX)

B := evecsX.(1/evecsbQ); simplify(X-B__.bQ.(1/B))

Matrix(%id = 36893490171955216380)

So putting it all together we have X = b*B.Q.(1/B) and b*B.Q.(1/B) = b*B.(S+L).(1/B) and b*B.(S+L).(1/B) = -b.exp(-(1/2)*k).Identity+b*B.L.(1/B)

X__L := -b*exp(-(1/2)*k)*IdentityMatrix(2)+b*B.L.(1/B)

simplify(X-X__L)

Matrix(%id = 36893490171909539588)

NULL

NULL

NULL

NULL

 

NULL

Download XintoL4.mw

Earlier less streamlined version:

Download XintoL3.mw

@WA573 So L and X are U_n and V_n of Eq 2.1 of the paper? But L is not the identity + another matrix so this can't be right. Your use of "constant" without any context in the original question was misleading, and it is still not clear what the rules for writing "X in terms of L" are. Please specify exactly what you mean. You mention eigenvalues, and if you mean the eigenvalues of X and L it might be possible, by considering these, to write X = A + B.L.B^(-1) or some similar expression, where A and B are matrices. Please check that the red k/2 in my answer above is correct and is not -k/2, which would likely make any transformations simpler.

@dharr Here's some progress that may also illustrate the difficulty. Since P factors into a cubic and quartic as does Q, you can find one solution by setting the cubic and quartic coefficients equal and solving, as I did here. But I could have set one cubic equal to half the other, and then doubled the correspomding quartic to get another solution. And if the quartic factors into a cubic and a linear factor (which it always does, though not always easily), then there is a choice of cubics to match up, giving some more freedom. I think you need to reformulate the problem to fix two variables (or change variables to a form with two less) in order to give Maple a problem without this freedom.

solve_test.mw

Here's an example where I set a[3]=b[3]=1 and now Maple easily finds a unique (?) solution

solve_testa3b3equal1.mw

Notice that the denominator is just the coefficient of y^3 in P, so if you had also chosen P normalized to have the coefficient of y^3 = 1, then the solution would be simpler.

@Steven_Huang With a symbolic parameter that is much harder. One strategy would be to rewrite your equations in terms of fewer variables, since some seem to be redundant. In your earlier problem, were your expecting it to return the coefficients you set in A1 and B1? That didn't happen, two parameters were arbitrary (Maple chose a[0] and a[3]) - is that what you expected?
 

@Khair Muhammad Saraz It's not very clear what you want. Here I generate the error plots as separate plots because they will be on a different scale; I haven't optimized the look of the plots.

DQMSixth.mw

@Carl Love Thanks. Yes, I usually forget 'nolist' as well. Since we only ever want one class representative, the code in ListTools:-Classify can be simplified to

classreps := proc(f, L);
  local class, e;
  for e in L do class[f(e)] := e end do;
  [entries(class, 'nolist')];
end proc

with your proc

f := g -> [seq](op(4, GraphTheory:-CanonicalGraph(g)))

Probably also one would want to change it is to always give generic vertices 1..n, though that could be delegated to the user for the input list.

@Carl Love I misinterpreted that statement as meaning you could have two different canonical orders even if the graphs were isomorphic. But I see your interpretation is correct. But the description is also correct; op(4,G) are the same for those graphs; the same problem also exists with @lcz's original list. The error has to do with a misplaced op(1,..) and perhaps 'nolist' for entries (always a problem!). One solution is 

op~(1,{entries}(ListTools:-Classify(
    g-> [seq](op(4, GraphTheory:-CanonicalGraph(g))),graph_list),'nolist'));

 

@Carl Love Yes, I use it when I can find such a procedure. In this case, it needs some work; the example is from the CanonicalGraph help page.

Download classify.mw

@lcz To add to @Carl Love 's explanation, op(4,G) is the same information you would get from Neighbors(G), but with generic labelling, and an Array of sets rather than a list of lists. See ?graph,details.

@mmcdara Very nice. Vote up.
@KIRAN SAJJAN Note that the paper has -Gr/(R*Re)*H' and not +Gr/(R*Re)*H'

@mmcdara Thanks! I was still bothered by not taking advantage of duplicate rows, and have now uploaded  above a version that can manage the L=10 case in about 1 s on my machine.

Since there is more likely not a match and that can be detected at the first mismatch, you want to try to not to generate all the sequence unless you have to (as in the McCarthy rules for 'and'). I don't think the timing tests are that reliable here. I suspect there are greater efficiencies to be found in the larger algorithm, so this is just a reply, not an answer.

tests.mw

@janhardo The error message had "unexpected option: [-50 .. 50, -5 .. 5, 0 .. 5] = [-50 .. 50, -5 .. 5, 0 .. 5]", where the option should have been "view = [-50 .. 50, -5 .. 5, 0 .. 5]", which suggests that the lhs view had unexpectedly obtained a value. That in turn suggests that the view parameter for complex_surface should be named something else. Another fix would be to have the complexplot3d option as ':-view' = view but I think renaming to something else avoids confusion.

@janhardo I fixed some small things. The error was fixed by using a different name than view for the complex_surface parameter. You weren't returning the generated plot because of the following printf statement (added to debug?).

restart;

complex_surface := proc(f, z_range::range(complex), plotview::list) local plt;
    
    printf("Complex function f(z) defined.\n");
    printf("Range of z values: %a\n", z_range);
    
    # Plot de complexe functie met complexplot3d
    printf("Plotting complex function...\n");
    plt := plots:-complexplot3d(f, z_range, 'view' = plotview, 'grid' = [50, 50]);
    
    printf("Plot generated.\n");
    
    plt   # return result
end proc:

# Definieer de complexe functie f(z) = ln(z)
f := z -> ln(z);

# Definieer het bereik van z-waarden voor de plot
z_range := -50-5*I .. 50+5*I;

# Roep de procedure aan om de plot te genereren
complex_surface(f, z_range, [-50..50, -5..5, 0..5]);

proc (z) options operator, arrow; ln(z) end proc

-50-5*I .. 50+5*I

Complex function f(z) defined.
Range of z values: -50-5*I .. 50+5*I
Plotting complex function...
Plot generated.

NULL

NULL

NULL

Download complexplot.mw

You might consider complexplot3d, e.g., 

plots:-complexplot3d(Zeta(z), z=-50-5*I..50+5*I,view=[-50..50, -5..5, 0..5],
                    grid=[50,50]);

gives

First 10 11 12 13 14 15 16 Last Page 12 of 72