Dr. David Harrington

5799 Reputation

21 Badges

19 years, 236 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

@acer But evalf(D(f)(1e-9)); and fdiff(f, [1], [1e-9]); gave very different values; it seems more is going on here.


1. Just change MakeFunction to unapply and it should work in older versions.

2. Correct.

3. Divide the expression by sigma__v*gamma, then it becomes D(f)(x)+f(x)/x, where x= Gamma. Plot this for x=0..10 and it is seen to be always positive, i.e., the second term is always larger than the first.

4. Consider (f(0)-f(infinity)*exp(D(f)(0)*x)*P(x) + f(infinity). For P(x)=1, this has the correct values at 0 and infinity and the correct slope at 0. It's not great. But write P(x) = 1 +b*x+c*x^2). If you evaluate this at two x values and set equal to the numerical value of f at those values, you can get two equations in b and c that fsolve can solve. This might be better. You can always add more terms to the polynomial and make it fit at more places. It might improve.

Edit: I think you want something simple like the above function, so you can see what is going on? But the derivatives are unlikely to be accurate, and oscillations can occur with interpolation. So you can do much better as @acer showed, but then why not just evaluate the function (or its derivative) numerically?

geom3d:-gtetrahedron can be used to make this solid.

@Mike Mc Dermott Here's one way that one might think about this, but you run into the same problem. If the rules of the game are that you are not allowed to know the network structure, then I thnk you cannot easily do this.


Take an equivalent resistance expression.

Z := (R[3]+R[4])*(R[1]+R[2])/(R[1]+R[2]+R[3]+R[4]); Rs := [indets(Z, name)[]]


[R[1], R[2], R[3], R[4]]

Define the conductances corresponding ro the resistances

Gs := subs(R = G, Rs); RtoG := `~`[`=`](Rs, `~`[`^`](Gs, -1)); GtoR := `~`[`=`](Gs, `~`[`^`](Rs, -1))

[G[1], G[2], G[3], G[4]]

[R[1] = 1/G[1], R[2] = 1/G[2], R[3] = 1/G[3], R[4] = 1/G[4]]

[G[1] = 1/R[1], G[2] = 1/R[2], G[3] = 1/R[3], G[4] = 1/R[4]]

Is this a parallel connection, i.e., a sum of terms involving conductances? First convert to an admittance in terms of conductances - find four parallel connections (not the two that we know is the right answer if we are allowed to know the network structure.

Y := expand(normal(eval(1/Z, RtoG), expanded)); type(Y, `+`)



Take the first term and go back to resistances - find 4 series connections, but two of them are not going to simplify further

op(1, Y); expand(normal(eval(1/%, GtoR), expanded))



What if we tried a series connection first, i.e., a sum of terms involving resistances?

Z := expand(Z); type(Z, `+`)



Convert first term to admittance and conductance - we have the same problem.

op(1, Z); expand(normal(eval(1/%, RtoG), expanded))




Download SeriesParallelReduction.mw

@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



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)


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)


Divide by first example

q := value(Z1p)


q2 := normal(Z2/q)


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)




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




`&||`(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. 


@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.


@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.




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).)


@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.


[Edited to correct mistakes as noted below]

4 5 6 7 8 9 10 Last Page 6 of 59