Problems with solve command.

1 次查看（过去 30 天）
mattyice
mattyice 2016-5-19
评论： John D'Errico 2016-5-26
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);

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2016-5-19
编辑：John D'Errico 2016-5-19
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. :)
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2016-5-25
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
John D'Errico 2016-5-26
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by