How to solve these complex equations for plot?

Question

I am trying to Plot this code but it is not showing any plot. I want a bistable type plot between "NL" and "P0" and there are six equations and six unknowns. The code is given below,

delta = 4;
g0 = 1;
Pb = 2;
KL = 1;
kb = 0.2;
ome1 = 1;
gma = 0.02;
nth = 2;
sol[P0_?NumericQ] := 
  NSolve[{I*P0*(Conjugate[a1] - a1) + 
      I*Pb*(Lf - Conjugate[Lf] + Conjugate[Bf] - Bf) - KL*NL - 
      kb*(Lf*a1 + Conjugate[Lf]*Conjugate[a1] + Bf*Conjugate[a1] + 
         Conjugate[Bf]*
          a1 + (NL - 2*Abs[a1]^2)*(b1 + Conjugate[b1])) == 0, 
    I*g0*((NL - 2*Abs[a1]^2)*(Conjugate[b1] - b1) + 
         Conjugate[Lf]*Conjugate[a1] - Lf*a1 + Conjugate[Bf]*a1 - 
         Bf*Conjugate[a1]) + 
      I*Pb*(Conjugate[Bf] - Lf + Conjugate[Lf] - Bf) - gma*Nb + 
      gma*nth == 
     0, -I*delta*Lf - I*ome1*Lf + 
      I*g0*((2*NL - Abs[a1]^2 - Nb + 2*Abs[b1]^2 + b1^2)*
          Conjugate[a1] - (Conjugate[Bf] + 2*Lf)*b1 - 
         Lf*Conjugate[b1]) - I*P0*b1 == 0, 
    I*delta*Bf - I*ome1*Bf + 
      I*g0*((2*NL - Abs[a1]^2 + Nb - 2*Abs[b1]^2)*
          a1 + (2*Bf - a1*b1 + Conjugate[Lf])*b1 + Bf*Conjugate[b1]) +
       I*P0*b1 + I*Pb*(NL + Nb + a1^2 + b1^2 + 1) - (KL + gma)*Bf/2 - 
      kb/2*((2*Bf - a1*b1 + Conjugate[Lf])*b1 + (Nb - 2*Abs[b1]^2)*
          a1 + Bf*Conjugate[b1]) == 0, 
    I*delta*a1 + I*delta*a1 + I*g0*(Conjugate[Lf] + Bf) + I*P0 + 
      I*Pb*(b1 + Conjugate[b1]) - KL*a1/2 - 
      kb*(Conjugate[Lf] + Bf)/2 == 0, 
    I*ome1*b1 + I*g0*NL + I*Pb*(a1 + Conjugate[a1]) - gma*b1/2 == 
     0}, {a1, b1, NL, Lf, Bf, Nb}];
Plot[{NL /. sol[P0]}, {P0, 0, 30}]
This code is showing error like this;
NSolve::ivar: 0.07` is not a valid variable. >>
ReplaceAll::reps: {NSolve[{(0. +0. I)-NL+2 I (Times[<<2>>]+Lf+Conjugate[<<1>>]+Times[<<2>>])-0.2 (Times[<<2>>]+Times[<<2>>]+Times[<<2>>]+Times[<<2>>]+Times[<<2>>])==0,0.04 -0.02 b1 Conjugate[b1]+I (Times[<<2>>]+Times[<<2>>]+Times[<<2>>]+Times[<<2>>]+Times[<<2>>])+2 I (Times[<<2>>]+Times[<<2>>]+Conjugate[<<1>>]+Conjugate[<<1>>])==0,(0. -0.000612857 I) b1-5 I Lf+I (Times[<<3>>]+Times[<<2>>]+Times[<<3>>])==0,<<6>>+I (Times[<<2>>]+<<1>>+Times[<<2>>])==0,(-0.035+0.560613 I)+2 I (b1+Conjugate[<<1>>])-(0.1 -1. I) (Bf+Conjugate[<<1>>])==0,(0. +0.28 I)-(0.01 -1. I) b1+I NL==0},{0.07,b1,NL,Lf,Bf}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>
If anyone can resolve this will be appreciated.

Answer 1

We can reduce system to the real form using substitutions eqs=ComplexExpand[sys/.{a1 -> a1x + I a1y, b1 -> b1x + I b1y, NL -> NLx + I NLy, Lf -> Lfx + I Lfy, Bf -> Bfx + I Bfy, Nb -> Nbx + I Nby}]. Then we can divide eqs into two part eq1 = eqs /. I -> 0; eq2 = eqs /. I -> 1;. Let us join equations Join[eq1, eq2] /. P0 -> p and define function (all parameters are numeric using data delta = 4;g0 = 1;Pb = 2;KL = 1;kb = 0.2;ome1 = 1;gma = 0.02;nth = 2;)

eqf[p_] := {0.8` a1x^2 b1x + 0.8` a1y^2 b1x - 0.4` a1x Bfx + 4 Bfy - 
    0.4` a1y Bfy - 0.4` a1x Lfx - 4 Lfy + 0.4` a1y Lfy - NLx - 
    0.4` b1x NLx + 2 a1y p == 0, 
  0.04` - 4 a1x^2 b1y - 4 a1y^2 b1y - 2 a1y Bfx + 4 Bfy + 2 a1x Bfy + 
    2 a1y Lfx + 4 Lfy + 2 a1x Lfy - 0.02` Nbx + 2 b1y NLx == 
   0, -a1x^2 a1y - a1y^3 + 3 a1y b1x^2 - 2 a1x b1x b1y + a1y b1y^2 + 
    b1y Bfx - b1x Bfy + b1y Lfx + 5 Lfy + 3 b1x Lfy - a1y Nbx + 
    a1x Nby + 2 a1y NLx - 2 a1x NLy + b1y p == 
   0, -4 a1x a1y + a1x^2 a1y + a1y^3 + 
    0.30000000000000004` a1x b1x^2 + 3 a1y b1x^2 - 4 b1x b1y + 
    2 a1x b1x b1y - 0.2` a1y b1x b1y + 0.1` a1x b1y^2 + a1y b1y^2 - 
    0.51` Bfx - 0.30000000000000004` b1x Bfx - b1y Bfx - 3.` Bfy - 
    3 b1x Bfy + 0.1` b1y Bfy - 0.1` b1x Lfx - b1y Lfx + b1x Lfy - 
    0.1` b1y Lfy - 0.1` a1x Nbx - a1y Nbx - 2 Nby - a1x Nby + 
    0.1` a1y Nby - 2 a1y NLx - 2 NLy - 2 a1x NLy - b1y p == 
   0, -(a1x/2) - 8 a1y - 0.1` Bfx - 1.` Bfy - 0.1` Lfx + 1.` Lfy == 
   0, -0.01` b1x - 1.` b1y - NLy == 0, 
  0.8` a1x^2 b1x + 0.8` a1y^2 b1x - 0.4` a1x Bfx + 4 Bfy - 
    0.4` a1y Bfy - 0.4` a1x Lfx - 4 Lfy + 0.4` a1y Lfy - NLx - 
    0.4` b1x NLx - NLy - 0.4` b1x NLy + 2 a1y p == 0, 
  0.04` - 4 a1x^2 b1y - 4 a1y^2 b1y - 2 a1y Bfx + 4 Bfy + 2 a1x Bfy + 
    2 a1y Lfx + 4 Lfy + 2 a1x Lfy - 0.02` Nbx - 0.02` Nby + 
    2 b1y NLx + 2 b1y NLy == 
   0, -a1x^3 - a1x^2 a1y - a1x a1y^2 - a1y^3 + 3 a1x b1x^2 + 
    3 a1y b1x^2 - 2 a1x b1x b1y + 2 a1y b1x b1y + a1x b1y^2 + 
    a1y b1y^2 - b1x Bfx + b1y Bfx - b1x Bfy - b1y Bfy - 5 Lfx - 
    3 b1x Lfx + b1y Lfx + 5 Lfy + 3 b1x Lfy + b1y Lfy - a1x Nbx - 
    a1y Nbx + a1x Nby - a1y Nby + 2 a1x NLx + 2 a1y NLx - 2 a1x NLy + 
    2 a1y NLy - b1x p + b1y p == 0, 
  2 + 2 a1x^2 - a1x^3 - 4 a1x a1y + a1x^2 a1y - 2 a1y^2 - a1x a1y^2 + 
    a1y^3 + 2 b1x^2 - 2.7` a1x b1x^2 + 3.3` a1y b1x^2 - 4 b1x b1y + 
    2.2` a1x b1x b1y + 1.8` a1y b1x b1y - 2 b1y^2 - 0.9` a1x b1y^2 + 
    1.1` a1y b1y^2 + 2.49` Bfx + 2.7` b1x Bfx - 1.1` b1y Bfx - 
    3.51` Bfy - 3.3` b1x Bfy - 0.9` b1y Bfy + 0.9` b1x Lfx - 
    1.1` b1y Lfx + 1.1` b1x Lfy + 0.9` b1y Lfy + 2 Nbx + 
    0.9` a1x Nbx - 1.1` a1y Nbx - 2 Nby - 1.1` a1x Nby - 
    0.9` a1y Nby + 2 NLx + 2 a1x NLx - 2 a1y NLx - 2 NLy - 
    2 a1x NLy - 2 a1y NLy + b1x p - b1y p == 
   0, (15 a1x)/2 - (17 a1y)/2 + 4 b1x + 0.9` Bfx - 1.1` Bfy + 
    0.9` Lfx + 1.1` Lfy + p == 0, 
  4 a1x + 0.99` b1x - 1.01` b1y + NLx - NLy == 0};

Now we use it as

vars = {a1x, a1y, b1x, b1y, NLx, NLy, Lfx, Lfy, Bfx, Bfy, Nbx, Nby};
sols[p_] := 
 FindRoot[eqf[p], Table[{vars[[i]], 1/10}, {i, Length[vars]}]] // 
  Quiet 

Finally plot

ListLinePlot[Table[{p, NLx /. sols[p]}, {p, 3, 30, .5}]]
{ListLinePlot[Table[{p, NLx /. sols[p]}, {p, 0, 3, .05}]], 
 ListLinePlot[Table[{p, NLy /. sols[p]}, {p, 0, 3, .05}]]}

Answer 2

Even after making corrections mentioned in the comments, NSolve is too slow for me and it times out trying to evaluate a single point. Since all equations are of the form $f_i(a_1, b_1, N_L, L_f, B_f) = 0$, we could instead recast the problem as a minimization. Since $f_i$ may be complex valued, we want to minimize the complex norm squared $\sum_i|f_i(a_1, b_1, N_L, L_f, B_f)|^2$ which should achieve a value close to zero near the solution if it exists.

delta = 4;
g0 = 1;
Pb = 2;
KL = 1;
kb = 0.2;
ome1 = 1;
gma = 0.02;
nth = 2;
(* original equations with '== 0' removed. )
equations = {
  IP0(Conjugate[a1] - a1) + IPb(Lf - Conjugate[Lf] + Conjugate[Bf] - Bf) - KLNL - kb(Lfa1 + Conjugate[Lf]Conjugate[a1] + BfConjugate[a1] + Conjugate[Bf]a1 + (NL - 2Abs[a1]^2)(b1 + Conjugate[b1])),
  Ig0((NL - 2Abs[a1]^2)(Conjugate[b1] - b1) + Conjugate[Lf]Conjugate[a1] - Lfa1 + Conjugate[Bf]a1 - BfConjugate[a1]) + IPb(Conjugate[Bf] - Lf + Conjugate[Lf] - Bf) - gmaNb +  gmanth,
 -IdeltaLf - Iome1Lf + Ig0((2NL - Abs[a1]^2 - Nb + 2Abs[b1]^2 + b1^2)Conjugate[a1] - (Conjugate[Bf] + 2Lf)b1 - LfConjugate[b1]) - IP0b1,
  IdeltaBf - Iome1Bf + Ig0((2NL - Abs[a1]^2 + Nb - 2Abs[b1]^2)a1 + (2Bf - a1b1 + Conjugate[Lf])b1 + BfConjugate[b1]) +IP0b1 + IPb(NL + Nb + a1^2 + b1^2 + 1) - (KL + gma)Bf/2 - kb/2((2Bf - a1b1 + Conjugate[Lf])b1 + (Nb - 2Abs[b1]^2)a1 + BfConjugate[b1]),
  Ideltaa1 + Ideltaa1 + Ig0(Conjugate[Lf] + Bf) + IP0 + IPb(b1 + Conjugate[b1]) - KLa1/2 - kb(Conjugate[Lf] + Bf)/2,
  Iome1b1 + Ig0NL + IPb(a1 + Conjugate[a1]) - gmab1/2
}//ComplexExpand//FullSimplify;
minm = Total[Abs[#]^2 & /@ equations];
sol[p_?NumericQ] := NMinimize[minm /. P0 -> p, {a1, b1, NL, Lf, Bf, Nb}]

Testing it out on a value gives:

sol[5.0]
(* {0.630615, {a1 -> 0.0718601, b1 -> -1.0705, NL -> 0.587601, Lf -> 1.48774, Bf -> -2.72366, Nb -> 0.511481}} *)

The error of 0.630615 is quite large, and while close to zero, this is not that close. Unfortunately, to get closer to the solution we may need NMinimize to work over the Complexes to pick complex values and this hasn't been implemented in Mathematica yet.

We can plot the error for all P0:

ListLinePlot[{#, First[sol[#]]} & /@ Range[0, 30, 0.25], AxesLabel -> {"P0", "MSE"}]

There may be solutions near 0, 11, 29, and 30. The best P0 was somewhere around 30.0 and if we evaluate the equations at the minimization solution we get a bit closer to zero for all of them:

With[{p = 30},
 equations /. Flatten[{Last[sol[p]], P0 -> p}]
]
(* result: {-0.00474041 + 0. I, -0.0430146, 
 0. + 0.00776613 I, -0.131341 - 0.00697846 I, -0.65231 - 0.0632559 I, 
 0.0304914 + 0.0276885 I} *)

So, all in all it's not good news. NMinimize is probably getting trapped in local minima, cannot utilize complex solutions, and there may be problems in the equations shown in Hugh's answer. It's not possible to draw an accurate NL,P0 plot until these problems are resolved, unless your equations allow some slack. With that said, here's what it looks like right now:

ListLinePlot[{#, NL /. Last[sol[#]]} & /@ Range[0, 30, 0.25], AxesLabel -> {"P0", "NL"}]

Answer 3

I suggest you don't rush to plot. Look at your equations in detail. Starting again:

delta = 4;
ome1 = 1;
Pb = P0*g0/(2*ome1);
g0 = 1;
KL = 1;
kb = 0.2;
gma = 0.02;
nth = 2;
eqns = {IP0(Conjugate[a1] - a1) + 
     IPb(Lf - Conjugate[Lf] + Conjugate[Bf] - Bf) - KLNL - 
     kb(Lfa1 + Conjugate[Lf]Conjugate[a1] + BfConjugate[a1] + 
        Conjugate[Bf]a1 + (NL - 2Abs[a1]^2)(b1 + Conjugate[b1])) ==
     0, Ig0((NL - 2Abs[a1]^2)(Conjugate[b1] - b1) + 
        Conjugate[Lf]Conjugate[a1] - Lfa1 + Conjugate[Bf]a1 - 
        BfConjugate[a1]) + 
     IPb(Conjugate[Bf] - Lf + Conjugate[Lf] - Bf) - 
     gmaConjugate[b1]b1 + gmanth == 
    0, -IdeltaLf - Iome1Lf + 
     Ig0((2NL - Abs[a1]^2 - Conjugate[b1]b1 + 2Abs[b1]^2 + b1^2)*
         Conjugate[a1] - (Conjugate[Bf] + 2Lf)b1 - 
        LfConjugate[b1]) - IP0b1 == 0, 
   IdeltaBf - Iome1Bf + 
     Ig0((2NL - Abs[a1]^2 + Conjugate[b1]b1 - 2Abs[b1]^2)*
         a1 + (2Bf - a1b1 + Conjugate[Lf])b1 + BfConjugate[b1]) + 
     IP0b1 + 
     IPb(NL + Conjugate[b1]b1 + a1^2 + b1^2 + 1) - (KL + gma)
      Bf/2 - 
     kb/2((2Bf - a1b1 + Conjugate[Lf])
         b1 + (Conjugate[b1]b1 - 2Abs[b1]^2)a1 + 
        BfConjugate[b1]) == 0, 
   Ideltaa1 + Ideltaa1 + Ig0(Conjugate[Lf] + Bf) + IP0 + 
     IPb(b1 + Conjugate[b1]) - KLa1/2 - 
     kb(Conjugate[Lf] + Bf)/2 == 0, 
   Iome1b1 + Ig0NL + IPb(a1 + Conjugate[a1]) - gmab1/2 == 0};
eqns1 = ComplexExpand[eqns] // Simplify
(* 
{1. a1^2 b1 == 0.5 a1 Bf + 0.5 a1 Lf + 1.25 NL + 0.5 b1 NL, 
 1. b1^2 == 2., a1^3 + 5 Lf + b1 (Bf + 3 Lf + P0) == 2 a1 (b1^2 + NL),
  0. + 0.2 a1 b1^2 - 0.51 Bf - 0.3 b1 Bf - 0.1 b1 Lf + 
   I (2. + 2 a1^2 - a1^3 + 4 b1^2 + 3. Bf - 2 a1 (b1^2 - NL) + 2 NL + 
      b1 (3 Bf + Lf + P0)) == 
  0, (0. + 1. I) a1 + (0.0311284 + 0.498054 I) b1 - (0.00466926 - 
      0.125292 I) Bf - (0.00466926 - 0.125292 I) Lf + (0.0077821 + 
      0.124514 I) P0 == 0, 0. + I (0. + 4 a1 + 1. b1 + NL) == 0.01 b1}


*)

I have put in ComplexExpand to get rid of all the Conjugates . This assumes that all your variables are real. Is this correct? We now have a much simpler set of equations.

We can now useReduce

sol = Reduce[eqns1, {a1, b1, NL, Lf, Bf}]
(*
((P0 == -63.7447 - 12.4296 I || P0 == -20.0536 + 84.515 I || 
     P0 == -1.84633 + 0.27489 I || P0 == 30.3551 - 105.254 I) && 
   a1 == (-1.1398110^-289 - 
       2.4956110^-289 I) ((6.3941710^287 - 
         1.6252410^288 I) + (8.2397410^286 - 
          2.3493910^287 I) P0 - (1.825410^284 - 
          5.208410^284 I) P0^2 + (2.3979210^282 + 
          9.1318210^282 I) P0^3) && b1 == 1.41421 && 
   NL == (5.2181110^-299 + 
       1.142510^-298 I) ((8.0670110^296 - 
         4.0052410^297 I) + (7.1993410^296 - 
          2.0527410^297 I) P0 - (1.5949210^294 - 
          4.5507610^294 I) P0^2 + (2.0951410^292 + 
          7.9787810^292 I) P0^3) && 
   Lf == (2.1597610^-299 - 
       3.9573610^-299 I) ((-3.9400410^297 - 
         4.7053710^297 I) - (3.8057210^297 + 
          6.3505610^297 I) P0 - (2.4874410^295 + 
          3.2516310^295 I) P0^2 + (1.5984810^294 + 
          1.7269610^294 I) P0^3) && 
   Bf == (8.0956810^-284 + 
       2.754910^-283 I) ((-9.6455510^281 + 
         5.2469210^282 I) - (6.1316910^280 - 
          4.0402810^281 I) P0 + (1.8984110^278 - 
          2.2175910^279 I) P0^2 + (1.6165310^277 + 
          4.3439810^278 I) P0^3)) || ((P0 == -2.40298 - 1.07568 I || 
     P0 == 3.98485 + 0.345156 I || P0 == 13.0124 - 13.6786 I || 
     P0 == 32.4428 - 0.156766 I) && 
   a1 == (-3.614910^-286 - 
       6.3902910^-287 I) ((-1.8431510^284 + 
         8.5912410^284 I) + (8.3755710^284 - 
          2.719210^284 I) P0 - (6.2339410^283 + 
          3.2811410^282 I) P0^2 + (1.6032810^282 + 
          2.4347510^281 I) P0^3) && b1 == -1.41421 && 
   NL == (8.2745910^-295 + 
       1.4627510^-295 I) ((1.3381610^294 + 
         1.2248910^294 I) + (1.463610^294 - 
          4.7517210^293 I) P0 - (1.0893610^293 + 
          5.7336910^291 I) P0^2 + (2.8016810^291 + 
          4.2546610^290 I) P0^3) && 
   Lf == (6.8118410^-294 + 
       9.7373510^-294 I) ((-5.9525710^291 + 
         2.8125210^291 I) + (3.9971710^292 - 
          7.2149410^292 I) P0 - (4.7136110^291 - 
          2.6823710^291 I) P0^2 + (8.0871510^289 + 
          1.8031610^289 I) P0^3) && 
   Bf == (6.5876610^-278 - 
       3.7265410^-277 I) ((-2.7313610^276 + 
         1.3333710^277 I) + (9.1611110^275 + 
          1.3768110^276 I) P0 - (3.7696210^274 + 
          3.1869910^275 I) P0^2 + (7.7561310^272 + 
          1.1070610^274 I) P0^3))
*)

This suggests that there are only four values for P0. The numerical coefficients are huge which is worrying. You can look at each of the solutions for your variables in terms of the four values of P0. The answers are complex which does go against the use of ComplexExpand. However, you need to look further at your equations.