@MaPal93 

Let me be clear about my goal. I am looking for a (i) closed-form, (ii) real, and (iii) positive solution. Even just one solution satisfying all three constraints would be fantastic.

Now, it seems that sol[1] (2 equations in the lambdas) and sol[5] (3 equations in the lambdas) do not look weird and are the most promising for further analysis.

sol[2..5] all involve solutions of high degree polynomials, which in general do not have solutions in terms of radicals. So when you put your parameters in, you will not likely get symbolic solutions. So sol[1], which is surprisingly just a quadratic is your best chance. 

However, for the case with all parameters =1, an exhaustive analysis proves that there are no positive solutions for any of the 5 cases:
triade.mw

My questions:

1) How did you find the unique and real but negative lambdas from sol[5]?

In each case, you solve the third equation for lambda__1, substitute it into the second equation and solve it for lambda_2 and then substitute these solutions into the first equation and solve for lambda__3 (this is the backsubstitution on the triangularized system). See details in the file above. 

2) You showed me one example in which also sol[1] produces a real but negative solution. In the same comment you mentioned "but real and positive is harder". Why is it the case? If mathematically real and positive solutions are not ruled out a priori, is there a more systematic way to find them? (even just a loop which tries many substitutions and stops as soon as we find one would do, I guess?)

I did it in detail on the attached. With parameters it will be harder, but some things can help in this sort of analysis, Descartes rule of signs etc.

3) It's not clear to me the benefit of dividing out the common factor in the system with all parameters normalized to 1 (since we already found solutions and, factorized or not, these solutions would be the same). The benefit would be clear if such removable common factor existed in the uncalibrated equations EqN := ((numer@evala@:-Norm@numer)'tilde'@eval)(MyEqs). How to verify this possibility?

I agree. I assume factorizing before solving speeds up the solution process. Factorizing may not work with parameters but you can try the factor() command to see.

4) Surely a naive question, but would a solution found as in 2) solve also the uncalibrated system? That is, would simplify(eval(EqN,sol)) still give me 0 for EqN := ((numer@evala@:-Norm@numer)'tilde'@eval)(MyEqs) instead of EqN := ((numer@evala@:-Norm@numer)'tilde'@eval)(MyEqs,P='tilde'1)?

Seems very unlikely.

5) Related to 4). My end goal is to study how lambda_1, lambda_2, and lambda_3 vary with my parameters. Is it legit to pick a real and positive lambda_2 and lambda_3 (found as in 2)) and plug them back into uncalibrated EqN[1] (and solve for lambda_1), then pick lambda_1 and lambda_3 and plug them back into uncalibrated EqN[2] (and solve for lambda_2), finally, pick lambda_1 and lambda_2 and plug them back into uncalibrated EqN[3] (and solve for lambda_3)? See my original script 160523_stylized_triade.mw: the first equation is for lambda_1, second for lambda_2, and third for lambda_3.

I'm not sure I understand this, but in general the solution for all parameters=1 is not going to generalize to the case with parameters.

Really thank you for helping me out with this.

You're welcome. There are some interesting symmetries that might help you, For example EqN[2] and EqN[3] are the same if lambda__1 and lambda__2 are exchanged (and this was true for sol[1].) Is this expected?

@MaPal93 I looked again at the solutions, and many will be roots of high-degree poynomials that have no closed form. Given that, I don't think backsubstitute=true will help. You can analyse to find when there might be real or positive roots. For example, the polynomial in the fifth solution has only one real solution, which is negative, so lambda_1 = lambda_2 =-0.41 and lambda_3 = -0.82. No positive solutions. Here is the solution with the common factor removed

factored.mw

@dharr For the EqN for all parameters 1, the three equations have a common factor, that can be divided out so that you could work on solving a simpler system. If this still works with parameters in, that could afford a considerable simplification. 

@MaPal93 In general, I expect the solvability to depend on the parameters. This seems like a very complicated problem with too many parameters - how about solving a simpler version of the problem with a smaller set of parameters first? But if I were doing this problem I would step back and analyze the matrices. Are lambda's eigenvalues? if so a lot can be said about their possible values depending on the structure of the matrix of which they are eigenvalues, e.g., we know many classes of matrices that are guaranteed to have one or more real eigenvalues.

@dharr For solution 1 you can find real solutions as follows (but real and positive is harder). 

Download EqNsol.mw

@MaPal93 Some answers:

1) what is exactly Eqn? could you load my solution instead of hard coding expressions?

I definitely did not reenter this by hand! This was extracted from your worksheet by opening a new worksheet with the same Maple engine - try EqN; in your own worksheet to check.

2) is this system indeterminate (infinitely many solutions)?

Yes.

3) would it be preferred if I chose backsubstitute = true?

Probably the answer is too complicated to analyse, and maybe too slow to achieve - but you could try.

4) what happened to my parameters P := indets(MyEqs, name) minus MyVars ? I was hoping to obtain solutions as functions of those.

I just looked at what you were passing to the solver, i.e., EqN. Looks like in the line EqN := ((`~`[((numer@evala)@(:-Norm))@numer])@eval)(MyEqs, `=`~(P, 1)) you set all the parameters to 1. 

5) is this solving approach (a) reliable and (b) preferable to engine=groebner or standard solve()?

I don't see why the answer wouldn't be reliable. It worked for the case I tried at random. I don't know much about different engines; solve probably passes to one or other of these.

6) how to enforce RealDomain?

If you want only real solutions, try RealTriangularize in the RegularChains package, which was a recent suggestion somewhere here on Mapleprimes. But I really don't know much about solving polynomial equations beyond the bivariate case.

@MaPal93 Here's my understanding. You chose backsubstitute=false, so you can carry out the backsubstitution yourself, e.g.

Download EqNsol.mw

@occipitalbaker Thanks.

I did something similar before here but for the case where the walk had to start and end at a particular vertex. (Perhaps that is the same problem?) There is a significant difference between just finding the length and finding the walk itself. Are you OK with just the length?

@Axel Vogt Vote up. The second one can be done without a change of variable.

Download int.mw

@rockyicer You're welcome. If you want to convert the .lib file to .mla, you can use the LibraryTools:-ConvertVersion command.

Yes, development is always going on - new integration methods were added in the latest (2023) version, see

https://www.maplesoft.com/support/help/Maple/view.aspx?path=updates%2fMaple2023%2fAdvancedMath

@mmcdara If you try Eigenvalues(SQRT),Eigenvalues(Q_L); you will see that two of the eigenvalues are negative of what they should be, i.e., the method has effectively chosen the wrong sqrt. I'm guessing that can be solved by a different starting guess, but not sure how exactly that could be done. Note that for the Eigenvalue method, you can work completely in reals (and more efficiently) by telling it C is symmetric:
 

LEC := proc(C)
  local L, E;
  L, E := LinearAlgebra:-Eigenvectors(Matrix(evalf(C),shape=symmetric)):
  E . DiagonalMatrix(sqrt~(L)) . E^(-1)
end proc:

But if you use a numerical version of my routine above (with shape = symmetric) then you can avoid the numerical inverse (just a transpose now) and it might be more efficient (but the normalization will slow it down)

Msqrtf:=proc(A::Matrix(square,datatype=float))
  local evs,evals,evecs,U,i,j,k,g,Lambda;
  uses LinearAlgebra;
  evs:=Eigenvectors(Matrix(A,shape=symmetric),output='list');
  evals:=table();evecs:=table();
  j:=0;
  for i in evs do
    if i[2] = 1 then
      j:=j+1;
      evals[j]:=i[1];
      evecs[j]:=Normalize(i[3][],'Euclidean');
    else
      g:=GramSchmidt(i[3],conjugate=false,normalized=true);
      for k to i[2] do
        j:=j+1;
        evals[j]:=i[1];
        evecs[j]:=g[k];
      end do;
    end if;
  end do;
  Lambda:=DiagonalMatrix(map(sqrt,convert(evals,list)));
  U:=convert(convert(evecs,list),Matrix);
  U.Lambda.U^%T;
end proc:

@sand15 Thanks - I remembered that and had a quick look for it but didn't find it. The emphasis there was on the numerical case, but here on symbolic. For symbolic use, the method here requires that the exact eigenvalues can be found. Here even though the matrices are large, the eigenvalues and eigenvectors are easily found (for M1..M3, the eigenvalues are actually integers). The method here avoids finding a matrix inverse, or more accurately the inverse is just the hermitian transpose so is easy to compute.

My recollection is that Maple's method can suceed in more cases, at the expense of being slow. The solution is presumably for MatrixPower and/or MatrixFunction to do more preselection of different algorithms depending on the matrix type.

convert(expr,elementary) is intended to convert hypergeoms etc to elementary functions, but doesn't seem to work on these cases.