MATLAB Answers

sinh solver and numerical solution

66 views (last 30 days)
Mustafa AYDIN
Mustafa AYDIN on 14 Dec 2020
Commented: Rik on 18 Dec 2020
This question was flagged by Walter Roberson
Hi,
I need solution of equation below, I didn't do it correctly. Can you help me, Please?
equation must be x= bla bla
and second one
numerical solution for
pi=3.14;
kb=1.38E-23;
h=6.626E-34;
hbar=h/(2*pi);
e=1.602E-19;
T0=4.2;
T=10;
A0=1.26E3;
A=1.24E3;
B=10;
equ1 = ((T*sinh((2*(pi^2)*kb*T0*x)/hbar*e*B))/(T0*sinh((2*(pi^2)*kb*T*x)/hbar*e*B)))==(A/A0)
Thank you

  2 Comments

John D'Errico
John D'Errico on 14 Dec 2020
By the way, NEVER do things like define pi as 3.14. pi already exists in MATLAB. All this does is introduce additional error into the problem.
Rik
Rik on 18 Dec 2020
Why did you edit away parts of your question? That is extremely rude.

Sign in to comment.

Answers (2)

Walter Roberson
Walter Roberson on 14 Dec 2020
format long g
Q = @(v) sym(v);
Pi = Q(pi);
kb = Q(1.38)*Q(10)^(-23);
h = Q(6.626)*Q(10)^(-34);
hbar = h/(2*pi);
e = Q(1.602)*Q(10)^(-19);
T0 = Q(10);
T = Q(4.2);
A0 = Q(1.26);
A = Q(1.24);
B = Q(10);
syms x
equ1 = ((T*sinh((2*(Pi^2)*kb*T0*x)/hbar*e*B))/(T0*sinh((2*(Pi^2)*kb*T*x)/hbar*e*B)))==(A/A0)
equ1 = 
char(equ1)
ans = '(21*sinh((55269*x*pi^3)/41412500000))/(50*sinh((1160649*x*pi^3)/2070625000000)) == 62/63'
I am having trouble finding additional solutions beyond the one near -50.

  3 Comments

Mustafa AYDIN
Mustafa AYDIN on 14 Dec 2020
I face with same problem while solving this equation as well. A0 and A can be multiplied with 10^3
Mustafa AYDIN
Mustafa AYDIN on 14 Dec 2020
at least if I find symbolic solution of this equation, it will be a good job
Walter Roberson
Walter Roberson on 15 Dec 2020
There is are roots at about 0 +/- 8244.2223439771955066317321720362051816482188828761 * 1i -- a purely imaginary number
If you make the substitution x = 1i*XI then the sinh() turns into plain sin() and you are dealing with a real-valued problem. A taylor series and solve() of that permits you to get approximate locations, and you can then ask for a numeric search near those to get the more precise value. There is a clean root near the above locations if you plot the graph.
There are potentially other solutions, but not many. If you substitute x = XR + 1i * XI and do a 3D plot of the imaginary() part of left side of the equation, then you can see that the imaginary part crosses zero in two lines
The flat plane is 0 . The line near XR = 0 suggested to me that I try for a purely imaginary solution, which gave the results I indicated above. The line near XI = 0 hints that by symmetry there just might be two roots with near-zero imaginary components.

Sign in to comment.


John D'Errico
John D'Errico on 14 Dec 2020
Edited: John D'Errico on 14 Dec 2020
One mistake in what you are both doing: Do not use == there to form equ1. Just subtract off the right hand side constant.
The == prevents you from ever plotting anything. And that prevents you from seeing the fundamental shape of what results. And that prevents you from thinking about the real problem in what you wish to solve.
Next, recognize the result is just a ratio of a pair of sinh terms, less a constant. So if I look at the resulting expression, we have
Q = @(v) sym(v);
Pi = Q(pi);
kb = Q(1.38)*Q(10)^(-23);
h = Q(6.626)*Q(10)^(-34);
hbar = h/(2*pi);
e = Q(1.602)*Q(10)^(-19);
T0 = Q(10);
T = Q(4.2);
A0 = Q(1.26);
A = Q(1.24);
B = Q(10);
syms x
equ1 = ((T*sinh((2*(Pi^2)*kb*T0*x)/hbar*e*B))/(T0*sinh((2*(Pi^2)*kb*T*x)/hbar*e*B))) - (A/A0)
equ1 = 
Again, a ratio of simple sinh terms. so the fundamental shape is symmetric in x. But perhaps it is better to clear that fraction. I'll also divide by T. As long as we are not interested in the trivial solution at x==0, even if such a solution exists, we can do this (We know that sinh has no other zeros besides at zero, so this changes nothing else in the problem.)
equ2 = ((sinh((2*(Pi^2)*kb*T0*x)/hbar*e*B)) - (A/A0/T)*(T0*sinh((2*(Pi^2)*kb*T*x)/hbar*e*B)));
And, as well, the scaling of this thing is miserable. so a simple transformation will help to make the numbers better.
syms y
equ3 = subs(equ2,x,y*1e5);
vpa(equ3,5)
ans = 
You should agree I have done nothing that changes the fundamental nature of the problem. But what I did was make it easier to work with, and easier to understand the next steps because we can understand these simple numbers a little better.
We are looking for the root(s) of that thing. Again, recognize this is symmetric in x (as well as y now.) y == 0 is clearly a solution, but a spurious one created when we cleared the fraction. So ignore that solution.
Do any other roots exist? In fact, the first sinh term will dominate the result, because sinh grows in magnitude so rapidly. One trick is to consider the Taylor series of equ3, remembering this old rule as we go:
vpa(taylor(equ3,'order',50),3)
ans = 
The point is to carefully look at the coefficients of every term in that series. I generated a lot of terms here, just to convince you that none of the coefficients are negative, at least out that far. So if we have a root at zero, and all coefficients of the Taylor polynomial are positive, then we can assert that no other root exists. (A useful fact to remember - The Taylor series for sinh(u) is convergent for all u. Therefore, if we can prove no root can exist for the Taylor polynomial in this case, then no root can exist for the original function either.) And since the result is a symmetric function in the unknown (there may only be odd order terms in the Taylor polynomial, think about it - this is just the difference of a pair of sinh terms), then no negative root exists either.
So next, I suppose we could spend some time to compute the general coefficient in the Taylor expansion of this function, then to prove it must always be positive. I think I can do that, based on the Taylor series expansion for sinh.
It is pretty simple really. The Taylor series of sinh(u) has only odd terms in it, and the nth term is just u^(2*n+1)/factorial(2*n+1). But since the sinh(4.1381*y) in our expansion will dominate the sinh(1.738*y) term, the difference of the two will always have a positive coefficient, for large n as well as small.
And again, as long as the resulting polynomial always has positive coefficients (with a zero constant term) then there are no positive roots, nor are there any negative roots due to symmetry.
One other comment: I see you think there was a root near -50? That was at best a spurious one, due to numerical problems since we have shown no roots can exist. When x is small in your original version equ1, you were dividing by a VERY small number when x is approximately -50. And that will cause numerical problems. That is not an issue in what I did. As long as EVERY coefficient in the Taylor series expansion of equ3 is uniformly positive (they are) then no root can exist other than zero.

  0 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!