dharr

Dr. David Harrington

5712 Reputation

21 Badges

20 years, 344 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

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

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

RemoveUnits.mw

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

restart

local gamma;

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;
isolate(%,lambda__1);
Eqq1:=simplify(eval(Eq1,%))*sigma__d/sigma__v;
solve(Eqq1,Lambda__1);

Lambda__1 = lambda__1*sigma__d/sigma__v

lambda__1 = Lambda__1*sigma__v/sigma__d

(-5*Lambda__1^3+2*Lambda__1)/(5*Lambda__1^2+2)

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};
ss:=solve(%,{lambda__2,lambda__3});
Eqq2:=numer(normal(eval(Eq2,ss),expanded));
Eqq3:=numer(normal(eval(Eq3,ss),expanded));

{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}

-(5*gamma^2*Lambda__2^3*sigma__d^2*sigma__v^2+8*gamma^2*Lambda__2^2*Lambda__3*sigma__d^2*sigma__v^2+5*gamma^2*Lambda__2*Lambda__3^2*sigma__d^2*sigma__v^2-2*gamma^2*Lambda__3*sigma__d^2*sigma__v^2+40*gamma*Lambda__2^3*Lambda__3*sigma__d*sigma__v+32*gamma*Lambda__2^2*Lambda__3^2*sigma__d*sigma__v-4*gamma*Lambda__2*Lambda__3*sigma__d*sigma__v-4*gamma*Lambda__3^2*sigma__d*sigma__v+80*Lambda__2^3*Lambda__3^2-8*Lambda__2*Lambda__3^2)*sigma__v

-(5*gamma^2*Lambda__2^2*Lambda__3*sigma__d^2*sigma__v^2+8*gamma^2*Lambda__2*Lambda__3^2*sigma__d^2*sigma__v^2+5*gamma^2*Lambda__3^3*sigma__d^2*sigma__v^2-2*gamma^2*Lambda__2*sigma__d^2*sigma__v^2+32*gamma*Lambda__2^2*Lambda__3^2*sigma__d*sigma__v+40*gamma*Lambda__2*Lambda__3^3*sigma__d*sigma__v-4*gamma*Lambda__2^2*sigma__d*sigma__v-4*gamma*Lambda__2*Lambda__3*sigma__d*sigma__v+80*Lambda__2^2*Lambda__3^3-8*Lambda__2^2*Lambda__3)*sigma__v

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;
sg:=isolate(%,gamma);
Eqqq2:=simplify(eval(Eqq2/sigma__v,sg));
Eqqq3:=simplify(eval(Eqq3/sigma__v,sg));

Gamma = gamma*sigma__d*sigma__v

gamma = Gamma/(sigma__d*sigma__v)

(-5*Gamma^2*Lambda__2-32*Gamma*Lambda__2^2-80*Lambda__2^3+4*Gamma+8*Lambda__2)*Lambda__3^2-8*(Gamma*Lambda__2^2+5*Lambda__2^3-(1/4)*Gamma-(1/2)*Lambda__2)*Gamma*Lambda__3-5*Gamma^2*Lambda__2^3

(-5*Gamma^2*Lambda__3-32*Gamma*Lambda__3^2-80*Lambda__3^3+4*Gamma+8*Lambda__3)*Lambda__2^2-8*(Gamma*Lambda__3^2+5*Lambda__3^3-(1/4)*Gamma-(1/2)*Lambda__3)*Gamma*Lambda__2-5*Gamma^2*Lambda__3^3

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

ans23:=SolveTools:-PolynomialSystem({Eqqq2,Eqqq3},[Lambda__2,Lambda__3],engine=triade);

{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}

ans32:=SolveTools:-PolynomialSystem({Eqqq2,Eqqq3},[Lambda__3,Lambda__2],engine=triade);

{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

ans23[3];
map(simplify,[allvalues(%)]);
ans32[3];
map(simplify,[allvalues(%)]);

{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)}]

NULL

Download NondimshortAllEqns.mw

 

3 4 5 6 7 8 9 Last Page 5 of 58