Numerical solutions of equation

Question

NILESH PANDEY 2018-7-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/411569-numerical-solutions-of-equation

评论： madhan ravi 2018-7-23

采纳的回答： Walter Roberson

在 MATLAB Online 中打开

I'm trying to solve the following equation with respect to b

I need to find the value of b code is

clear all
E=1.5e8;
syms b
 a0=-4e7/400;
T0=600;
T=300;
temp=1/(2*a0*(T0-T));
X=5050700286489334129256038400000/(4349164374701853503175341475971*b) + exp(erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2)^2)*exp(-(2*b*erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2) - b^2*1i)^2/(4*b^2))*(50507002864893341292560384/(4349164374701853503175341475971*b) + (50507002864893341292560384*(erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b)))/(4349164374701853503175341475971*b*((erf(pi/b) - (101014005729786682585120768*E + 4040560229191467303404830720000000)/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2))) + exp(erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2)^2)*exp(-(2*b*erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2) + b^2*1i)^2/(4*b^2))*(50507002864893341292560384/(4349164374701853503175341475971*b) + (50507002864893341292560384*(erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b)))/(4349164374701853503175341475971*b*((erf(pi/b) - (101014005729786682585120768*E + 4040560229191467303404830720000000)/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2))) + (5050700286489334129256038400000*(erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b)))/(4349164374701853503175341475971*b*((erf(pi/b) - (101014005729786682585120768*E + 4040560229191467303404830720000000)/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)) + 7737125245533627/9671406556917033397649408
 vpasolve(X==temp,b)

but matlab(2016 version) gives the following error:

ans =
Empty sym: 0-by-1

equation also contain imaginary part which is written 1i in the expression of X

please help me for the same code for the calculation of X given below please have a look maybe this can help, code for the calculation of X is:

   clear all
syms  E b;
N0=1e6;
E0=1e8;
Eeff=-0.4e8;
a=19.3567;
p=0.1;
ed=0.4e-9;
efe=4.5*1e6*8.854e-12;
x=2*efe*(E-Eeff)/(a*b*p*sqrt(pi));
xm=erf(pi/b);
delta=0.05;
H_f=x-0.5*(x-xm+sqrt((x-xm).^2+4*delta));
theta=b*erfinv(H_f);
N1=N0*p;
J=N1*erf(pi/b);
p1=0.5*N1*exp(-b*b*0.5);
P_diel=2*ed*(E-E0);
p2=erf((2*theta+i*b*b)/(2*b))+erf((2*theta-i*b*b)/(2*b))+N1*erf(theta/b);%+P_diel;
P=p1+p2+P_diel;
X= diff(P,E);

5 个评论
显示 3更早的评论隐藏 3更早的评论

NILESH PANDEY 2018-7-23

Thanks for your time,I'll try some other way maybe series expansion. If I found the solution I'll share it.

madhan ravi 2018-7-23

Yes sure do ,good luck.

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2018-7-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/411569-numerical-solutions-of-equation#answer_329894

在 MATLAB Online 中打开

%clear all
syms  E b;
N0=1e6;
E0=1e8;
Eeff=-0.4e8;
a=19.3567;
p=0.1;
ed=0.4e-9;
efe=4.5*1e6*8.854e-12;
x=2*efe*(E-Eeff)/(a*b*p*sqrt(pi));
xm=erf(pi/b);
delta=0.05;
H_f=x-0.5*(x-xm+sqrt((x-xm).^2+4*delta));
theta=b*erfinv(H_f);
N1=N0*p;
J=N1*erf(pi/b);
p1=0.5*N1*exp(-b*b*0.5);
P_diel=2*ed*(E-E0);
p2=erf((2*theta+1i*b*b)/(2*b))+erf((2*theta-1i*b*b)/(2*b))+N1*erf(theta/b);%+P_diel;
P=p1+p2+P_diel;
X= diff(P,E);
Eval=1.5e8;
a0=-4e7/400;
T0=600;
T=300;
temp=1/(2*a0*(T0-T));
eqn = subs(X, E, Eval) - temp;
B = linspace(-10,10);
Y = double( subs(eqn, b, B) );
plot(B, Y)
vpasolve(eqn, b, [6 7])

answer is approximately 6.688061088729162561767129937816

2 个评论
显示无隐藏无

NILESH PANDEY 2018-7-23

在 MATLAB Online 中打开

I tried to run your code but Matlab gives the following error:

    Error using symengine
DOUBLE cannot convert the input expression into a double array.
Error in sym/double (line 613)
        Xstr = mupadmex('symobj::double', S.s, 0);
Error in sol (line 33)
Y = double( subs(eqn, b, B ) );

there is an error in line 33:

Y = double( subs(eqn, b, B ) );

is it because of Matlab version I'm using Matlab-2016?

NILESH PANDEY 2018-7-23

Thanks Sir for answer now I see above error is for the plotting, so I commented above lines in Matlab now matlab is able to solve above equation. Thanks again for your help and time!

请先登录，再进行评论。