Dr. David Harrington

5973 Reputation

21 Badges

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

Simple and resonably accurate


f := u -> RootOf(8*_Z^4 + 12*u*_Z^3 + (5*u^2 - 4)*_Z^2 - 4*u*_Z - u^2):
f0 := 1/sqrt(2):
finf := 1/sqrt(5):
Df0 := -1/4:

fapprox:=(f0-finf)*exp(-x/8)/(1+(5/6)*x)+finf; # 5/6 could be 0.837 for D(fapprox)(0) = -0.25002




Maximum relative error for 0 to infinity <2%



Download MaPal93approx.mw

@MaPal93 I remain confused about what you actually want:

1. a simple function (a few parameters) that gives good agreement for the function at 0 and infinity and for the derivative at 0, but is not necessarily that accurate.

2. an accurate approximation over a limited range.

3. an accurate approximation over the full range.

My objective was (1), because you asked for a 1-line function and you seemed to want to interpret it easily. That was also what I did in that paper. But now you are comparing as though you want (2), which you did ask for at some point and @acer gave a nice answer to. For (3) you can just use the function itself.

Since the comparison you provided is for (2) then a fair comparison requires the same number of parameters, which is 8 here in each case (though maybe 9 for the rational function; I'm not sure). (That's more that I was originally thinking of for (1).) Now plot them out to 200 and you will see that mine is OK at infinity (by design), but has the oscillation. The second one is a polynomial and will become large at infinity. The rational function can go to a constant value at infinity if the degree of the numerator and denominator polynomials are the same, and @acer chose [4,4], presumably for this reason. As @acer pointed out, the rational functions are usually better than the polynomials anyway, and you see that here.

The numapprox routines are focussed on accuracy over the whole range and they do well for that. They should also do well for derivatives at zero since they are based on some series expansions around that point. I'm assuming they are based on evaluating the function at evenly spaced points across the chosen interval or some criterion that treats all the interval evenly. If you get to choose the interpolation points then you can do better than if they are evenly spaced. In Gaussian quadrature, you optimize this, so (from memory) an (n/2)-degree polynomial with optimally chosen points can do as well as an n-degree polynomial at evenly spaced points. This is not quadrature, so I don't think the Laguerre points are optimal, but the basic idea is that they should spread out with more near the origin (where the function is changing) and fewer at large values (where it isn't changing much). But that again is assuming you want to approximate out to infinity. And then why not (3)? 

@MaPal93 Nonlinear equations can be tricky. DirectSearch, an external package from the Maple Applications Centre can solve this. You need some interpolation points closer to the origin, and then spreading out with fewer later where the function is featureless. One possibility is to use roots of Laguerre polynomials, which are spread out in this way, and are used in Guass-Laguerre quadrature for functions in a semi-infinite range; see the end of the file. But that is just a guess.


@MaPal93 That's a nice fit (0.2%) over that range. Df0 looks good. If you want more than the probe accuracy, you can use numapprox:-infnorm. You have oscillations outside the cutoff. If you want it to work over a longer range, you will have to spread the interpolation points out. High-degree polynomials can oscillate - you can likely remove the oscillations with a lower degree polynomial and still get reasonable accuracy. But maybe it is fine as it is.



  1. Should I consider looking into [15,16,17] to try and find a similarly accurate and simple approximation for my function? Would that be challenging yet worthy? Would you help?
    Those refs are about fitting to experimental data and not relevant. I developed that approximation by playing around using series and asympt and some intuition. In fact a referee asked me to explain how I had "derived it", but I was unable to answer that question and so only the vague description is is the paper. For your case the series/asympt expressions are complicated, so I don't think they will help; another approach will be required. Notice that the approximation was not accurate (up to 4%); the errors only had to be small compared to the experimental errors. I'll take a look to see if anything is obvious.
  2. You also mention "relative errors (with respect to Eq. (2))" and "systematic error in the parameters can be estimated by individually varying the parameters to find the minimum in the residual sum of squares". I think it could be interesting to quantify the errors for my approximation as well.
    The relative error is just from a plot of (fapprox-fexact)/fexact. The other errors are related to the fit to experimental data.
  3. Talking about interpolation instead, you mention "Exact value and derivative at zero preclude any of the things that fit (as opposed to interpolate) arbitrary functions unless they are carefully designed not to disrupt the exact values" and "c1*x term will mess up the derivative at zero". Which replacement term would preserve the derivative at 0?​​
    The simplest would be to replace the Df0 exponential decay constant with a parameter a, differentiate the whole thing and set its value at zero equal to Df0 as another equation to be solved.

@Preben Alsholm Vote up. Nice analysis.

@MaPal93 That is more or less what I meant, but you interpolated the polynomial not the overall function fapprox. Howerer, I should have given it a bit more thought before suggesting it. I suggested the exp form with P(x) = 1 because I wanted something that has the correct behaviour at infinity. I don't like the cutoff idea for that reason, and like to hack things that work well at both ends (for a real example see Eq. (13) in doi:10.1016/j.elecom.2019.106566). I also wanted a good derivative at 0. Exact value and derivative at zero preclude any of the things that fit (as opposed to interpolate) arbitrary functions unless they are carefully designed not to disrupt the exact values. Other methods work hard to match values and derivatives at one point (like the numapprox routines with polynomials) but need a cutoff. Since a decaying exponential dominates a polynomial at infinity, I suggested multiplying the exp by a polynomial. I made the constant 1 to keep the value at 0 correct, but forgot that the c1*x term will mess up the derivative at zero.

So there are tradeoffs about what you want to work well, the simplicity of the function and the ease of implementation.

@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(" &Vert; "), 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.

First 6 7 8 9 10 11 12 Last Page 8 of 62