@acer thank you. I didn't know about RootFinding:-Parametric.

1. Why the plot doesn't change if I also include -1<rho and rho<1 among the equations to solve? I understood that regions 1, 4 and 5 have 0 solutions anyway but if I add the constraints on rho I should expect 2 regions instead of 5, right?
2. I don't understand SampleSolutions(m,2)=0.56 and SampleSolutions(m,3)=0.51 (even after reading help page). How are these numeric values found?
3. Looking at the big picture, how to reconcile this CellPlot with A) my plot3d of Eq and B) @mmcdara Sturm's analysis. I am having a hard time putting all together.

Thanks again for following up on this thread.

@mmcdara I am finally back to this. Since it's my fault not having specified Gamma>0 as a pre-requisite of this question, I "fixed" your analysis worksheet myself to consider only Gamma>0: rho-analysis_mmcdara_Gammapositive_MaPal.mw.

1. Would you please double check whether I did everything correctly? I highlighted in light blue my changes.
2. What's the conclusion of Sturm's analysis? My understanding: for rho=0.6 there are 4 distinct real roots when 0<Gamma<max(4 roots-green bar), 2 real roots when Gamma is in yellow bar, and 2 real roots when Gamma is in red bar. Am I correct?
3. We now confirmed the existence of 4 real Lambda roots when Gamma values are in the green bar. How to isolate them? That is, what is the closed form of these 4 roots? How do I formally prove that just one out of four is positive (since this is what plot3d(Eq) for Gamma>0 and -1<rho<+1 suggests, see plot attached below)? 
4. What about the two real roots in the yellow rectangle and the two real roots in the red rectangle? Why plot3d shows just 4 roots for any Gamma in the range 0<Gamma<10? 
   Thank you for addressing my remaining doubts @mmcdara .

@acer sure, no problem. Or perhaps there's an easier way to infer the signs of these partial derivatives given that I know the signs of l1, l2 in their specified form and of their partial derivatives wrt to their variables Gamma_1 and Gamma_2 (thanks to @dharr's nondimGamma1Gamma2.mw, see the plots at the bottom). 

Note also that Gamma_1=gamma*sigma_d*sigma_v1>0 and Gamma_2=gamma*sigma_d*sigma_v2>0 because all the underlying parameters are strictly positive. These are the partial derivatives for L1, L2, B1, B2, A1, A2, A2S and A1S when l1 and l2 are left unspecified:

Download derivatives_signs.mw

@dharr and @acer thank you for your comments. Do you agree that both approaches are equivalent?

Moreover, doubling the digits doesn't solve the jumps...I also just noticed that the partial derivatives of L1 and L2, while still convoluted, are significantly simpler than those for B1, B2, A1, A2, A2S and A1S. Maple takes forever to find the derivatives for the latter 6 functions (the diff() step right before plot3d(preprocess())...). For example, figuring out just dB1dSd and its 3D plot is taking Maple more than 45 minutes...

So perhaps "manipulation to a nicer form" for both l1 (the RootOf()) and l2 (rational function of greatest l1 root) is needed here, but I understand that it is not easy. @dharr you mention continued fractions...is it virtually equivalent to an approximation that is accurate for an upper bounded range of positive Gamma_1 and Gamma_2 values? @acer do you have some ideas on how to manipulate or, if easier, directly approximate the bivariates l1 and l2 (and measure the accuracy)?  

If needed, I can spin off the question about the manipulation/approximation of l1 and l2 into a separate question (since it's not about non-dimensionalization anymore) but moderators probably don't want that.

@dharr that makes sense, it's just a scaling. Just replacing T with Gamma__2/Gamma__2 confirmed what you say. Thank you.

In a comment above you wrote:
"The problem may be close roots not being properly resolved by the default numerical routines used by plot3d. Also, the order of the RootOfs if there is more than one positive root is to have the smallest one first. The solution is to be careful about always choosing the most positive root for Lambda__1."

I think a similar issue emerges when I try to plot3d partial derivatives of some functions of Lambda__1(Gamma__1,Gamma__2)=RootOf(degree-10-polynomial) and Lambda__2(Gamma__1,Gamma__2) (where Lambda__2 should be interpreted as a rational function of the greatest root of Lambda__1 as in your nondimGamma1Gamma2.mw above):
derivatives_jumps.mw

In the attached script derivatives_jumps.mw you see jumps/discontinuities emerging in the derivatives of L1 and L2, seemingly for Gamma__1=Gamma__2~1...I am not sure whether these jumps really exist though, given that all parameters are strictly positive and in nondimGamma1Gamma2.mw above you found the partial derivatives of Lambda__1 and Lambda__2 wrt to Gamma__1 and Gamma__2 to be either strictly positive of strictly negative...

then how to have a "better"/smoother numerical evaluation in these more complicated partial derivative cases and for the ranges Gamma__1=0.01..5,Gamma__2=0.01..5? Quoting you, I know I should "make a function that reliably returns the numeric value of the largest root of Lambda_1 first and then change the RootOf form just to Lambda__1 (similar to if we had solved without backsubstitution) and use the numerical solution for Lambda__1 to find Lambda__2" but where to exactly implement this in derivatives_jumps.mw without messing up the function assignments? Perhaps when I specify the forms for l1 and l2 and their respective aliases? Moreover, for the derivatives I guess I would need the denoms of Lambda__1 and Lambda__2...

@dharr great, thanks. It's good to see that the positive Lambda_1 and Lamnda_2 are unique. By the way, how did you explicitly work out the nondims and corresponding Eqq2 and Eqq1 (output (3) in nondimGamma1Gamma2.mw) using Hubert's theory? Is it an adaptation of Case_with_radicals.mw (but without the additional equation)? I can't get to the same Eqq2 and Eqq1 using that...

@dharr thank you for getting back to me on this.

In your worksheet you mention "Lambda__1 is the largest root of a degree 10 polynomial" and while answering to my point 1. in the comment above you mention "If the RootOf for Lambda__1 is consistently interpreted as the most positive real root, and this value is used to calculate Lambda__2 then both Lambda__1 and Lambda__2 are consistently positive".

What's the meaning of "most positive" in this context? and above all:

• Given strictly positive Gamma_1 and Gamma_2, is the positive Lambda_1 (and thus Lambda_2) we use unique? In other words, can I prove that the remaining 9 roots of the RootOf(_Z^10-polynomial) are either negative or not consistently positive for strictly positive Gamma_1 and Gamma_2?

@dharr great catch for spotting the mismatching notation! 

@dharr I am back to this thread. Your conclusion here was that numerical analysis showed that Lambda_1 is actually strictly positive (for Sigma = 0.01 .. 5) and Lambda_2 as well (for Sigma = 0.03 .. 5). For both Lambda_1 and Lambda_2 the plots are not just strictly positive but, importantly, they also show no weird cliffs as in the plots for the analytical/symbolic-form analysis.

You tried two different non-dimensionalizations and concluded that in these two cases analytical/symbolic-form analysis leads to a Lambda_1 with some weird cliffs and a Lambda_2 that is NOT strictly positive, i.e., the surface has a boundary. In other words, you couldn't find a symbolic form for Lambda_1 and Lambda_2 that is consistent with the numerical solutions (at least for the first root of the degree-10 polynomial and for T and Sigma both not too close to 0).

Questions :

1. (Assuming that it's tough to guess a non-dimensionalization Lambda_1 and Lambda_2 leading to two symbolic-forms whose plots are consistent with the well-behaved numerical solutions) is it possible to reverse-engineer a (even approximative) symbolic form from the numerical solutions?
2. Even if we managed to reverse-engineer such symbolic form, would I encounter issues if I wanted to take partial derivatives of Lambda_1 and Lambda_2 wrt to sigma_v1 and sigma_v2 (since T and Sigma "entangles" both sigma_v1 and sigma_v2)? Specifically, you transitioned from a Gamma_1=sigma_v1*sigma_d*gamma and Gamma_2=sigma_v2*sigma_d*gamma non-dimensionalization to a T=sigma_v2/sigma_v1 and Sigma=gamma*sigma_d*sqrt(sigma_v1^2+sigma_v2^2) non-dimensionalization just to get rid of the square roots, but wouldn't be more "intuitive" to work with Gamma_1 and Gamma_2? (it's easier to interpret a 3D plot with Gamma_1 and Gamma_2 as x- and y- axes)
3. Perhaps a stupid question: worst case scenario, are there ways to calculate partial derivatives when no analytical expressions are available? 

@mmcdara I am not sure I know how to make use of the previous answer to reply this question, which is very different. It would take me some time to go through all the details but, if I understood correctly, you prove that for Gamma>0 and -1<rho<+1 there is only one real and positive root. (Basically you provided a very detailed and formal analysis, which I really appreciate, of the existence of the unique positive root shown simply with plot3d for the ranges of Gamma and rho above, or am I wrong? I need to read more about Sturm's method).

In this question I am just interested in the sign of the partial derivatives of such unique positive root.

@mmcdara deserves best answer for the impressive analysis. It'll take me some time to go through it and digest it. Thanks a lot!

@acer thanks. I think I also figured out query 1). I think I just made a bad choice about my values for Gamma and rho when I tested the zeros of my function.

Since I am only interested in Gamma>0, s2 is not relevant to me and I can just play around with Gamma values > or < than the Gamma__0 given by s3, for which there will always be just one positive root. This is then consistent with the unique positive zero given by implicitplot3d with Explore().  

Download always_just_one_root.mw

@dharr wow that's great! and it seems that apart from the region for small Sigma and T the surfaces are quite smooth for both lambda_1 and lambda_2...

@dharr thank you so much for all the effort!! I really appreciate your help. 

I am a bit confused though, perhaps I need to review your three worksheets again...Your conclusion is that we can't work out an expression for the positive portion of lambda_2 because we can't exactly identify its boundaries? 

@C_R thank you for the detailed answer and @mmcdara for your comments!