@acer thanks. I successfully installed the package and ran SolveEquations() with and without the AllSolutions option. Below my results:

DirectSearch_nonlinear_system_calibration.mw

1) Is my implementation of SolveEquations() correct?

2) How do I interpret my results in the two cases?

While I don't have ballpark figures for my mu(s) and lambda(s), I would expect them to be larger given my calibration.

3) How to (hopefully in a minimal way) alter my calibration to get larger solutions? Alternatively, is there a way to just output the largest of all solutions?

Thanks

@dharr and @mmcdara thanks for your comments. 

Unfortunately, I can't find a solution with my calibration of free parameters as below:

NOT_WORKING_nonlinear_system_calibration.mw

However, with an alternative but unrealistic calibration, I can find more than one solution (check the loop):

WORKING_nonlinear_system_calibration.mw

How can I modify my script to obtain some results with my calibration? If roots are really not available, I am okay with using optimisers. However, I had issues installing DirectSearch as outlined in here:

Evaluate min and max in Excel - MaplePrimes

Thank you.

@dharr perfect. I adapted your routine to my script. All's good now.

Thank you again for the continued support.

@dharr thank you. Then, since I need to stick to a numerical solution, I guess it makes sense to look into the convergence of my solutions. I was thinking of verifying the convergence by computing basic statistics (mean, std, and so on...) of my 8 solutions ( the 6 lambda(s) and 2 mu(s) ) for an arbitrary number of fsolve() runs.

• Is this a reasonable approach? (If not, please let me know your recommended alternative).

Ideally, for all runs combined, I would like to see a very low std (so that I am confident that the algorithm used by fsolve() isn't simply stuck in a local minimum) and use the 8 means (computed over all the number of runs, which can be ~100) for each of my variables to perform further analysis.

In the example below, I start with just 4 runs but I would like to choose the number of runs arbitrarily and compute said statistics in the most efficient way:

num_sol_generic_stats.mw

First of all, I encounter a problem with my data structure. How can I map the fsolve() outputs to a matrix so that each row corresponds to a run and each column to a solution in the set of solutions returned by the corresponding fsolve() run? For example, for 100 runs I would expect a results matrix of dimensions 100 rows x 8 columns. The reason why I'd like to have a matrix M rather than a dataframe, table or similar is that the means of my solutions would simply be the means of all 100 values for each column (related to this: how do I extract the numerical value of a solution given, e.g., as mu_jk = -0.2279422049?), which are simply given by the command Mean(M). 

• Moreover, I'd like to have a procedure which allows me to flexibly amend the values of the free parameters in Eqp ([q__0jk = 0, q__0ki = 0, Sigma__0jk = 1, Sigma__0ki = 1, rho__XX__23 = 0, sigma__uji = 1, sigma__ujk = 1, sigma__uki = 1]) to check how the nonlinear system responds.
• Finally, given my values / means for the 6 lambda(s) and the 2 mu(s), I'd need to pin down the beta(s) and the alpha(s) defined at the beginning of the script.

Thank you again for your extensive supportive. I very much appreciate Maple responsive community!

@dharr thank you for the updated file and the explanation. I am quite new to Maple, would you mind explaining:

• what the ^+ does at the SMI step (4th line)
• why in natlq2 and natlq3 steps SMI2 and SMI3 are not transposed (see your initial comment) before being multiplied by Y

Thank you again for your continued support.

@dharr as you see, even if I don't compute the inverse explicitly, the execution still gets stuck, this time at the LinearSolve() step.

@dharr thank you for the reply. I am quite new to Maple and I am not sure I am following your script. How do I avoid finding out explicitly the inverse in my case?

Please see below the full script which shows the follow-up steps to find, eventually, natlq_1 and natlq_2. Thanks!

matrix_inverse_2.mw

@Carl Love (and @acer and whoever else can help)

As you can see in my screenshot below, I encounter a "Kernel Commection Lost" error whenever I try to solve the system analytically, according to your suggested routine. Therefore, I decided to opt for numerical approaches for the time being. (However, I haven't given up on finding analytical solutions).

The values I get using fsolve() are very tiny, which does not make sense for the problem I am trying to solve.

1. How do I make sure that the numerical solution is accurate?
2. How do I check whether there is a clear convergence in the numerical solution (error is virtually zero)?
3. How do I alter fsolve() to find:
   • More solutions
   • All solutions
   • A relevant set of solutions (e.g. the largest solutions among all or those satisfying some specific chosen properties)

The attached picture displays my (too tiny) lambdas solutions. The free parameters in the system, i.e. those variables excluded from the set of variables lambdas I am trying to solve for, were calibrated according to some logic (either set to 0 or to 1). Please note that I do not have any ballpark figure for the (positive) values these lambdas should take. 

@Carl Love thank you again for the details. I followed your suggestions. kernelopts(numcpus) returns 256. Unfortunately, I still lose the connection to the Kernel after hours of running...

You previously mentioned that you wouldn't be surprised if the solver procedure takes several days. Does it mean that such error is somehow unexpected or avoidable for the problem at hand?

@Carl Love thank you for your kind support and continued effort. Your approach seemed promising, but my execution has been running for more than 7 hours now. Given your experience, is it reasonable that the execution is stuck at the evala@:-Norm step for such a long time?

After about an hour of trying with engine=groebner I switched to engine=triade and re-started the execution, as both you and @acer seemed to prioritise it.

I use Linux. I anticipate that I do not have much experience with parallel processing. For both engine=groebner and engine=triade, I noticed that the memory is not deployed beyond 30% of its capacity. Throughout all these hours, memory usage has been mostly sable at ~310GB of 1.1TB (Cache ~774GB). "Swap" 's pie chart shows 4.1GB of 4.1GB (100%). Some (but few) CPUs are working at 100% of their capacity, while most of the others are not engaged at all (0%) and just a few inefficiently engaged (<40%). The command kernelopts(numcpus) should/can be run only once I stop executing the block right?

Do you guys have further insights, given the above?

Thank you again.

@acer ok, I understand now. I am a new member of Maple community, so thanks for clarifying.

I did find some answers of yours to questions related to multicore processing and parallel computing, but I could not find anything I could direclty implement to my script. Would you kindly look into it?

Thank you all for the continued support!

@acer 

Why was it removed? Rather than a follow-up, the problem I now have has changed and I thought (and still think) it deserves another question.

@Carl Love

After some simplifying assumptions about the parameters and manipulations (exclusion of redundancies), I re-wrote my nonlinear system as shown in the second last document block (the first solve() out of the two consecutive solve()) in the following:

newsystem.mw

The number of variables are still 6, but also the equations are now 6. While the system is not overdetermined anymore (thus not almost certainly inconsistent anymore), the execution of solve() still gets stuck. The equations are all algebraic but still quite complicated. I now wonder how I can leverage parallel computing on Maple to find an analytical solution (before trying numerical procedure through fsolve()).

My machine has 128 cores (256 threads) and over 1TB of RAM. GPU is an NVIDIA Tesla A100. What's the most effective way on Maple to maximise the usage of such computational power? How can I implement solve() with the best parallel computing procedure for my problem at hand?

Thank you for the continued support!

@Carl Love Thanks for trying it independently. I encounter the same issues.

I hope someone else can help me thinking of some alternative ways to check consistency and, if consistent, if any (even just a single) solution can be found. I can't vote you up because I don't have enough reputation.

@Carl Love thank you for the very detailed and clear answer. Much appreciated. 

However, I get stuck at step 3.: EqP:= eval(Eq, indets(EqN, name)=~ 1):

I get the error "Error, numeric exception: division by zero" regardless of the simple exact value I choose (1, 0 or anything else). What's a workaround for this? I tried to follow the instructions in the following:

division by zero - Maple Help (maplesoft.com)

and by redefining how to handle exceptions does not really help me: EqP's value ("View EqP" from the "Variables" palette) includes some undefined, which I think will corrupt the procedure to test whether my system is inconsistent.

Would you mind directly implementing your check procedure on my script and sending it back to me? I attach here my script (cleared of the massive outputs) with the Eq and V sets for you to execute.

systemhasolutions.mw

Thank you for your continued support.