dharr

Dr. David Harrington

8345 Reputation

22 Badges

21 years, 6 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

@Mike Mc Dermott Once you go from the network/graph to the effective resistance, it's hard to go back, unless you have a well-defined step-by-step algorithm.

On the other hand, if you are allowed to know about the graph, then it is easier. Is the expression (R1 + R2)||(R3 + R4) at the beginning of the second example deliberately the resistance of the first network without RB, or could it potentially have something to do with some completely diferent network?

If the main part of the exercise is combining series and parallel resistors successively to simplify a network, that could be done something like the following hack, which is intended just to give the flavor of what might be done. But I don't konw how you would finish off the second example, because n3/d3 doesn't seem obviously related to a known network.

See mapleprimes

First example"

I1 0 1
R1 1 2
R2 2 0
R3 1 3
R4 3 0

restart

with(GraphTheory)

From Carl Love:

`print/&||` := proc (A, B) Typesetting:-mrow(Typesetting:-Typeset(A), Typesetting:-mo(" ‖ "), Typesetting:-Typeset(B)) end proc; `value/&||` := proc (a, b) options operator, arrow; normal(a*b/(a+b)) end proc

Construct graph from nodes and resistances (edges), Assume all edges have different resistance names

R := table(); VertexInOut := {0, 1}; edge[1] := {1, 2}; R[edge[1]] := R1; edge[2] := {0, 2}; R[edge[2]] := R2; edge[3] := {1, 3}; R[edge[3]] := R3; edge[4] := {0, 3}; R[edge[4]] := R4; R := [entries(R, 'pairs')]; edges := convert(edge, list); verts := convert(`union`(edges[], VertexInOut), list); G := Multigraph(verts, edges); DrawGraph(G, size = [200, 200], layout = tree)

[{0, 2} = R2, {0, 3} = R4, {1, 3} = R3, {1, 2} = R1]

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

Find equiv resistance (routine in startup code)

Z1 := Resistance(G, VertexInOut, R)

(R3+R4)*(R1+R2)/(R1+R2+R3+R4)

Simplify series connections (pairwise for now). Find vertices of degree 2 that aren't the in or out vertices

verts := Vertices(G); map2(Degree, G, verts); seriesverts := `minus`({verts[[ListTools:-SearchAll(2, %)]][]}, VertexInOut)

{2, 3}

Eliminate these vertices and edges and add new edges with resistance the sum of these.
(DeleteEdges not properly handled here for a multigraph.)

"edgesdel:={}:  for i,vert in seriesverts do    nv[i]:={Neighbors(G,vert)[]};    deledges:={seq({vert,j},j in nv[i])}:     eval(deledges,R);    R:=[R[],nv[i]=add(`%`)] ;    R:=remove(has,R,deledges);    DeleteEdge(G,deledges);    edgesdel:=edgesdel union deledges;  end do:  nv:=convert(nv,list):   R;"

[{0, 1} = R1+R2, {0, 1} = R3+R4]

Update the graph - the 2 indicates 2 parallel edges

newedges := remove(has, edges, [edgesdel[]]); newedges := [newedges[], nv[]]; G := DeleteVertex(G, seriesverts); G := Multigraph(Vertices(G), newedges); DrawGraph(G, size = [200, 200])

Now parallel cases.
Just do this one edge manually

e := {0, 1}; if EdgeMultiplicity(G, e) > 1 then R, par := selectremove(proc (x) options operator, arrow; evalb(x = e) end proc, R); par := `&||`(map(rhs, par)[]); R := [R[], e = par] end if; R

[{0, 1} = `&||`(R1+R2, R3+R4)]

Z1p := rhs(op(R))

`&||`(R1+R2, R3+R4)

Second example
I1 0 1
R1 1 2
R2 2 0
R3 1 3
R4 3 0
RB 2 3

R := table(); VertexInOut := {0, 1}; edge[1] := {1, 2}; R[edge[1]] := R1; edge[2] := {0, 2}; R[edge[2]] := R2; edge[3] := {1, 3}; R[edge[3]] := R3; edge[4] := {0, 3}; R[edge[4]] := R4; edge[5] := {2, 3}; R[edge[5]] := RB; R := [entries(R, 'pairs')]; edges := convert(edge, list); verts := convert(`union`(edges[], VertexInOut), list); G := Multigraph(verts, edges); DrawGraph(G, size = [200, 200], layout = tree)

[{0, 2} = R2, {0, 3} = R4, {1, 3} = R3, {2, 3} = RB, {1, 2} = R1]

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

Find equiv resistance (routine in startup code)

Z2 := Resistance(G, VertexInOut, R)

(((R3+R4)*R2+(R4+RB)*R3+RB*R4)*R1+((R4+RB)*R3+RB*R4)*R2)/((R2+R4+RB)*R1+(R3+RB)*R2+(R4+RB)*R3+RB*R4)

Divide by first example

q := value(Z1p)

(R3+R4)*(R1+R2)/(R1+R2+R3+R4)

q2 := normal(Z2/q)

(R1*R2*R3+R1*R2*R4+R1*R3*R4+R1*R3*RB+R1*R4*RB+R2*R3*R4+R2*R3*RB+R2*R4*RB)*(R1+R2+R3+R4)/((R1*R2+R1*R4+R1*RB+R2*R3+R2*RB+R3*R4+R3*RB+R4*RB)*(R3+R4)*(R1+R2))

This limits to 1, but numerator and denominator each tend to infinity, so the coeffs of RB must be the same

limit(q2, RB = infinity); n2 := collect(numer(q2), RB); d2 := collect(denom(q2), RB)

1

(R1*R3+R1*R4+R2*R3+R2*R4)*(R1+R2+R3+R4)*RB+(R1*R2*R3+R1*R2*R4+R1*R3*R4+R2*R3*R4)*(R1+R2+R3+R4)

(R1+R2+R3+R4)*(R3+R4)*(R1+R2)*RB+(R1*R2+R1*R4+R2*R3+R3*R4)*(R3+R4)*(R1+R2)

cRB := coeff(n2, RB)*RB; n3 := map(simplify, collect(n2/cRB, RB)); d3 := map(simplify, collect(d2/cRB, RB))

1+(((R3+R4)*R2+R3*R4)*R1+R2*R3*R4)/((R3+R4)*(R1+R2)*RB)

1+(R2+R4)*(R1+R3)/(RB*(R1+R2+R3+R4))

Z1p*n3/d3

`&||`(R1+R2, R3+R4)*(1+(((R3+R4)*R2+R3*R4)*R1+R2*R3*R4)/((R3+R4)*(R1+R2)*RB))/(1+(R2+R4)*(R1+R3)/(RB*(R1+R2+R3+R4)))

``

Download SeriesParallelReduction.mw

@MaPal93 Good catch!

@MaPal93 You should be able to rely on limits directly from RootOf, and it is disappointing that some are not correct. Converting to radicals is only possible in some cases and if it is possible, the expression may be too complicated to get the limit, as you found. In general, Maple is not always reliable in limits of complicated expressions, so one can always do a numerical check. So it looks like various tests show the correct limits in the end - in one case the limit at infinity is a constant even if it looked initially like zero. 

caseAcaseB.mw

@MaPal93 I agree you can't tell much from those "instabilities", especially on infinity plots where successive x-axis values can be quite different. The plot routine uses hardware precision and the radical expression is complicated, so at extreme values higher accuracy (high Digits) may be needed for good numerics.

@MaPal93 Thanks for the appreciation.

The limit can be undefined if approaching from the left gives a different answer from appoaching from the right. Plotting for sigma__d from -20..20 siggests this is true; negative values give no plot (complex values). In principle,

limit(convert(lambda__2, radical), sigma__d=0, right)  assuming positive;  

should then return infinity, but I gave up after some time. Another way to check is to look at the first term of 

series(convert(lambda__2, radical), sigma__d, 3) assuming positive;  

which shows it does go to infinity - the other terms go to constant or zero values.

@jalal use M^+ to transpose a Matrix M.

S3TabProportionnAléatoire.mw

@MaPal93 Yes, I meant that plot. The values were found by taking limits, but RootOf is buggy here, so you have to convert to radicals first. For sensitivity analysis I assume you want the derivatives. Limits with respect to parameters can also be found.

MaPal93.mw

@jalal 

S3TabProportionnAléatoire.mw

The /~ operator does this. To enter it in 2D, you need: space \ / ~ space.

@dharr Since I think you only are interested in positive solutions, and now all the solutions returned by solve have been analysed, this is my summary of the situation, for the assumptions that you are interested in

For any positive values of the parameters, there exists a unique positive solution (i.e., with positive lambda__1, lambda__2 and lambda__3), in which

lambda__1 = sqrt(2/5)*sigma__v/sigma__d

lambda__2 = lambda__3 = f(Gamma)*sigma__v/sigma__d

where Gamma = gamma*sigma__v*sigma__d and f(Gamma) is a (known) complicated function of Gamma that lies between 1/sqrt(10) = 0.316 and 1/3 = 0.333.

I left out the trivial solution lambda__1 = lambda__2 = lambda__3 =0 becuse as you originally formulated the equations there is a division by zero error. After they are in normal form (single numerator / single denominator) it is a solution.

(The complicated solution is from converting the RootOf to radical form. If instead of solving the "redox"  equation for Lambda_2 you solve for Gamma, you get a simpler form for the inverse function Gamma = g(Lambda__2).)

QED?

@MaPal93 Thanks for noticing that. I meant "positive Lambda__2 and negative Lambda_3" in both cases, since I am only talking about the first solution that solve returns (first two roots) 

Showing that they are superimposed means

(1) that the statement about the maximum Gamma applies to both Lambda__2 and Lambda__3. I should have checked that algebraically but I was lazy.

(2) They could of course look wildly different, so the fact that when you put all the roots together they look the same means they is some interesting symmetry. I thought at first that it was just when you exchanged Lambda__2 and Lambda__3 the sign changed but it is not that simple. Some permutation of the roots turns one into the other.  (I tried fully factoring delta, but the result was an unenlightening mess.)

@MaPal93 That's a tough question, since a lot depends on what you want to do with the answers. It's always useful to try simplify, evalc, or convert(.., radical) to see if there is a closed form, which can be found for some higher degree polynomials.

Sometimes the complicated closed form solution is harder to work with than a RootOf.

And you can do things with RootOfs that give useful information. As an example consider the first solution with the 6th order polynomial. Plots show that for positive gamma there are solutions with positive Lambda__2 but negative Lambda__3, which may be enough to know that this is not useful to you. But suppose you wanted to know the maximum Gamma value for which it is possible to find a positive Lambda__2. You can differentiate and solve to find this value is Gamma_max = 0.29407, so  for any sigma__v and sigma__d and gamma, their product has to be less than 0.29407 to give a real solution (with a positive Lambda__2 and negative Lambda__3).

Some analysis of this type here.

NondimshortAllEqns.mw

[Edited to correct mistakes as noted below]

@MaPal93 Turns out you can only eliminate 3 of the 13 variables.

difficult_non-dimensionalization.mw

@dharr Here is the plot - the positive one is consistently the first (red) one.

solve(Eqq111, Lambda__2):
allvals:=map(simplify, [allvalues(%)]):
plot(allvals,Gamma=0..3,color=[red,blue,black,magenta]);

Nondimshort.mw

@acer In this case, the question is about non-dimensionalizing a more complicated set of equations. That is certainly not apparent from the title of the thread. I non-dimensionalized the simpler set by inspection as an aside to answering the question about the complexity, but this one will likely require the machinery of Hermite normal forms.

@MaPal93 

  • The only real and positive solution is for lambda_2 = lambda_3, as shown explicitly in Nondimshort.mw. Is this solution the same as ans23[2] and ans32[2] in NondimshortAllEqns.mw (i.e., Lambda_2 = Lambda_3 = beta)?
    Yes, if you look at the definition of beta in the alias, it is the same quartic as before.
  • In Nondimshort.mw you type "Seems to always have one positive root - the curve Jumps at three places (these places can be found with solve(...,parametric).)". How to do so?
    solve(Eqq111, Lambda__2, parametric);
    Moreover, this curve is Lambda_2 as a function of positive Gamma...how to obtain the original lambda_2?
    use the relationship lambda__2 = Lambda__2*sigma__v/sigma__d.
    If I type solve(Eq2redox, lambda__2); map(simplify, [allvalues(%)]); I obtain 4 roots...which one is the positive one?
    You cannot tell from looking at the solution, since it will depend on the values of the parameters. You have a better chance with solve(Eqq111, Lambda__2); map(simplify, [allvalues(%)]); because there is only one parameter Gamma, but even there it is very complicated to analyse. In principle, solve(Eqq111, Lambda__2, parametric,explicit); should give you the different solutions for different ranges of Gamma, but it only gives values at special values of Gamma, presumably because this case is too hard to solve. You can probably plot different solutions for the different ranges of Gamma between the critical values to figure this out.
    Is it unique?
    Descarte's rule of signs says there is one positive root, but which one it is will depend on the value of Gamma.
First 32 33 34 35 36 37 38 Last Page 34 of 87