NSolve cannot solve system of equations in Mathematica

Question

I have the following two equations:

Ieq1=-1 - (0.763546 I1 (1/(-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2) - (
18 (1 - 9/I1 + 10/I2))/(
I1 (-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2)^2)) X^1.3)/((1.36- 
I1) (I1/(-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2))^2.3) + (
0.587343 X^1.3)/((1.36- 
I1) (I1/(-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2))^1.3) + (
0.587343 I1 X^1.3)/((1.36- 
I1)^2 (I1/(-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2))^1.3);

Ieq2=1 + (0.763546 I1^2 ((20 (-8 + 10/I2))/I2^2 - (20 (1 - 9/I1 + 10/I2))/
I2^2) (-8 + (1 - 9/I1 + 10/I2)^2) X^1.3)/((1.36- 
I1) (I1/(-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/
  I2)^2))^2.3 (-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2)^3) - (
 0.587343 I1 ((20 (-8 + 10/I2))/I2^2 - (20 (1 - 9/I1 + 10/I2))/
I2^2) (-8 + (1 - 9/I1 + 10/I2)^2) X^1.3)/((1.36- 
I1) (I1/(-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/
  I2)^2))^1.3 (-(-8 + 10/I2)^2 + (1 - 9/I1 + 10/I2)^2)^2);

and then I use the NSolve function:

 NSolve[{Ieq1 == 0, Ieq2 == 0}, {I1, I2}]

because I want to solve each one as a function of X. But unfortunately it just gets stuck while running. Do you have any suggestions to overcome this problem?

Answer 1

This will not give you the complete solution for your problem (which appears not quite trivial), but rather some startup to get used to a wider range of tools Mathematica gives you. If you don't understand what some function used in my code does, you can always move text cursor to its name and press F1 in Mathematica notebook interface to see documentation.

First, to ease work with your big functions, let's define them as functions instead of just symbols (this will be necessary for Plot3D, for instance):

Clear[Ieq1,Ieq2]
Ieq1[X_,I1_,I2_]=YourFirstFunction;
Ieq2[X_,I1_,I2_]=YourSecondFunction;

Now let's see how these two functions look with an example X=0.2:

With[{X = 0.2},
Plot3D[
{
    Log[Abs[Ieq1[X, I1, I2]]],
    Log[Abs[Ieq2[X, I1, I2]]]
}, {I1, -10, 10}, {I2, -10, 10},
 PlotRange -> {-10, 2}, PlotPoints -> 100, PlotStyle -> {Red, Green}]
]

We used Log[Abs[...]] here to see zeroes of each function. You can see that there's plenty of zeroes for each function.

But, let's see whether your functions have any simultaneous roots (I use 500 plot points here to see fine details of the function; if your machine is not too fast this may take several minutes to render and appear sluggish when trying to rotate):

With[{X = 0.2}, Plot3D[
{
    Log[Abs[Ieq1[X, I1, I2]] + Abs[Ieq2[X, I1, I2]]]
}, {I1, -10, 10}, {I2, -10, 10}, PlotRange -> {-10, 2}, 
PlotPoints -> 500, PlotStyle -> {Red, Green}, Exclusions -> None]
]

So, we have a candidate for a root. Let's try it with FindRoot with nearby starting values:

res=With[{X = 0.2}, Chop@FindRoot[
{
    Ieq1[X, I1, I2], Ieq2[X, I1, I2]
}, {I1, 0.5}, {I2, 0.5}, MaxIterations -> 1000, WorkingPrecision -> 100]
]

{I1 -> 1., I2 -> 0.354532}

You'll get a FindRoot::lstol warning, but you do get some values. Let's check whether they are good enough:

Ieq1[0.2, I1, I2] /. res
Ieq2[0.2, I1, I2] /. res

-8.93913*10^-8

-4.26765*10^-7

So, OK, we've found the roots for X=0.2. You'll have now to use Table to make final list of {X,{I1,I2} elements, but:

NOTE Using FindRoot you always have to check if the roots found are indeed the roots, especially if you get some warnings. This may make your experience with automatizing your computations with Table considerably harder. If you get wrong values, try perturbing starting values. Also, always try to visualize the functions you work with if it's possible (i.e. they have not too many arguments).