@mmcdara Good catch. Great that you spotted that and I appreciate all the details! I leave a reference to an old question of mine who was best answered by @Carl Love for anyone who is interested in learning more about these commands: https://www.mapleprimes.com/questions/236473-How--To-Extract-The-Numerators-Of-Equations

You write: "Maybe this correction will raise new problems because it seems the system  EqN := numer( normal~(Eqs_cal) ) has only one trivial root for the data you provide (lambda__n=0, n=1, 2, 3,  but this gives undefined values for Eqs_cal)." The data I provided are absolutely random (within the minimal requirements that gamma, sigma_e and sigma_v are all real and positive and that p is a real, positive number between 0 and 1) and I wouldn't be surprised if they produce weird solutions. That's why I think it's worthy for me to keep trying to obtain a non-calibrated analytical solution first and then perhaps do some further analysis to check whether the solution fails for some specific values of the parameters...

What would you suggest to do next?

@dharr thanks for 'connecting the dots'. I now fully understand that while the two solvers output the same solution, as @mmcdara also showed, PolynomialSystem is preferable in some cases.

You write: "had there been restrictions on lambda__2 and lambda__3, then fishing out the required solutions from all of them probably would have been easier starting from the PolynomialSystem, backsubstitute=false form, as was the case in your earlier problem". What restrictions were imposed on lambda_2 and lambda_3 in my earlier problem? Which earlier problem?

If we are referring to the same earlier problem, everything should be the same (if I recall correctly) except for the fact that before the system of three equations always collapsed to one of two equations if we allowd for two of the three lambdas to be equal. The significant simplification also occurs in the current problem, as I show in the original script in this thread (see EqN_redox part), but only if probability p=0.5. However, it's important for my problem at hand to consider p other than 0.5, so I will need to keep working on the system with all three equations. If, instead, you were referring to other restrictions on the lambdas, please tell me what you meant. 

Thank you.

@mmcdara First of all, it seems that mapleprimes.com was down for the whole day yesterday. Sorry for not being able to get back to you earlier. Voted up of course! Your analyses always go above and beyond, but was thinking of waiting a couple of days more or so to give best answer (so that more users are incentivized to contribute if they feel like).

Somehow I missed your second version of the script, which clearly answered my clarifying questions 1. and 3. already. I guess there was an issue with attachments on the mobile version of mapleprimes.

Now to your third and final version 141123_Problem_NoCorrelation_mmcdara_3. You came up with the idea of taking the limits. This is interesting.

1. First of all, I changed the values of the parameters to some reasonable numbers. As I mentioned before, they can't take random values, as you assumed. Moreover, depending on the values of such parameters, one of the three equations converges to some rational expression, while the other two always to -lambda_2 and -lambda_3: 161123_Problem_NoCorrelation_mmcdara_Limits.mw. What does it mean? 
2. Now, as far as I know, if the solution to a system of nonlinear rational equations nullifies both numerators and denominators of the equations, it means that the solution is a common factor of the original equations. Moreover, if the limit of these rational equations is finite, it means that these equations are not undefined. In other words, the fact that the limit is finite implies that the solution is not a pole (non-removable singularity) of the equations, which would otherwise make them undefined. Is this good news or not? How should I interpret such solution? Is it a valid and/or meaningful solution or is it not (since it's not unique and the fact that lambda_1 is pinned down precisely by arbitrary lambda_2s and lambda 3s is just a trivial consequence of my setup rather than an actual solution) The fact that such solution is a common factor should be because of the intrinsic symmetry in my model.

In short, I am not sure what I should make of the solution obtained.

Thanks for further clarifying on the above.

@mmcdara ok, I'll dig deeper into PolynomialSystem.

thanks for explaining all total degrees.

What about my previous two clarifying questions 1. and 3.? Both refer to *your* second script.

(ignore the 'redox' snippet for p=.5, as I don't even use it)

@mmcdara thanks, as usual, for the great details. I understand that for each lambda the max exponent is just 3, but then there are convoluted products of them. I am not sure I fully grasp the interpretation of AllTotalDegrees (is it just the sum of the degrees for each lambda term?). In any case, you convinced me that solve() can't tackle my Eqs, but I don't plan to solve them numerically. I totally understand that a closed-form solution is a really hard sell here, but if I were to solve them numerically I would rather just stick to the numerical solutions I already obtained for more complicated versions of my problem...

1. In the second script of yours (where you try to solve for the numerators), if my understanding of it is correct, you show me that lambda_2 and lambda_3 are basically free to choose and lambda_1 can be pinned down from them. Is it true for all positive combinations of lambda_2 and lambda_3 that the roots of the numerators are also those of the denominators of my equations?
2. If you look at my original script with SolveTools:-PolynomialSystem, a constraint appears (check Solution[1][2]). What can I make of it? How to use Solution[1][2] to find the roots of the polynomial Solution[1][1]? Would backsubstituting some lambda_2s and lambda_3s on Solution[1][1] lead me to the same expression for lambda_1 that you found by using solve()? More in general, you didn't discuss much the interpretation of what I obtained from using SolveTools:-PolynomialSystem. 
3. Finally, in your answer you mention three alternative ways to do the 'verification' in my case, but none of these work with my solution structure. Since you write: "... or something else depending on the structure of sol", what way would work then?

Always appreciate your help.

@C_R 

If I understood your script correctly, you basically show me (as you mentioned) that solve() gets confused. Thanks for pointing that out. Then, what's the correct way to make solve() work with my assumptions? They are not that many actually, just 8 (and mostly very simple positivity requirements).

You also write "I also feel a bit uneasy about the indexed subscripts because of typos that may not be visible." I understand, but what do you suggest?

Thanks a lot!

Both answers clarify my doubts (thank you @acer), but @Carl Love 's one goes above and beyond with all the refreshers on the mathematical notions and the digressions. Hence the best answer.

@Carl Love thanks.

1. How to do this replacement within fsolve() command?
2. Since you brought up convergence...what's the standard (implicit) convergence of fsolve? How to output the error in case I want to have a better idea of the magnitudes compared to the numerical solutions given (apart from eval(equations,solutions))?

@tomleslie amazing answer! I am learning a lot. Thank you for all the details. I see that you provided fsolve() with specified starting points in all calibrations, not just those with "missing" values. Is this a good practice more in general? Should I always provide starting guesses or only if I encounter issues from some iteration onwards?

@Carl Love thank you for teaching me about this.

@Carl Love oh okay I understand now. Thanks!

@Carl Love the 1 plot per file solution (2nd option) seems like what I need. 

Since I have 8 consecutive plots:-display( seq( plot( [$ 20], series1of8[1..20,j], color=cols1[j]), j=1..3)); (one for each one of the 8 series), does it mean I first need to combine them all together as a plot sequence? (Is this what was implicit in Plots:= <seq(plot([$20], 'color'= cols1[j]), j= 1..3)>: ?)

@tomleslie in your first comment/answer you wrote "you have to ensure that all the data is available for all the plots whihc you want!"

I found missing data for data series 7 and 8: someplots_cal7_cal8_bounds_issue.mw

I found that series 7 stops at run 33, while series 8 at run 52. Of course if I only plot the first 20 elements or so there's no issue...but why do I have substantially less data for these two calibrations?

Separately, I noticed that your fsolve command for lam did not include the positive range for the lambdas, e.g. lam := fsolve(eval(MyEqs_cal, ncal3)) instead of lam := fsolve(eval(MyEqs_cal, ncal3), {lambda__1 = 0 .. infinity, lambda__2 = 0 .. infinity, lambda__3 = 0 .. infinity});...did you do this deliberately? For what reason?

@dharr thank you so much for the efforts. Interesting. I didn't think about it...I am learning a lot in this thread.

I was trying to look at a much simpler system: Analysis_quintic_solution.mw. Is there a way to find a closed-form real and positive solution to the quintic (with quadratic term missing) I find when I don't make any assumptions on the positivity of the two variables? Descartes rule of sign suggests that one such solution may be attainable, but how to find it explicitly? In the script are my attempts.

The system is unrelated to the one discussed over this thread, so feel free to ignore this last analysis, if you wish. Thank you for your help so far.

@dharr thank you for the exhaustive analysis. I checked a few times. EqN_rhou in the startup code seems to be the same I fed into the solver. However, as long as the solution solves the quadratic equation which was output by the solver, then it must also be a solution to the system...am I wrong?

I re-run that single calibration (took just 1 minute with the reduced 2-equations system). The quadratic equation produced (solver worksheet) is exactly the same as the one loaded in the analysis worksheet, but strangely the size of the solution differs slightly between the two worksheets. Is it the command sols_rhou := eval(SOLUTION_LAMBDAS_cal_rhou, '='tilde(%, {rho__u, lambda__2, lambda__3})) in the analysis worksheet that slightly reduces the length?

SOLVER WORKSHEET:

ANALYSIS WORKSHEET:

In any case, I reproduced the analysis calculations in the solver worksheets, directly on the solution output by the solver, and the lambdas still do not solve EqN_rhou. Why? How come that they solve the quadratic polynomial but not the original system?