Problems with solve command.

Hello,

I am trying to solve the following code for x. According to the attached plot, x should be ~1.3

however, using the solve command gives me this error:

Warning: The solutions are parameterized by the symbols: k, z1.
 To include parameters and conditions in the solution, specify the 'ReturnConditions' option. 
> In solve>warnIfParams (line 500)
  In solve (line 356)
  In HW5 (line 27) 
Warning: The solutions are valid under the following conditions: exp(log(z1) + k*pi*2i) ~=
-661055968790248598951915308032771039828404682964281219284648795274405791236311345825189210439715284847591212025023358304256/2969614242875447
& in(k, 'integer') & (z1 == root(z^3 +
(661055968790248598951915308032771039828404682964281219284648795274405791236311345825189210439715284847591212025023358304256*z^2)/2969614242875447
-
(7389015366435323967908908022371910771811603653921098701335406985236788968973433234508946305610994222893417207445031565100274626525888048356861860564629277155285935890388881077320961676529455407039532535682776566997321416812576672328618134497918976*z)/8552739147866811329799841636241
-
5572448313541411075572754442691595320934361841848741683192517222247028915535681815974691113703119944719040419346258432307987094672556923829465808318357201252183187273097626875632198907520424901383961794447304103842444201755391556919461542396321333248/8552739147866811329799841636241,
z, 1) | z1 == root(z^3 +
(661055968790248598951915308032771039828404682964281219284648795274405791236311345825189210439715284847591212025023358304256*z^2)/2969614242875447
-
(7389015366435323967908908022371910771811603653921098701335406985236788968973433234508946305610994222893417207445031565100274626525888048356861860564629277155285935890388881077320961676529455407039532535682776566997321416812576672328618134497918976*z)/8552739147866811329799841636241
-
5572448313541411075572754442691595320934361841848741683192517222247028915535681815974691113703119944719040419346258432307987094672556923829465808318357201252183187273097626875632198907520424901383961794447304103842444201755391556919461542396321333248/8552739147866811329799841636241,
z, 2) | z1 == root(z^3 +
(661055968790248598951915308032771039828404682964281219284648795274405791236311345825189210439715284847591212025023358304256*z^2)/2969614242875447
-
(7389015366435323967908908022371910771811603653921098701335406985236788968973433234508946305610994222893417207445031565100274626525888048356861860564629277155285935890388881077320961676529455407039532535682776566997321416812576672328618134497918976*z)/8552739147866811329799841636241
-
5572448313541411075572754442691595320934361841848741683192517222247028915535681815974691113703119944719040419346258432307987094672556923829465808318357201252183187273097626875632198907520424901383961794447304103842444201755391556919461542396321333248/8552739147866811329799841636241,
z, 3)).
 To include parameters and conditions in the solution, specify the 'ReturnConditions' option. 
> In solve>warnIfParams (line 507)
  In solve (line 356)
  In HW5 (line 27)

-----------------------------------------------------------------------------------------------------------------

Here is the code:

syms x

Econv=1.60218*10^-19; %%J/eV

Eg=1.11             %%Bandgap energy in eV

k=(1.38064852*10^-23)/Econv;   %%Bolzmann Constant in eV/K

m0=9.109*10^-31;     %%Mass of electron

mn=1.1*m0;           %%Mass of e carrier

mp=0.58*m0;           %%Mass of e hole
Eion=0.045;
hbar=1.054571800*10^-34;   %Planck's constant in Js
T=50;                  %Temperature range in K

Nd=((10^15)*(100)^3);      %%#Donors/m^3

Nc=2.*(mn*Econv*k.*T./(2*pi*hbar^2)).^(3/2);

Nv=2.*(mp*Econv*k.*T./(2*pi*hbar^2)).^(3/2);

Eiv=Eg/2+(3/4)*k.*T.*log(mp/mn);

ni=((Nc.*Nv).^(1/2)).*exp(-Eg./(2*k.*T));

p=ni.*exp(-x./(k.*T)).*(exp(Eiv./(k.*T)));

Ndion=Nd./(1+exp(x./(k.*T)).*exp((Eion-Eg)./(k.*T)));

LHS=p.*(p+Ndion);

Efv=solve(LHS==(ni.^2),x);

Answer 1

[John D'Errico]

Looks like you defined T in the comments.

Anyway, time to learn what vpa does for you.

vpa(Efv)
ans =
1.0706432917519708822549379645966 - 4.6193018557455934635744555230166e-41i

You probably need to get better at reading plots too. Your guess of 1.3 is not that close. :)

[mattyice]

right. I guess the fact that I was getting an imaginary component was throwing me off. Any idea why that's there? There shouldnt be any imaginary numbers here

[Walter Roberson]

double(Efv) does not show any imaginary component.

Your form is a cubic; the general solution to a cubic involves sqrt(-1) that might or might not cancel out. The roots of the cubic are in the order of 10^108, so you have a lot of potential for round-off error to affect the result such that the sqrt(-1) is not exactly balanced. Indeed, until you use at least 111 digits of computation, you do not even get the right integer portion of the third root (which is about -754); using only the default 32 digits the third root comes out about -10^77.

[John D'Errico]

The analytical solution to a polynomial uses complex arithmetic, because, well, there are often complex roots. A numerical solver won't generate a complex root, but then it won't be analytical.

The fact that the imaginary term had an attached power of 10^-41 should be a good hint that it really is not there. Just an artifact of floating point arithmetic, an illusion, a phantom of mathematical computation.