dharr

Dr. David Harrington

3120 Reputation

17 Badges

17 years, 347 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 replies submitted by dharr

@Carl Love Thanks. Edge {3,6} is missing on the 2 ohm one - still planar after adding it, 2 ohms between vertices 1 and 15.

@Tokoro No-one solved your first one-ohm problem. As I said before, it seems hard to search for these solutions, so if you have a method, you should share it.

@Tokoro The is a nice answer, and quite unexpected. Across which vertices is the equivalent resistance equal to one? Please post your code creating the graph so we can play with it. All vertex degrees are three except two...

I think that the rules should have: "all currents are different and none of them are zero", otherwise my second answer (and adding a pendent edge to your answer) is stlll valid.

I didn't want to do much work here since it wasn't really clear to me that it wasn't a trick question. But after seeing it I wondered how easy it would be to search for a solution by generating connected graphs with iteratively more vertices, with all combinations of edges to the previous vertices (this generates duplicates and might not be the most efficient). However, after 8 vertices it seemed like it would never get to 9. So the real question is - did you find this with Maple?

@evanhowington Hmm - very interesting, I thought the strict inequality case might be clearer. At least in Maple 15:

is(1+I < 2+I);
FAIL

This is what I would expect:  "<" not defined for complex. Now

is(1+I < 1+I);
false

I would still expect FAIL here - if you can't define it, then you can't answer the question. But it seems that Maple decides that if they are equal, then a decision can be made. Then the same for the symbolic case:

is(x<x) assuming x::complex
false

 

@Earl You're welcome. I think plotting is always useful. But also, as @tomleslie pointed out, letting Maple do the definite integral gives immediately the correct answer - Maple is smart about understanding discontinuities and dealing with them. You can also use value() if you have an unevaluated integral (gray integral sign). For example:. 

int1 := ArcLength(Ellipse, theta = 0 .. 2*Pi, output = integral);
value(int1);

 

@gtbastos Here are some thoughts.

There are some rejection tests you can do, for example if N is a power of M then they will commute, so if they don't commute there is no solution. Probably inspecting the eigenvalues is the best way, but if they aren't in the field then there is extra work as @Carl Love pointed out - he has a nice post about that here: https://www.mapleprimes.com/questions/203977-How-To-Find-Roots-Of-Polynomial-In-Finite#answer215097

You can do it more easily if both M and N happen to have only one root in the field -  then you can inspect all 5^9 powers of the eigenvalue of M to see which matches with the eigenvalue of N. If you can't factor at all, you can do the same thing with determinants - see attached - this takes only 15 seconds on my laptop. For determinants you may have to test more or revert to the full eigenvalue method. Is the maximum order of a matrix in the product group 5^9-1?

Both these methods can give false positives, so you then could test if the power you found actually works. A power of n or higher, where n is the matrix size, can be expressed as a linear combination of n-1 and lower powers, so that can be an efficient test (circumventing the discrete log problem to some extent). To test M^i, note the polynomial p_i(x) =q(x)p_M(x)+r(x), where p_i(x)=x^i, p_M(x) is the characteristic polynomial of M and q(x) and r(x) are the quotient and remainder when dividing x^i by p_M(x). If you take x=M in this polynomial, then p_M(M)=0 by Cayley-Hamilton, and M^i = r(M), so evaluate r(M) and see if it is N.

restart;

with(LinearAlgebra:-Generic):

Shorthand for matrix multiplication in field G

`&*`:= (A::Matrix,B::Matrix)->LinearAlgebra:-Generic:-MatrixMatrixMultiply[G](A,B):

Set up field. Extension ex is polynomial of degree n that is irreducible mod p. GF will randomly choose one if not supplied, but supplying one makes t global and a little easier to work with.

p:=5;n:=9;q:=p^n;
ex:=Nextprime(t^n,t) mod p;
G:=GF(p,n,ex);

5

9

1953125

t^9+t^2+2*t+3

_m143781824

Matrix may be entered in integer form and then converted to internal form

iM:=Matrix(3,3,[[98548,360142,206555],[841913,414734,1740986],[155301,596971,414440]]);
M:=G:-input~(iM);

iM := Matrix(3, 3, {(1, 1) = 98548, (1, 2) = 360142, (1, 3) = 206555, (2, 1) = 841913, (2, 2) = 414734, (2, 3) = 1740986, (3, 1) = 155301, (3, 2) = 596971, (3, 3) = 414440})

M := Matrix(3, 3, {(1, 1) = modp1(ConvertIn(t^7+t^6+t^5+2*t^4+3*t^3+t^2+4*t+3, t), 5), (1, 2) = modp1(ConvertIn(4*t^7+3*t^6+t^4+t^3+3*t+2, t), 5), (1, 3) = modp1(ConvertIn(2*t^7+3*t^6+t^5+2*t^3+2*t^2+t, t), 5), (2, 1) = modp1(ConvertIn(2*t^8+3*t^6+4*t^5+2*t^4+t^2+2*t+3, t), 5), (2, 2) = modp1(ConvertIn(t^8+t^6+2*t^5+3*t^4+2*t^3+4*t^2+t+4, t), 5), (2, 3) = modp1(ConvertIn(4*t^8+2*t^7+t^6+2*t^5+2*t^3+4*t^2+2*t+1, t), 5), (3, 1) = modp1(ConvertIn(t^7+4*t^6+4*t^5+3*t^4+2*t^3+2*t^2+1, t), 5), (3, 2) = modp1(ConvertIn(t^8+2*t^7+3*t^6+t^5+3*t^2+4*t+1, t), 5), (3, 3) = modp1(ConvertIn(t^8+t^6+2*t^5+3*t^4+2*t^2+3*t, t), 5)})

Try N = M^9

N:=M &* M: #M^2
N:=N &* N: #M^4
N:=N &* N: #M^8
N:=N &* M; #M^9

N := Matrix(3, 3, {(1, 1) = modp1(ConvertIn(t^8+t^7+2*t^6+t^5+t^2+t+1, t), 5), (1, 2) = modp1(ConvertIn(3*t^8+4*t^7+2*t^4+t^3+4*t^2+2*t+1, t), 5), (1, 3) = modp1(ConvertIn(3*t^7+2*t^6+2*t^5+2*t^4+t^3+t+2, t), 5), (2, 1) = modp1(ConvertIn(2*t^8+4*t^7+3*t^6+2*t^5+2*t^3+2*t+2, t), 5), (2, 2) = modp1(ConvertIn(3*t^8+t^7+4*t^5+4*t^3+t^2+4*t+3, t), 5), (2, 3) = modp1(ConvertIn(t^8+2*t^7+4*t^6+2*t^4+4*t^3+4*t+3, t), 5), (3, 1) = modp1(ConvertIn(t^8+4*t^7+4*t^6+3*t^5+3*t^2+1, t), 5), (3, 2) = modp1(ConvertIn(4*t^7+4*t^6+t^5+t^2+4*t, t), 5), (3, 3) = modp1(ConvertIn(3*t^8+3*t^7+4*t^6+3*t^5+3*t^4+3*t^3+3*t^2+4*t+4, t), 5)})

If they are powers then they should commute

EqualEntries(N &* M, M &* N);

true

Generate the determinants

dM:=Determinant[G](M);
dN:=Determinant[G](N);

modp1(ConvertIn(3*t^8+t^6+3*t^4+t^3+t+1, t), 5)

modp1(ConvertIn(t^8+2*t^6+4*t^5+4*t^4+t^3+2*t^2+2*t, t), 5)

Check all q-1 powers of the determinant of M to see if they match the determinant of N

st:=time():
prod:=dM:
for i to q-1 do;
  if G:-`=`(prod,dN) then print(i); end if;
  prod:=G:-`*`(prod,dM);  
end do:
time()-st;

9

14.094

 

 

Download GF.mw

@Carl Love OK, I see your point about the hard solve. Thinking about the eigenvalues sheds some light on how hard it will be to find a solution, even if it is not computationally efficient to evaluate them. Of course there are many ways to (partly) locate and decide the types of eigenvalues, which may mean you can decide there is no solution before you get to the hard solve. A trivial example would be that there will be no solution if M is singular but N is not. Not sure how Gershgorin disks and the like carry over for finite fields.

One ends up with more equations than unknowns. Perhaps it is easier to decide if the equation are incompatible and there is no solution than is is to try to solve them.

 

@Carl Love I'm not really serious about those, but it would be nice to get some clearer rules before I spend more time on this. If they are not OK I have a vague idea that 1 = 1/2 + 1/4 + 1/8 ... or similar series might form the basis of an infinite network solution.

@Joe Riel  But then the other three resistors have the same (zero) current. I suppose having one of the four unconnected with zero current is OK. But perhaps the problem doesn't insist on four resistors?

In version 2015 it just returns unevaluated. But if you assume all variables positive it returns something similar to Mathematica's answer and presumably equivalent to it.

int.mw

@escorpsy Could you please be more specific about what you want. CauchyPrincipalValue is for singularites in the integration range, https://en.wikipedia.org/wiki/Cauchy_principal_value. Perhaps you want to do contour integration, in which case singular and residue can de used. Or perhaps your question is unrelated to the integral you gave, and is about the principal branch of the log function. The help page for ln says:

For complex-valued expressions x,  ln(x) = ln(|x|) + I*argument(x), where  -Pi<argument(x)<=Pi. Throughout Maple, this computation is taken to be the definition of the principal branch of the logarithm.
 

@escorpsy I don't understand. evalc(ln(x+I*n)) (i.e. assuming x and n real) returns (1/2)*ln(n^2+x^2)+I*arctan(n, x), which explains where all the arctans are coming from. I don't see simplify getting rid of imaginary parts. In terms of the integral, if you try integrating from 0..p you will be asked about various assumptions, so I'm not sure Maple can correctly do this integral without knowing more about the symbolic parameters.

As @Axel Vogt points out there are problems with the integral to infinity for nonzero p. You can try singular(integrand,q) to find the singularities. 

@escorpsy 

II:=int(q*(ln((-(w+t)^2+(q+p)^2-I*n)/(-(w+t)^2+(q-p)^2-I*n))
    +ln((-(-w+t)^2+(q+p)^2-I*n)/(-(-w+t)^2+(q-p)^2-I*n))), q=0..infinity,'CauchyPrincipalValue'):
simplify(II) assuming positive;

Leads to infinity plus a long expression.

@lcz Very nice. Looking through the paper, it seems that rather than take the pseudoinverse of L, you can take the regular inverse of L+(1/n)*J, where J is the matrix of 1's. (The wikipedia page seems to have missed this point.) Since you are not making use of GraphTheory's LaplacianMatrix, and making your own weighted version, you may as well generate the above matrix directly (called LJ in the code). Its diagonals are chosen so the row/column sums are 1. I'm assuming the regular inverse is more efficiently calculated than the Moore-Penrose pseudoinverse, though I didn't check. This leads to the code below (works for weighted and unweighted graphs).

ResistanceMatrix:=proc(G::Graph)
  local W,M1,LJ,n,X;
  uses GraphTheory;
  if IsDirected(G) then error "Must be an undirected graph" end if;
  if IsWeighted(G) then W:=WeightMatrix(G) else W:=WeightMatrix(MakeWeighted(G)) end if;
  n:=upperbound(W)[1];
  M1:=Matrix(n,n, (i,j)-> if W[i,j]<>0 then 1/n-1/W[i,j] else 1/n end if,'shape'='symmetric');
  LJ:=M1+LinearAlgebra:-DiagonalMatrix(Vector(n,1)-MTM:-sum(M1,2));
  X:=LinearAlgebra:-MatrixInverse(LJ);
  Matrix(n,n,(i,j)->X[i,i]+X[j,j],'shape'='symmetric')-2*X;
end proc:

weightedResistanceMatrix.mw

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