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

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

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


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

solve(Eqq111, Lambda__2):
allvals:=map(simplify, [allvalues(%)]):


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


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

@Nicole Sharp Here is a worksheet with both your examples in which eval(q, Unit =1) works. Since you didn't provide the actual worksheet, it's hard to further diagnose. Or perhaps I have misunderstood your question?


@dharr Here's the full case that I think better answers your question about the apparent asymmetry that arises from the way the equations are solved, but really the answers are not that different.


local gamma;


# Consider a system of three equations:  

eq1 := sigma__v^2/(lambda__1*(sigma__v^2/(2*lambda__1^2) + (5*sigma__d^2)/4)):

eq2 := (gamma*sigma__v^2 + 2*lambda__3)*sigma__v^2/(((2*gamma*sigma__v^2 + 8*lambda__2)*lambda__3 + 2*gamma*lambda__2*sigma__v^2)*(2*(gamma*sigma__v^2 + 2*lambda__3)^2*sigma__v^2/((2*gamma*sigma__v^2 + 8*lambda__2)*lambda__3 + 2*gamma*lambda__2*sigma__v^2)^2 + lambda__2^2*(gamma*sigma__v^2 + 4*lambda__3)^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + gamma^2*sigma__v^4*lambda__3^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__2)*lambda__3 + 2*gamma*lambda__2*sigma__v^2)^2 + sigma__d^2)):

eq3 := (gamma*sigma__v^2 + 2*lambda__2)*sigma__v^2/(((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)*(2*(gamma*sigma__v^2 + 2*lambda__2)^2*sigma__v^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + lambda__3^2*(gamma*sigma__v^2 + 4*lambda__2)^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + gamma^2*sigma__v^4*lambda__2^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + sigma__d^2)):

Eq1, Eq2, Eq3 := eq1 - lambda__1, eq2 - lambda__2, eq3 - lambda__3:

"non-dimensionalize" Eq 1 by hand. New non-dim variable Lambda__1. End up with no parameters in the equation

Lambda__1 = lambda__1*sigma__d/sigma__v;

Lambda__1 = lambda__1*sigma__d/sigma__v

lambda__1 = Lambda__1*sigma__v/sigma__d


0, (1/5)*10^(1/2), -(1/5)*10^(1/2)

We can find the non-dimensionalization here by a systematic process, but let's guess Lambda__2 and Lambda__3 are  defined similarly to Lambda__1

{Lambda__2 = lambda__2*sigma__d/sigma__v,Lambda__3 = lambda__3*sigma__d/sigma__v};

{Lambda__2 = lambda__2*sigma__d/sigma__v, Lambda__3 = lambda__3*sigma__d/sigma__v}

{lambda__2 = Lambda__2*sigma__v/sigma__d, lambda__3 = Lambda__3*sigma__v/sigma__d}



We can divide through by sigma__v, which is in every term, and we see that Gamma = gamma*sigma__v*sigma__d

Gamma = gamma*sigma__d*sigma__v;

Gamma = gamma*sigma__d*sigma__v

gamma = Gamma/(sigma__d*sigma__v)



alias(alpha=RootOf((Gamma^2 + 4)*_Z^2 + 4*Gamma*_Z + Gamma^2),
beta=RootOf(40*_Z^4 + 36*Gamma*_Z^3 + (9*Gamma^2 - 4)*_Z^2 - 4*Gamma*_Z - Gamma^2),
delta=RootOf((200*Gamma^2 + 800)*_Z^6 + 560*Gamma*_Z^5 + (-24*Gamma^2 - 80)*_Z^4 + (-52*Gamma^3 - 64*Gamma)*_Z^3 + 16*Gamma^2*_Z^2 + 24*Gamma^3*_Z + 5*Gamma^4));

alpha, beta, delta

Solve in both orders - we see that we have the Lambda__2=Lambda__3 solution that we saw before. Lambda__2 and Lambda__3 are swapped, i.e. the apparent complexity depends on the solving order


{Lambda__2 = (1/5)*(-900*Gamma^7*delta^4+8400*Gamma^6*delta^5+9*Gamma^8*delta-876*Gamma^7*delta^2-9328*Gamma^6*delta^3+41520*Gamma^5*delta^4-37600*Gamma^4*delta^5+1330*Gamma^7-700*Gamma^6*delta+152*Gamma^5*delta^2+46464*Gamma^4*delta^3-120960*Gamma^3*delta^4-268800*Gamma^2*delta^5-5360*Gamma^5-11088*Gamma^4*delta+14144*Gamma^3*delta^2+20160*Gamma^2*delta^3+12800*Gamma*delta^4+64000*delta^5+800*Gamma^3+2240*Gamma^2*delta-1920*Gamma*delta^2-6400*delta^3)/(Gamma^4*(27*Gamma^4+148*Gamma^2+160)), Lambda__3 = delta}, {Lambda__2 = beta, Lambda__3 = beta}, {Lambda__2 = -(Gamma^2*alpha+4*Gamma+4*alpha)/(Gamma^2+4), Lambda__3 = alpha}, {Lambda__2 = 0, Lambda__3 = 0}


{Lambda__2 = delta, Lambda__3 = (1/5)*(-900*Gamma^7*delta^4+8400*Gamma^6*delta^5+9*Gamma^8*delta-876*Gamma^7*delta^2-9328*Gamma^6*delta^3+41520*Gamma^5*delta^4-37600*Gamma^4*delta^5+1330*Gamma^7-700*Gamma^6*delta+152*Gamma^5*delta^2+46464*Gamma^4*delta^3-120960*Gamma^3*delta^4-268800*Gamma^2*delta^5-5360*Gamma^5-11088*Gamma^4*delta+14144*Gamma^3*delta^2+20160*Gamma^2*delta^3+12800*Gamma*delta^4+64000*delta^5+800*Gamma^3+2240*Gamma^2*delta-1920*Gamma*delta^2-6400*delta^3)/(Gamma^4*(27*Gamma^4+148*Gamma^2+160))}, {Lambda__2 = beta, Lambda__3 = beta}, {Lambda__2 = alpha, Lambda__3 = -(Gamma^2*alpha+4*Gamma+4*alpha)/(Gamma^2+4)}, {Lambda__2 = 0, Lambda__3 = 0}

Consider the third solution where one solution looks simpler than the other. But after taking allvalues  we see that finally the complexity is the same although it didn't look like it. For this third solution the Lambda's are different (complex conjugates) and complex


{Lambda__2 = -(Gamma^2*alpha+4*Gamma+4*alpha)/(Gamma^2+4), Lambda__3 = alpha}

[{Lambda__2 = -(I*Gamma+2)*Gamma/(Gamma^2+4), Lambda__3 = (I*Gamma-2)*Gamma/(Gamma^2+4)}, {Lambda__2 = (I*Gamma-2)*Gamma/(Gamma^2+4), Lambda__3 = -(I*Gamma+2)*Gamma/(Gamma^2+4)}]

{Lambda__2 = alpha, Lambda__3 = -(Gamma^2*alpha+4*Gamma+4*alpha)/(Gamma^2+4)}

[{Lambda__2 = (I*Gamma-2)*Gamma/(Gamma^2+4), Lambda__3 = -(I*Gamma+2)*Gamma/(Gamma^2+4)}, {Lambda__2 = -(I*Gamma+2)*Gamma/(Gamma^2+4), Lambda__3 = (I*Gamma-2)*Gamma/(Gamma^2+4)}]


Download NondimshortAllEqns.mw


@dharr For the case where lambda__2 = lambda__3, there is always a positive root. The same analysis could also be done for the case where the two lambda's are not set equal.

Edit - actually Descartes rule of signs tells us there must be one positive root no matter what the sign of the term with coefficient 4-9*Gamma^2




2. In general I understand, but I thought my commands l1-remove(has, l1, I), l2-remove(has, l2, I), and l3-remove(has, l3, I) would return the complex solutions for the respective lambda, if any exist. Since they all return zero, I thought this meant all solutions were real. Am I wrong? 

has(...,I) is not a good test. Firstly, an expression may have some I's that cancel if the expression is simplify using evalc. more importantly, everything depends on the parameters. Sor example sqrt(1-sigma__v) might be real or complex depending on the values of sigma__v even though there are no explicit I's.


3. I am not sure I got you. In your example the two specular equations generate specular elements for a solution. Don't look at your av variable but at your ans variable (i.e., implicit form) instead. In particular, look at your third one (with RootOf):

They are not exactly the same, but the RootOf term is the same (a quadratic). In short, the difference is just a sign and the -1. Instead, I get:

Check output 16 of the new script I already attached to see more clearly.

It is not the convolutedness that surprises me, but the fact that lambda_3 is a quartic in _Z, while lambda_2 is made of some weird combination of (actually the same) quartic and with some randomly looking leading coefficients. My question: is this "asymmetry" plausible from two perfectly symmetric equations? In particular, how can the significant difference in the size of the implicit forms for lambda_2 and lambda_3 be explained?

Respone: I'm not sure which file has the relevant eq 16 in; not the last one you posted?. I expect that if you have the quadratic formula of the  quartic formula, it might look different from the quartic formula on the quadratic formula, even if they are actually the same. (I would solve the quadratic first becuase is is easier to see what the right root might be.) To see if they are the same, choose some random values of the parameters and see if they actually evaluate to the same number.

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