1

I need to solve a $25\times 25$ very non-linear system of equations, I have very superficial knowledge of numerical methods, and I've been struggling in finding a real solution to my system.

For what I know, in order to find a solution first initial guess values are crucial, but I don't really know how to get good initial guesses for my system since the "most" non-linear equations avoid closed form expressions for some variables, then I'd be very grateful if you may refer me to some literature/content (anything) that may help me understand how to get good initial guess values for my problem.

What I've done so far (before posting the question) is writing the system $F(\boldsymbol{x})$ in a function in Matlab (the system is $25\times 25$), and then try to solve for $F(\boldsymbol{x})=0$ using fsolve() and using 1 as initial guess for all the variables. The result is a very long time of computation, and after reducing (by augmenting considerably the tolerance) the computation time, in order to at least check where the solver is going, I get as a (inaccurate) solution a non-real valued vector (non-zero imaginary part).

I can reduce $n$ to a problem of around $20\times 20$ (very cumbersome) equations by substituting variables with closed form expressions. Uniqueness of solution I would not be very sure, but I'd dare to say that aside from possible trivial solutions, there's an unique solution, or at least a real-valued solution that satisfies some restrictions (like being positive for some variables or similar) would work fine.

Besides, most non-linearities arise from negative and fractional exponents, sums raised to such exponents. Not sure if it helps but there's a sample of three equations:

$$ a\left( \frac{bx_{4}}{c\left(x_{5}^{1-\phi}x_{3}+x_{6}^{1-\phi}x_{3}\right)^{1/(2-\phi)}} \right)^{\frac{\phi}{1-\phi}} = x_{7}x_{8}-dx_{9}$$

$$ dx_{4} = x_{9}\left(x_{2}\left(1-gx_{10}^{1+\theta}\right)\right)^{-1}\left(x_{9}\left[1-hx_{10}^{1+\theta}\right]-kx_{2}x_{10}^{\theta}\right) $$

$$ x_1 = x_{3}\left(1+l\left[ \frac{x_{3}}{x_{11}}\right]^2+m\left[ \frac{x_{11}}{x_{3}}\right]^2-2\sqrt{x_{15}x_{16}}\right) + x_{12}\left(1+l\left[ \frac{x_{12}}{x_{11}}\right]^2+m\left[ \frac{x_{11}}{x_{12}}\right]^2-2\sqrt{x_{15}x_{16}}\right) + qx_{13} + rx_{14} $$

Where $x_i$ are variables (most with restrictions of being strictly positive, aside of course from the other equations), latin alphabet letters are positive constants, and greek letters are also given constant between ranges: $0<\phi<1$, $\theta>1$ (of course when I'm solving the system I know the values of all those constants).

Another possibly relevant detail is that I emailed the author of a paper that solves a system with similar non-linearities (since what this problem is about is finding a steady state), he didn't give me much detail but mentioned that he was using "homotopy algorithms".

Thanks in advance!

PS: In addition, maybe if this question is too "engineery" I'd be grateful for a recommendation of a forum where it would fit best. Thanks.

manifold
  • 417
  • 1
    Depending (very much!) on the specific nonlinearities, you might be able to take linear approximations to the functions and then solve the corresponding linear system to arrive at an initial guess. – Greg Martin Jul 14 '21 at 00:28
  • 2
    This greatly depends on the specific problem, unfortunately. – Thomas Andrews Jul 14 '21 at 00:29
  • Thanks for both responses! @GregMartin approximation sound a good first approach to try. Nonetheless, in order to keep this a good question, which specific details of the problem you recommend me to provide? – manifold Jul 14 '21 at 00:55
  • @GregMartin Also, one thing that comes to my mind, is around which point should I try this first order approximation? – manifold Jul 14 '21 at 00:57
  • 1
    How large is $n$ ? Do you think that the solution is unique ? – Claude Leibovici Jul 14 '21 at 03:19
  • @ClaudeLeibovici I can reduce $n$ to a problem of around $20\times20$ (very cumbersome) equations by substituting variables with closed form expressions. Uniqueness of solution I would not be very sure, but I'd dare to say that aside from possible trivial solutions, there's an unique solution, or at least a real-valued solution that satisfies some restrictions (like being positive for some variables or similar) would work fine. I'll put those details in the answer. – manifold Jul 14 '21 at 03:33
  • 1
    As a general rule, your options are going to depend on $n$. If $n$ is roughly in the range of 1-100, then you can often find a good solution using simulated annealing, which is much less sensitive to the initial condition than other common methods (e.g. gradient descent). From your comment, it sounds like $n$ is either 20 or 400. If it is 20, you can probably find a solution in a reasonable amount of time if one exists. If 400, this is less likely. Of course, in both cases, this depends on the problem and the actual density of solutions. – sasquires Jul 14 '21 at 03:48
  • 1
    (If this wasn't clear, then for simulated annealing you would use the following energy or cost function. If your system is $f_i(x)=0$ for $i=1,\dots,n$, use $|f(x)|$ for some vector norm.) – sasquires Jul 14 '21 at 03:53
  • 1
    For even smaller $n$ (e.g. $n<10$) or functions that are analytically easy to handle, then I would recommend looking at nonlinear programming techniques. For larger $n$ (e.g. $n>100$), this is typically a very hard problem. Often people just pick random initial conditions and use gradient descent to try to find a solution and repeat until one is found. – sasquires Jul 14 '21 at 03:57
  • @sasquires First thanks for the edit recommendation, I'm still an English learner. Thanks for the method recommendation, I'll check it and try to implement it on Matlab. Also to be more specific (I'll put that on the original question) my full system would be no more than 90 equations, those of which around only (maximum) 25 would be involved in non-linearities, the rest either are dependant on parameters only or dependant on those 25 non-linear but with close form. – manifold Jul 14 '21 at 03:59
  • I didn't try it for your system, but before I've successfully solved a bunch of tough systems with https://en.wikipedia.org/wiki/Numerical_algebraic_geometry#Homotopy_continuation method. Briefly, you start with a system which has obvious solution and gradually move toward the original system. If you meet no critical points then you're done. Otherwise try to choose another path. – uranix Jul 14 '21 at 10:01

1 Answers1

2

The problem is a bit too vague.

You have $n$ equations in the form $$f_i(x_1,x_2,\cdots,x_n)=0$$ Assuming that in a known domain (may be defined by equality/inequality constraints on the variable) the solution is unique, I should try the minimization of the norm $$\Phi(x_1,x_2,\cdots,x_n)=\sum_{i=1}^n \Big[f_i(x_1,x_2,\cdots,x_n)\Big]^2$$ subject to the possible constraints.

If $\Phi_{\text{min}} \sim 0$, you are done.

If the equations are really nasty, do not try to provide the analytical Jacobian/Hessian. Let the optimizer working with numerical derivatives.

This is a very standard procedure.

At least, make a try and, please, let me know.

Claude Leibovici
  • 260,315
  • Thanks! I'll try to try that as soon as possible, also will try making some changes to the original answer in order to be more explicit in what I'm looking for. – manifold Jul 14 '21 at 04:03
  • Hi. Just was able to try quadratic minimization approach, but with no luck, local minima is not zero. Worth mentioning that I used ones for initial guess, the good thing is that at least I got a real solution using fmincon(), since as I was using fsolve() despite not having complex valued function at the initial guess, Matlab just stated introducing non-real numbers to the solution tryings. – manifold Jul 19 '21 at 18:07