I am trying to solve the given equation in the code and want to plot 'N1' with 'P1'. But the code is not giving me any output and still showing running. What is the problem with it? If anyone can figure it out will be appreciated.
dela = -1.5;(*Subscript[\[CapitalDelta], L] ; Laser-Lower polariton \
detuning*)
g1 = 0.0003; (* Subscript[g, 0] ; vacuum optomechanical strength *)
del0 = -0.006; (* \[Delta] ; exciton-photon detuning*)
ome = 3.014; (* \[CapitalOmega] ; Rabi splitting *)
omem = 0.125; (* Subscript[\[CapitalOmega], m] ; mechanical \
resonator's frequency *)
kexc = 0.002; (* Subscript[\[Kappa], exc ; ]exciton decay rate*)
kk = .2; (* \[Kappa]; cavity decay rate *)
gma = 0.00001; (* \[CapitalGamma] ; phonon dcay rate *)
nth = 0;
kal = (kk + kexc)/2 - del0*(kk - kexc)/(2*ome);
g0 = g1*(1 - del0/ome)/2;
kb = 0.1;
eq1 = -I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*(Conjugate[a1] - a1) -
I*(P1*g1/(2*Sqrt[2]*ome))*(alab - Conjugate[alab] +
Conjugate[betab] - betab) + kal*N1 +
kb*(alab*a1 + Conjugate[alab]*Conjugate[a1] +
betab*Conjugate[a1] +
Conjugate[betab]*
a1 + (N1 - 2*(Abs[a1])^2)*(b1 + Conjugate[b1])) == 0;
eq2 = I*dela*alab + I*omem*alab -
I*g0*((2*N1 - Abs[a1]^2 - Np + 2*2*(Abs[b1])^2 + b1^2)*
Conjugate[a1] - (Conjugate[betab] + 2*alab)*b1 -
alab*Conjugate[b1]) + I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 - Np + (Conjugate[a1])^2 -
b1^2) + (kal + gma)*alab/2 +
kb/2*(Conjugate[betab]*b1 +
alab*Conjugate[b1] + (Np - 2*(Abs[b1])^2)*Conjugate[a1] +
2*alab*b1 - b1^2*Conjugate[a1]) == 0;
eq3 = -I*g0*((N1 - 2*(Abs[a1])^2)*(Conjugate[b1] - b1) +
Conjugate[alab]*Conjugate[a1] - alab*a1 + Conjugate[betab]*a1 -
betab*Conjugate[a1]) -
I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab] - alab +
Conjugate[alab] - betab) + gma*Np - gma*nth == 0;
eq4 = -I*dela*betab + I*omem*betab -
I*g0*((2*N1 - (Abs[a1])^2 + Np - 2*(Abs[b1])^2)*
a1 - (2*betab - a1*b1 + Conjugate[alab])*b1 +
betab*Conjugate[b1]) - I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 + Np + a1^2 + b1^2 + 1) + (kal +
gma)*betab/2 +
kb/2*((2*betab - a1*b1 + Conjugate[alab])*
b1 + (Np - 2*(Abs[b1])^2)*a1 + betab*Conjugate[b1]) == 0;
eq5 = -I*dela*a1 - I*g0*(Conjugate[alab] + betab) -
I*(P1/Sqrt[2]*(1 - del0/(2*ome))) -
I*(P1*g1/(2*Sqrt[2]*ome))*(b1 + Conjugate[b1]) + kal*a1/2 +
kb*(Conjugate[alab] + betab)/2 == 0;
eq6 = I*omem*b1 - I*g0*N1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(a1 + Conjugate[a1]) + gma*b1/2 == 0;
sol[P1_?NumericQ] := FindRoot[{eq1, eq2, eq3, eq4, eq5, eq6}, {N1, alab, Np, betab, a1,
b1}];
Plot[{Evaluate[N1] /. sol[P1]}, {P1, 0, 20}]
And it will really help if you write the ODE's each in separate equation, then use these later inside NDSolve instead of writing everything all in one command. No need to do that. Makes it hard for you and others to see things.
The same here. Better write each equation on its own (sayeq1
,eq2
etc..).then use these inside theSolve
commands to make it more clear. $\endgroup$P1
. Are there any known constraints on{N1, alab, Np, betab, a1, b1}
? $\endgroup$sol
without including its argument.N1 /. sol[P1]
$\endgroup$...initialization...;tbl=Table[{P1,N1}/.NMinimize[Abs[..]+Abs[..]+Abs[..]+Abs[..]+Abs[..]+Abs[..],{N1,alab,Np,betab,a1,b1}][[2]],{P1,0,20,1}];ListPlot[tbl,Joined->True]
pretty quickly shows a plot of P1 versus N1 for the range of P1 with no warnings or errors when I substitute in the left hand side of your six equations inside each of thoseAbs[..]
and make no other changes. Adjust the step size of P1 as desired. Looking at the full result from the NMinimize can show you how well it approximated zero. $\endgroup$