How to solve a large system of non-linear equations

I would like to solve Saha equations on Fe, which have 27 unknowns. Saha equations basically give you the ratios of the unknowns, but the unknowns are subject to constraints.

What function shall I use to solve this large system of equations?

Because each equation is only different by a parameter value, I do not want to type them one by one. So: Is there any easier way in Mathematica?

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Feb 15 '15 at 17:45
• Please show what you have tried so far, preferably with an example with many fewer equations. – bbgodfrey Feb 15 '15 at 17:46
• I imagine that you could use NSolve or LinearSolve, depending on the details of your equations. Search Mathematica.StackExchange for ideas and examples. – bbgodfrey Feb 16 '15 at 2:28
• @bbgodfrey I am wondering is there any easier way to key in equations to NSolve, say if I have 27 variables,ni i from 0 to 26, they all obey similar equations, can I just type in the general form of the equation and the constraint which is the sum of ni being some constant and ask NSolve to solve ni for me? – Shadumu Feb 19 '15 at 17:09
• Please add at least a sample of your equations to the question. It is very difficult to offer advice otherwise. – bbgodfrey Feb 19 '15 at 17:13

The generic equations to be solved are in the OP's Comments. To answer the second question first, the equations can be written compactly as

imax = 26;
eqs = Table[n0 n[i + 1] == C[i] n[i], {i, imax - 1}];
var = Table[n[i], {i, imax}];
eqnrm = Sum[n[i], {i, imax}] == 1;
eqsum = Sum[i n[i], {i, imax}] - n0 == 0;


The first question is answered as follows. In the absence of the last of these equations, the rest form a linear system which Mathematica can solve rapidly

slin = Solve[Join[eqs, {eqnrm}], var]


Not surprisingly, all elements of var are found to have a common denominator,

Denominator[slin[[1, 1, 2]]];


which can be used to simplify the remaining equation to a polynomial of order imax in n0 before it is solved.

eqn0 = Simplify[First[First[eqsum] /. slin] Denominator[slin[[1, 1, 2]]]] == 0;
Solve[eqn0, n0]


With n0 determined, the var follow immediately by back substitution. The calculation takes a few minutes for imax = 14 and progressively longer for larger values of imax.

In answer to an additional question posed by the OP in a comment, back substitution is accomplished as follows.

The process above yields imax solutions for n0. Presumably, physical considerations will guide the selection of which should be selected. Suppose it is solution 2. Then, with sn0 the list of solutions from the final Solve above,

j = 2; sn0j = sn0[[j]];
Table[Simplify[slin[[1, i, 2]] /. sn0j], {i, imax}]


accomplishes the back substitution.

Note, however, that this set of solutions rapidly becomes unmanageable in symbolic form for large imax. Therefore, numerical values of C probably should be substituted into the equations from the outset, and the equations solved using NSolve rather than Solve. If the C are functions of time, it may be better to determine the solutions for specific values of time and then, if desired, generate 'InterpolatingFunction`s rather than to attempt to obtain solutions symbolically as a function of time. Only experimentation can answer which approach is superior.

• how to do back substitution? – Shadumu Feb 20 '15 at 12:39
• Please provide C (i), so that I can try to reproduce what you are doing. – bbgodfrey Feb 20 '15 at 17:05
• Now, I have 26 equations only, I removed one constraint. – Shadumu Feb 20 '15 at 17:40
• L = {7.9024, 16.1878, 30.652, 54.8, 75.0, 99.1, 124.98, 151.06, 233.6, 262.1, 290.2, 330.8, 361.0, 392.2, 457, 489.256, 1266, 1358, 1456, 1582, 1689, 1799, 1950, 2023, 8828, 9277.69} T = 10^7 lambda = (6.63*10^-34)/(2*3.1415926*9.10938291*10^-31*1.3806488*10^-23*T)^0.5 K = 2/lambda^3 ne = 10^11 imax = 26; eqs = Table[ ne n[i + 1] == K*n[i] Exp[(L[[i]] - L[[i + 1]])*1.60217657/(T*1.3806488*10^-4)] , {i, imax - 1}]; var = Table[n[i], {i, imax}]; eqnrm = Sum[n[i], {i, imax}] == 10^6; slin = Solve[Join[eqs, {eqnrm}], var]; – Shadumu Feb 20 '15 at 17:41
• I got error msg : RowReduce::luc: Result for RowReduce of badly conditioned matrix .... may contain significant numerical errors. >> – Shadumu Feb 20 '15 at 17:44