Numerical solutions of equation

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);

[madhan ravi]

It’s not an error, maybe there is no solution or matlab is not able to compute .

[NILESH PANDEY]

Thanks for the suggestion. I'm giving the code by which I have calculated X please have a look maybe this will help Code for the X

    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);

[madhan ravi]

This looks fine but when I tried to solve for b it didn’t show the explicit solution , in this case vpasolve was used it is the best numerical solver in matlab so unfortunately no appropriate answer for this equation:/

[NILESH PANDEY]

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

[madhan ravi]

Yes sure do ,good luck.

Answer 1

[Walter Roberson]

%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

[NILESH PANDEY]

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]

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!