vpasolve and fimplicit not returning anything

I am trying to use vpasolve to numerically solve a equation, which is given as (I have explicitly written out the langevin function, since I could not find an inbuilt function for this). I have to solve for y, given a range of values of and . To make sure things work, I first fix , . The vpasolve function returns Empty sym.

I have also tried to use the fimplicit function to plot the equation, and the result is correct albeit only for larger values of x. I expect a continous curve starting from .

My attempt is as follows

ebar=10;
syms y
x=0.5;
solve(coth(12.*ebar./x.*(y.^(-5)-y.^(-11)))- ...
     1./(12.*ebar./x.*(y.^(-5)-y.^(-11)))- x./y == 0,y,1);


%Try fimplicit function
f = @(x,y) coth(12.*ebar./x.*(y.^(-5)-y.^(-11)))- ...
    1./(12.*ebar./x.*(y.^(-5)-y.^(-11)))- x./y;
fimplicit(f,[0 2 0.5 2])

Answer 1

[Walter Roberson]

ebar=10;
syms y
x=0.5;
eqn = coth(12.*ebar./x.*(y.^(-5)-y.^(-11))) - ...
     1./(12.*ebar./x.*(y.^(-5)-y.^(-11)))- x./y == 0

eqn = 

eq1 = matlabFunction(lhs(eqn))

eq1 = function_handle with value:
    @(y)coth(1.0./y.^5.*2.4e+2-1.0./y.^11.*2.4e+2)-1.0./(1.0./y.^5.*2.4e+2-1.0./y.^11.*2.4e+2)-1.0./(y.*2.0)

ysoln = fzero(eq1, 3)

ysoln = 3.5454

ysol = vpasolve(eqn, y, 3)

ysol = 

3.5453801763200780562815540493667

%fplot([lhs(eqn)-rhs(eqn), 0], [-4 4])
%xticks(-5:5)

[Walter Roberson]

ebar=10;
syms y
x=0.5;
eqn = coth(12.*ebar./x.*(y.^(-5)-y.^(-11))) - ...
     1./(12.*ebar./x.*(y.^(-5)-y.^(-11)))- x./y == 0

eqn = 

fplot([lhs(eqn)-rhs(eqn), 0], [0.01 4], 'MeshDensity', 7)
xticks(0:4)
ylim([-5 5])

The crossing at 1 is a discontinuity, the real crossing is near 3

[Walter Roberson]

ebar=10;
syms y
x=0.1;
eqn = coth(12.*ebar./x.*(y.^(-5)-y.^(-11))) - ...
     1./(12.*ebar./x.*(y.^(-5)-y.^(-11)))- x./y == 0

eqn = 

fplot([lhs(eqn)-rhs(eqn), 0], [0.8 10], 'MeshDensity', 7)
xticks(0:10)
ylim([-1 1])

vpasolve(eqn, y, 5)

ans = 

7.952510786446735090886989392305

[Walter Roberson]

Work in progress, not complete

ebar=10;
syms x positive
syms y
eqn = coth(12.*ebar./x.*(y.^(-5)-y.^(-11))) - ...
     1./(12.*ebar./x.*(y.^(-5)-y.^(-11)))- x./y == 0

eqn = 

eqn1 = lhs(eqn) - rhs(eqn);
eqn2 = subs(eqn1, y, 100001/100000)

eqn2 = 

solve(eqn2)

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

ans = 

0.048952698891728760676529060961778

eqn5 = subs(eqn1, x, 0.05)

eqn5 = 

vpasolve(eqn5, y, 1.001)

ans = 

1.0000104331446060728615627131999

[Walter Roberson]

As x gets above 1, the corresponding y gets very close to 1 . If you are working in numeric form, the round-off error cannot distinguish the value from the singularity at y = 1

As x reduces, such as getting less than 0.1, the y solution gets larger.

You can substitute any given large y in to the equation and solve the result for x, demonstrating that for y sufficiently large, there is a small x that solves the equation -- also implying that if your x is sufficiently small that your y could be arbitrarily large.

For x = 1, there are four solutions, +/- 2.44177 and +/- 1.04609

[Mohammad Ansari]

At at value looks right. I am actually trying to reproduce the following figure, so for I expect the value of to be near 1, and not large. The original equation which the paper mentions is , where (the Langevin function) and is its inverse. I have rearranged this written it out in the original question post, but writing it out explicitly to make sure that any numerical finesse is not lost in inversion of this equation.

Thank you for helping!

and

[Walter Roberson]

Look at the green line. Near 0.65 it has reached 3 and clearly is going up relatively steeply .

I showed above that for 0.1 it reaches close to 8.

By substituting in large y such as 1000 you can solve for x and that x will be in range 0 to 1. It follows that as x approaches 0 that the y solution grows arbitrarily large .

[Mohammad Ansari]

I see what you mean now.

But shouldn't y have another solution (as seen from the figure I attached), which is closer to y=1? eg at y=1, you showed that 1.04 and 2.44 are both solutions (and this can be seen from the figure I attached), can I get the smaller solution for ?

[Mohammad Ansari]

I was able to get the near 1 solution by properly choosing the initial guess as 1.0001. I guess this was the part I was missing, choosing the initial guess as 1 returns Empty Sym, while choosing 1.0001 does the job. Thank you so much for your help.