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?

  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Feb 15, 2015 at 17:45
  • $\begingroup$ Please show what you have tried so far, preferably with an example with many fewer equations. $\endgroup$
    – bbgodfrey
    Commented Feb 15, 2015 at 17:46
  • $\begingroup$ I imagine that you could use NSolve or LinearSolve, depending on the details of your equations. Search Mathematica.StackExchange for ideas and examples. $\endgroup$
    – bbgodfrey
    Commented Feb 16, 2015 at 2:28
  • $\begingroup$ @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? $\endgroup$
    – Shadumu
    Commented Feb 19, 2015 at 17:09
  • $\begingroup$ Please add at least a sample of your equations to the question. It is very difficult to offer advice otherwise. $\endgroup$
    – bbgodfrey
    Commented Feb 19, 2015 at 17:13

1 Answer 1


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.

  • $\begingroup$ how to do back substitution? $\endgroup$
    – Shadumu
    Commented Feb 20, 2015 at 12:39
  • $\begingroup$ Please provide C (i), so that I can try to reproduce what you are doing. $\endgroup$
    – bbgodfrey
    Commented Feb 20, 2015 at 17:05
  • $\begingroup$ Now, I have 26 equations only, I removed one constraint. $\endgroup$
    – Shadumu
    Commented Feb 20, 2015 at 17:40
  • $\begingroup$ 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]; $\endgroup$
    – Shadumu
    Commented Feb 20, 2015 at 17:41
  • $\begingroup$ I got error msg : RowReduce::luc: Result for RowReduce of badly conditioned matrix .... may contain significant numerical errors. >> $\endgroup$
    – Shadumu
    Commented Feb 20, 2015 at 17:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.