You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
I have a nonlinear equation ,, which i didn't know how to put it in matlab
4 views (last 30 days)
Show older comments
the equation
x exp x = [(segma / KB * T ) ^ 2 - x ]* (tr/ttr) * exp(E0 - Ea)/KB * T
while :
segma = 13 , KB =8.617*10^-5 , tr= 250 , ttr= 0.027 , E0=1.185 , Ea=E0+0.073
T=0:300
and the other one is
E(T)= E0 -(0.00000048 * T^2) /(270 + T) - x * KB *T
while
E0=1.185 , KB =8.617*10^-5 and also T =0:300
can any one help me with it please
2 Comments
beso ss
on 1 Mar 2018
i need to give a result for each x ,, that is means for every T there is an x for it ,,, if i but the other elemnts to the equation it will be
x exp x = [ ( 13 / 8.617*10^-5 * T )^2 - x ] * ( 250/0.027 ) * exp[ ( 0.073 ) / 8.617*10^-5 * T ]
there is still x and T but T change from 0 to 300
Answers (6)
John D'Errico
on 2 Mar 2018
Edited: John D'Errico
on 2 Mar 2018
You responded to my question without explaining what E(T) has to do with anything.
So is it clear? Not fully. It seems you have written TWO equations, one for E(T) which seems to be used for nothing.
The first equation can be solved (IF a solution exists) using fzero. Thus, for EACH value of T, you will solve the problem independently. A loop would suffice, and there is no real reason to do anything more sophisticated than a loop.
T = 0:300;
X = zeros(size(T));
syms x
for n = 1:numel(T)
Eqn1 = x*exp(x) == [(segma / KB * T(n) ) ^ 2 - x ]* (tr/ttr) * exp(E0 - Ea)/KB * T(n);
X(n) = double(vpasolve(Eqn1,x));
end
The above solution could also have been trivially written to use fzero.
14 Comments
beso ss
on 2 Mar 2018
E(T) is the energy depend on the T temperature ,,, varshni law ... and it is depend on x for each value of x in any temperature then E(T) also will change
i but the same what you write
T = 0:300; X = zeros(size(T)); syms x for n = 1:numel(T) Eqn1 = x*exp(x) ==(( 13 / 8.617*10^-5 * T(n) )^2 - x ) * ( 250/0.027 ) * exp ( 0.073 ) / 8.617*10^-5 * T(n) ;
X(n) = double(vpasolve(Eqn1,x));
end
but this massage appear ; Undefined function 'syms' for input arguments of type 'char'.
Error in xxxxxxxxxx (line 3) syms x
i don't know whats wrong ?????????
beso ss
on 2 Mar 2018
i write ver and the answer was
MATLAB Version 8.0 (R2012b) Aerospace Blockset Version 3.10 (R2012b) Aerospace Toolbox Version 2.10 (R2012b) Bioinformatics Toolbox Version 4.2 (R2012b) Communications System Toolbox Version 5.3 (R2012b) Computer Vision System Toolbox Version 5.1 (R2012b) Control System Toolbox Version 9.4 (R2012b) Curve Fitting Toolbox Version 3.3 (R2012b) DO Qualification Kit Version 2.0 (R2012b) DSP System Toolbox Version 8.3 (R2012b) Data Acquisition Toolbox Version 3.2 (R2012b) Database Toolbox Version 4.0 (R2012b) Datafeed Toolbox Version 4.4 (R2012b) Econometrics Toolbox Version 2.2 (R2012b) Embedded Coder Version 6.3 (R2012b) Filter Design HDL Coder Version 2.9.2 (R2012b) Financial Instruments Toolbox Version 1.0 (R2012b) Financial Toolbox Version 5.0 (R2012b) Fixed-Point Toolbox Version 3.6 (R2012b) Fuzzy Logic Toolbox Version 2.2.16 (R2012b) Global Optimization Toolbox Version 3.2.2 (R2012b) HDL Coder Version 3.1 (R2012b) HDL Verifier Version 4.1 (R2012b) IEC Certification Kit Version 3.0 (R2012b) Image Acquisition Toolbox Version 4.4 (R2012b) Image Processing Toolbox Version 8.1 (R2012b) Instrument Control Toolbox Version 3.2 (R2012b) MATLAB Builder EX Version 2.3 (R2012b) MATLAB Builder JA Version 2.2.5 (R2012b) MATLAB Builder NE Version 4.1.2 (R2012b) MATLAB Coder Version 2.3 (R2012b) MATLAB Compiler Version 4.18 (R2012b) MATLAB Report Generator Version 3.13 (R2012b) Mapping Toolbox Version 3.6 (R2012b) Model Predictive Control Toolbox Version 4.1.1 (R2012b) Model-Based Calibration Toolbox Version 4.5 (R2012b) Neural Network Toolbox Version 8.0 (R2012b) OPC Toolbox Version 3.1.2 (R2012b) Optimization Toolbox Version 6.2.1 (R2012b) Parallel Computing Toolbox Version 6.1 (R2012b) Partial Differential Equation Toolbox Version 1.1 (R2012b) Phased Array System Toolbox Version 1.3 (R2012b) RF Toolbox Version 2.11 (R2012b) Real-Time Windows Target Version 4.1 (R2012b) Robust Control Toolbox Version 4.2 (R2012b) Signal Processing Toolbox Version 6.18 (R2012b) xPC Target Version 5.3 (R2012b) xPC Target Embedded Option Version 5.3 (R2012b)
iam sorry this is the first time i work with matlab so i don't know ,,, how i install it
Torsten
on 2 Mar 2018
As the symbolic toolbox is not listed as licensed product, you'll have to buy and/or install it.
beso ss
on 2 Mar 2018
it's works thank you so muchh but i have another two equations please
the first one :
E = 1.185 - (0.00000048 * T^2) /(270 + T)
while T change from 0 to 300 ,,,,,, T = 0:300 I want E in every time T change ,, and i want to plot it
the other on is related with the both equations E= 1.185 - (0.00000048 * T^2) /(270 + T) - x * 8.617*10^-5 * T
thats means every time T will change also x will change and then E will change and also i want to plot it
im sooo sorry but im really i don't know anything about matlab , and i am thank you sooo much .
John D'Errico
on 3 Mar 2018
Please don't add answers just to comment. Moved your answer to a comment by @beso...
for the first equation
E = 1.185 - (0.00000048 * T^2) /(270 + T)
i type this :
T=linspace(0,300,300) for i= 1 : 300
E(i)=1.185-(0.00000048*T(i)^2)/270+T(i);
and when i want to plot it i type :
plot(T,E(i))
but the figure was sooo different from what i have in the paper ,,, i don't know what is the wrong can anyone answer me please
and i didn't know how i type the other one ;
E= 1.185 - (0.00000048 * T^2) / (270 + T) - x * 8.617*10^-5 * T
John D'Errico
on 3 Mar 2018
Edited: John D'Errico
on 3 Mar 2018
So you need to learn to use plot. You need to learn basic tools of MATLAB. (Which means you need to read the getting started tutorials.)
T=linspace(0,300,300);
E = 1.185-(0.00000048*T.^2)./(270+T);
plot(T,E)
That is all you need, and nothing more. There is no need to use a loop.
Note my use of the .^ operator for element-wise power. Sometimes you will also need to use .* and ./ for element-wise multiplications or divisions between two vectors. As you can see, I used ./ there too.
You use those operators whenever you need to perform an operation on the elements of vectors or arrays.
And, once you have solved for the vector of values x, one for every element of T, you do the same operations on the real line for E.
E = 1.185 - (0.00000048 * T.^2) ./ (270 + T) - 8.617e-5*X.*T;
As you can see here, I used the ./ operator, because I was dividing the elements of one vector into another vector of the same shape and size. I used .^ to square the elements of the vector T. I used X.*T to multiply the elements of those vectors. This works because X and T are the same size, so both are row vectors.
Also note my use of a semi-colon at the end of the line.
Finally, there is no need to use .* to multiply be a scalar value, because * works already for that.
Finally, I did write
8.617*10^-5
as
8.617e-5
The latter is more efficient since there is no need to raise 10 to a power, as well as easier to read and write.
beso ss
on 3 Mar 2018
i tried this sultion before but i thought may be a have to use a loop
but what about this equation E= 1.185 - (0.00000048 * T^2) /(270 + T) - x * 8.617*10^-5 * T it is contain x here ,, what we do it before how can i mix them together ???
beso ss
on 3 Mar 2018
Edited: Walter Roberson
on 4 Mar 2018
T=linspace(0,300,300);
E1 = 1.185-(0.00000048*T.^2)./(270+T);
plot(T,E1)
T = 0:300;
syms x
for n = 1:numel(T)
Eqn1 = x*exp(x) ==(( 13 / 8.617*10^-5 * T(n) )^2 - x ) * ( 250/0.027 ) * exp ( 0.073 ) / 8.617*10^-5 * T(n) ;
X(n) = double(vpasolve(Eqn1,x));
E2 = 1.185 - (0.00000048 * T.^2) ./ (270 + T) - 8.617e-5*X.*T;
plot(T,E2)
end
is it correct ???
the value of E1 and E2 take the same result for all T ,,, i mean the answer in each one of them is 1.185 for T=0:300
i don't know why
John D'Errico
on 4 Mar 2018
Sigh. Are you still stumbling on this?
You defined E2 inside the loop. But you don't want to do that.
You alsos seem to have defined T twice here, once using linspace, once using colon. But one of them will be of length 300, the other is of length 301.
You are just throwing random lines of code around without thinking about what you are doing now.
T = 0:300;
E1 = 1.185-(0.00000048*T.^2)./(270+T);
syms x
X = zeros(size(T));
for n = 1:numel(T)
Eqn1 = x*exp(x) ==(( 13 / 8.617*10^-5 * T(n) )^2 - x ) * ( 250/0.027 ) * exp ( 0.073 ) / 8.617*10^-5 * T(n) ;
X(n) = double(vpasolve(Eqn1,x));
end
E2 = 1.185 - (0.00000048 * T.^2) ./ (270 + T) - 8.617e-5*X.*T;
Yes. E2 is essentially 1.185, or very close to that value. But LOOK AT THE NUMBERS IN YOUR VECTORS. THINK ABOUT WHAT YOU ARE DOING. Don't just randomly type lines of code and expect a useful result.
What is X?
X
X =
Columns 1 through 13
0 2.6008e-12 2.0571e-11 6.8653e-11 1.6093e-10 3.1089e-10 5.3141e-10 8.3484e-10 1.233e-09 1.7372e-09 2.3583e-09 3.1067e-09 3.9924e-09
Columns 14 through 26
5.0249e-09 6.2136e-09 7.5671e-09 9.0941e-09 1.0803e-08 1.2701e-08 1.4796e-08 1.7095e-08 1.9605e-08 2.2334e-08 2.5287e-08 2.8471e-08 3.1891e-08
Columns 27 through 39
3.5555e-08 3.9466e-08 4.3631e-08 4.8055e-08 5.2743e-08 5.77e-08 6.2931e-08 6.8439e-08 7.423e-08 8.0308e-08 8.6677e-08 9.3341e-08 1.003e-07
Columns 40 through 52
1.0757e-07 1.1514e-07 1.2302e-07 1.3121e-07 1.3972e-07 1.4855e-07 1.5771e-07 1.6718e-07 1.7699e-07 1.8713e-07 1.976e-07 2.0841e-07 2.1956e-07
Columns 53 through 65
2.3105e-07 2.4288e-07 2.5506e-07 2.6759e-07 2.8047e-07 2.937e-07 3.0729e-07 3.2124e-07 3.3555e-07 3.5022e-07 3.6525e-07 3.8064e-07 3.9641e-07
Columns 66 through 78
4.1254e-07 4.2904e-07 4.4592e-07 4.6317e-07 4.8079e-07 4.9879e-07 5.1717e-07 5.3593e-07 5.5507e-07 5.746e-07 5.945e-07 6.1479e-07 6.3547e-07
Columns 79 through 91
6.5654e-07 6.7799e-07 6.9984e-07 7.2208e-07 7.4471e-07 7.6773e-07 7.9115e-07 8.1496e-07 8.3917e-07 8.6378e-07 8.8879e-07 9.1419e-07 9.4e-07
Columns 92 through 104
9.6621e-07 9.9282e-07 1.0198e-06 1.0473e-06 1.0751e-06 1.1033e-06 1.1319e-06 1.161e-06 1.1904e-06 1.2203e-06 1.2506e-06 1.2813e-06 1.3123e-06
Columns 105 through 117
1.3439e-06 1.3758e-06 1.4081e-06 1.4408e-06 1.474e-06 1.5076e-06 1.5416e-06 1.576e-06 1.6108e-06 1.646e-06 1.6817e-06 1.7178e-06 1.7543e-06
Columns 118 through 130
1.7912e-06 1.8285e-06 1.8663e-06 1.9045e-06 1.9431e-06 1.9821e-06 2.0215e-06 2.0614e-06 2.1017e-06 2.1424e-06 2.1835e-06 2.2251e-06 2.2671e-06
Columns 131 through 143
2.3095e-06 2.3524e-06 2.3956e-06 2.4393e-06 2.4835e-06 2.528e-06 2.573e-06 2.6184e-06 2.6642e-06 2.7105e-06 2.7572e-06 2.8043e-06 2.8519e-06
Columns 144 through 156
2.8999e-06 2.9483e-06 2.9971e-06 3.0464e-06 3.0961e-06 3.1463e-06 3.1968e-06 3.2478e-06 3.2993e-06 3.3512e-06 3.4035e-06 3.4562e-06 3.5094e-06
Columns 157 through 169
3.563e-06 3.617e-06 3.6715e-06 3.7264e-06 3.7818e-06 3.8376e-06 3.8938e-06 3.9504e-06 4.0075e-06 4.0651e-06 4.123e-06 4.1814e-06 4.2403e-06
Columns 170 through 182
4.2995e-06 4.3593e-06 4.4194e-06 4.48e-06 4.541e-06 4.6025e-06 4.6644e-06 4.7268e-06 4.7895e-06 4.8528e-06 4.9164e-06 4.9805e-06 5.0451e-06
Columns 183 through 195
5.11e-06 5.1755e-06 5.2413e-06 5.3076e-06 5.3744e-06 5.4415e-06 5.5092e-06 5.5772e-06 5.6457e-06 5.7147e-06 5.7841e-06 5.8539e-06 5.9242e-06
Columns 196 through 208
5.9949e-06 6.066e-06 6.1376e-06 6.2097e-06 6.2822e-06 6.3551e-06 6.4284e-06 6.5023e-06 6.5765e-06 6.6512e-06 6.7263e-06 6.8019e-06 6.8779e-06
Columns 209 through 221
6.9544e-06 7.0313e-06 7.1087e-06 7.1865e-06 7.2647e-06 7.3434e-06 7.4225e-06 7.5021e-06 7.5821e-06 7.6626e-06 7.7435e-06 7.8249e-06 7.9067e-06
Columns 222 through 234
7.9889e-06 8.0716e-06 8.1547e-06 8.2383e-06 8.3224e-06 8.4068e-06 8.4917e-06 8.5771e-06 8.6629e-06 8.7492e-06 8.8359e-06 8.923e-06 9.0106e-06
Columns 235 through 247
9.0986e-06 9.1871e-06 9.2761e-06 9.3654e-06 9.4553e-06 9.5455e-06 9.6362e-06 9.7274e-06 9.819e-06 9.9111e-06 1.0004e-05 1.0097e-05 1.019e-05
Columns 248 through 260
1.0284e-05 1.0378e-05 1.0473e-05 1.0568e-05 1.0664e-05 1.076e-05 1.0856e-05 1.0953e-05 1.1051e-05 1.1149e-05 1.1247e-05 1.1346e-05 1.1445e-05
Columns 261 through 273
1.1544e-05 1.1645e-05 1.1745e-05 1.1846e-05 1.1948e-05 1.205e-05 1.2152e-05 1.2255e-05 1.2358e-05 1.2462e-05 1.2566e-05 1.267e-05 1.2775e-05
Columns 274 through 286
1.2881e-05 1.2987e-05 1.3093e-05 1.32e-05 1.3307e-05 1.3415e-05 1.3523e-05 1.3632e-05 1.3741e-05 1.3851e-05 1.3961e-05 1.4071e-05 1.4182e-05
Columns 287 through 299
1.4293e-05 1.4405e-05 1.4517e-05 1.463e-05 1.4743e-05 1.4857e-05 1.4971e-05 1.5085e-05 1.52e-05 1.5315e-05 1.5431e-05 1.5548e-05 1.5664e-05
Columns 300 through 301
1.5782e-05 1.5899e-05
So all elements of X lie between 0 and 1.6e-5. They are tiny numbers.
What happens when for T = 0:300, you substitute a REALLY small value for X in the expression for E2?
T is not that large. S square it, then multiply it by 0.00000048, and divide by (270+T). It will be tiny.
Next, what is X*T? Multiply that by another tiny number, thus 8.617e-5. What do you expect there?
So, to compute E2, you add two very small numbers to 1.185. The result is very close to 1.185. Why are you even remotely surprised?
John D'Errico
on 4 Mar 2018
Edited: John D'Errico
on 4 Mar 2018
They Are NOT the same result!!!!!!!
format long g
E1 = 1.185-(0.00000048*T.^2)./(270+T);
E1(1:5)
ans =
1.185 1.18499999822878 1.18499999294118 1.18499998417582 1.1849999719708
E2(1:5)
ans =
1.185 1.18499999822878 1.18499999294117 1.18499998417581 1.18499997197075
Are they the same? No. Are they close as hell to each other? Of course. Because those other terms are SO small.
Consider E1.
It was computed from
E1 = 1.185-(0.00000048*T.^2)./(270+T);
So plot that second term in E1.
plot(T,(0.00000048*T.^2)./(270+T))
Look carefully at the y axis. Do you see on top the 10^-5 there?
That means all the elements in that plot are on the order of 8*10^-5, or SMALLER. So if we add 1.185 to a number that is SMALLER than 0.00001, and usually MUCH smaller than that, what do you expect to see?
How about E2? E2 is the same as E1. Except is has that term with X in it.
E2 = 1.185 - (0.00000048 * T.^2) ./ (270 + T) - 8.617e-5*X.*T;
Plot that last term now. Remember to look carefully at the y axis. Remember to look at the exponent attached to that axis. Do you see the 10^-7 there? That means that every number in that plot was smaller than 8e-7, but that the elements on the far left are FAR smaller than that.
plot(T,8.617e-5*X.*T)
Here are the first 5 elements of that term.
8.617e-5*X(1:5).*T(1:5)
ans =
0 2.24110912198165e-16 3.54526286107169e-15 1.77473849967693e-14 5.5470798356365e-14
Now, when you subtract a number of that size from 1.185, what do you expect to see?
Think abut what you are doing. Look carefully at the numbers involved.
beso ss
on 4 Mar 2018
can anyone answer me please
2 Comments
beso ss
on 4 Mar 2018
im sorry :( ,,,, thank you sooooo much ,,, some times i am so rush and i didn't git the hall answer in the same time ,, it comes separate ,,,
thank you all again very much
beso ss
on 17 Mar 2018
Edited: Walter Roberson
on 17 Mar 2018
Hi , it is me again
can i use another way to have numerically solve for the equation
syms x
X = zeros(size(T));
for n = 1:numel(T)
Eqn1 = x*exp(x) ==(( 13 / 8.617*10^-5 * T(n) )^2 - x ) * ( 250/0.027 ) * exp ( 0.073 ) / 8.617*10^-5 * T(n) ;
X(n) = double(vpasolve(Eqn1,x));
end
becouse in the plot there is a mistake i dont know way
the original equation is
x*exp(x) = (((sigma/(KB*T(i)))^2-x)*(tautr/taur)*exp(deltaE/(KB*T(i)))
where
sigma = 13*10^-3; %eV
deltaE = -0.073; %eV
E0 = 1.185; %eV
tautr = 0.027;
taur =250;
KB= 8.617*10^-5; %eV/K
T= 0:300
54 Comments
beso ss
on 17 Mar 2018
Edited: Walter Roberson
on 17 Mar 2018
i use another way to solve it but also the plot not correct i donnt know what is the wrong with it
syms sigma deltaE E0 tautr taur KB theta alfa
sigma = 13*10^-3; %eV
deltaE = -0.073; %eV
E0 = 1.185; %eV
tautr = 0.027;
taur =250;
KB= 8.617*10^-5; %eV/K
theta = 270; %K
alfa = 0.00048; %eV/K
T=0:300;
X=zeros(1,numel(T));
for i=1:numel(T)
syms x
X(i)=vpasolve(((sigma/(KB*T(i)))^2-x)*(tautr/taur)*exp(deltaE/(KB*T(i)))-x*exp(x)==0,x);
clear x
end
Walter Roberson
on 17 Mar 2018
You can switch to using a numeric solver based upon one side of the equation minus the other side.
There appear to be two roots, one in the negative hundreds and the other one extremely close to 0, unless T(i) is negative or is above roughly 10^368. Either way you are going to have a lot of difficulty with loss of numeric precision due to the exp(deltaE/(KB*T[i])) with deltaE being negative and KB being small, which together give you exp() of a large negative number divided T(i)
beso ss
on 17 Mar 2018
but when i plot it ,, there is somthing wrong i don't know where is the wrong in my solution
Walter Roberson
on 17 Mar 2018
Is it still the case that T=0:300 ?
It would help if you were to describe what you perceive as being "wrong" with your solution.
beso ss
on 17 Mar 2018
Edited: Walter Roberson
on 17 Mar 2018
i want spacily x because there are two equations the first one
E1=E0 -( alfa * T(i)^2 )/(theta + T(i))
and the other one is ;
E2=E0 -( alfa * T(i)^2 )/(theta + T(i)) - X * KB * T(i)
the difference between them is the part with x and that's make the plot different between them
it have to be like this
Walter Roberson
on 17 Mar 2018
Then your equations must be wrong. Those equations do not lead to anything symmetrical like you show.
With the equations you have, you should plot out to about T = 4600.
beso ss
on 17 Mar 2018
when i put T=3000 it works with x correct , but it didn't work when i plot E2 ???
:(
Walter Roberson
on 17 Mar 2018
I do not have values for alfa or theta .
Your X is a vector, so your E2 would be a vector for each value of i. Should the reference to X be a reference to X(i) ?
It would help if you were to present full current code.
beso ss
on 17 Mar 2018
Edited: Walter Roberson
on 17 Mar 2018
syms sigma deltaE E0 tautr taur KB theta alfa
sigma = 13*10^-3; %eV
deltaE = -0.073; %eV
E0 = 1.185; %eV
tautr = 0.027;
taur =250;
KB= 8.617*10^-5; %eV/K
theta = 270; %K
alfa = 0.00048; %eV/K
T=0:300;
X=zeros(1,numel(T));
for i=1:numel(T)
E1=E0 -( alfa * T(i)^2 )/(theta + T(i))
syms x
X(i)=vpasolve(((sigma/(KB*T(i)))^2-x)*(tautr/taur)*exp(deltaE/(KB*T(i)))-x*exp(x)==0,x);
clear x
E2=E0 -( alfa * T(i)^2 )/(theta + T(i)) - X * KB * T(i);
end
plot(T,E2)
hould the reference to X be a reference to X(i) ?yes
Walter Roberson
on 18 Mar 2018
Try this code:
syms sigma deltaE E0 tautr taur KB theta alfa
sigma = 13*10^-3; %eV
deltaE = -0.073; %eV
E0 = 1.185; %eV
tautr = 0.027;
taur =250;
KB= 8.617*10^-5; %eV/K
theta = 270; %K
alfa = 0.00048; %eV/K
T = linspace(0,4600);
X = zeros(size(T));
E2 = zeros(size(T));
syms x
for i=1:numel(T)
E1(i) = E0 -( alfa * T(i)^2 )/(theta + T(i));
X(i) = vpasolve(((sigma/(KB*T(i)))^2-x)*(tautr/taur)*exp(deltaE/(KB*T(i)))-x*exp(x)==0,x);
E2(i) = E0 -( alfa * T(i)^2 )/(theta + T(i)) - X(i) * KB * T(i);
end
subplot(1,3,1)
plot(T,E1)
title('E1')
subplot(1,3,2)
plot(T,E2)
title('E2')
subplot(1,3,3)
plot(T,X)
title('X')
Notice that the magnitude of E1 is much larger than the magnitude of X, and KB * T is less than 1 until T somewhere larger than 11500, so subtracting off X * KB * T is nearly the same as subtracting 0 as far as a plot can tell.
Walter Roberson
on 20 Mar 2018
Whatever. Change the 4600 to 300 then. But in your drawing of what the plot should look like, you want your plot to rise to a peak and then fall again, and with the equations that you gave, the fall cannot happen until T is larger than 300.
beso ss
on 21 Mar 2018
Edited: Walter Roberson
on 21 Mar 2018
now i try to change all the input ,, and when i run the matlab it is abear that x=0 in every point
syms sigma deltaE E0 tautr taur KB theta alfa
sigma = 1.1*10^-3; %eV
deltaE = -80; %eV
E0 = 3.184; %eV
tautr = 0.08;
taur =1;
KB= 8.617*10^-5; %eV/K
theta = 630; %K
alfa = 0.0011; %eV/K
T = linspace(0,300);
X = zeros(size(T));
E2 = zeros(size(T));
syms x
for i=1:numel(T)
E1(i) = E0 -( alfa * T(i)^2 )/(theta + T(i));
X(i) = vpasolve(((sigma/(KB*T(i)))^2-x)*(tautr/taur)*exp(deltaE/(KB*T(i)))-x*exp(x)==0,x);
E2(i) = E0 -( alfa * T(i)^2 )/(theta + T(i)) - X(i) * KB * T(i);
end
subplot(1,3,1)
plot(T,E1)
title('E1')
subplot(1,3,2)
plot(T,E2)
title('E2')
subplot(1,3,3)
plot(T,X)
title('X')
Walter Roberson
on 21 Mar 2018
deltaE/KB is -928397.354067541 so exp(deltaE/(KB*T(i))) underflows to 0.
beso ss
on 21 Mar 2018
exp(deltaE/(KB*T(i) put there is T change from 0 to 300 so the answer have to change in every point
Walter Roberson
on 21 Mar 2018
Unfortunately even if you switch everything to use symbolic constants for higher precision, vpasolve() has difficulty working with numbers that small. There are solutions, but they are on the order of 1E-133056 to 1E-1348
beso ss
on 21 Mar 2018
if i didn't use vpasolve() ,,, what else i can use to solve this equation ????
where T=0:300 and KB = 8.617*10^-5
in the figer (peak posission mean E(T))
Walter Roberson
on 21 Mar 2018
"exp(deltaE/(KB*T(i) put there is T change from 0 to 300 so the answer have to change in every point"
exp(deltaE/(KB*T(i)) is the same as exp((deltaE/KB)/T(i)), so you can compute deltaE/KB as a single constant, c, getting an expression of exp(c/T(i)) . With the values you have for deltaE and KB, that leads to exp(-928397.354067541/T(i)) . Over the range T = 0 to T = 300, that is exp(-928397.354067541/0) to exp(-928397.354067541/300) which is exp(-inf) to exp(-3094.6578468918) . exp(-inf) is 0 by definition. exp(-3094.6578468918) is 0 to within the available accuracy of double precision, as it is 1.016655160e-1344 . Even if you switch your computation to use symbolic numbers, the symbolic engine decides that the expression involving 1.016655160e-1344 is round off error on 0 and ends up ignoring that term, effectively leading to x*exp(x) == 0 which has a solution of x = 0.
beso ss
on 21 Mar 2018
this is the units for KB I choose the one in the middele becouse it has the same units as the other values , eV.K^-1
beso ss
on 21 Mar 2018
here how it must be look like , 1 it means E1 and 2 mean E2 you see how is the curve for E2 ,,, this is becouse of x
Walter Roberson
on 22 Mar 2018
Near T=1/100, X is approximately 1E-40319770 and by T=300 it has barely risen to 1E-1348 . That range of value has no numeric effect on E2.
beso ss
on 28 Mar 2018
if i want to solve this equation by using (bisection) how can i do it ??? also if i have a condition in x which is 0<x< (sigma/ (KB * T )^2)
Walter Roberson
on 28 Mar 2018
"if i want to solve this equation by using (bisection) how can i do it "
You cannot solve this by bisection. The range of values involved is too small for MATLAB to deal with, including in the Symbolic Toolbox. For example near T = 1/100, X is a number that has more than 40 million decimal places before the first digit.
Even if you were able to solve to that many decimal places, you would still have the fundamental problem that X is not large enough to have any visible effect on E2. If you were to magnify the plot to be the diameter of the Universe, then X would move the plot by far far less than one planck distance.
beso ss
on 28 Mar 2018
my prof told me that he want me to solve it by using (bisection) becouse it gives an accurate answer whould you help me to solve it in bisection lets try to do it with bisection but i dont know how to use bisection in the equation whould you help me please
Walter Roberson
on 28 Mar 2018
No, I will not help you solve it by bisection. It will not give you an answer that is any better than you can obtain with fsolve() [numeric] or vpasolve() [symbolic] .
John and I have independently demonstrated at length that for those equations and constants that your output cannot possibly be anything like you expect. Not for that range of T values. If you let T go over roughly 1500 then you can start to get a shape closer to what you are expecting, but with a notably different fall-off.
You should go back over the equations and the constants very carefully to see if you have made a typing mistake somewhere.
beso ss
on 28 Mar 2018
i spend with this equation more than 4 weeks i chick it every day ,, i get tired from it ,, i mager in physics not computer scince ,,
my prof solve it with bisection ,, and he said it give him an accurate answer ,,, i dont know what to do
beso ss
on 28 Mar 2018
The bisection method can be applied for solving nonlinear equations like f (x) = 0, only in the case where we know some interval [a, b] on which f (x) is continuous and the solution uniquely exist
beso ss
on 28 Mar 2018
(the best solving method is "bisection" Try to find accurate results !!)
that what he write for me ..
Walter Roberson
on 29 Mar 2018
Edited: Walter Roberson
on 29 Mar 2018
One of the following must be true:
- you made a mistake in copying one or more of the constants
- you made a mistake in forming the equation implementation leading up to the vpasolve
- the paper is wrong
- your professor is wrong
- both John and I have independently mis-analyzed it and two very different programs are making major major arithmetic errors in calculation.
beso ss
on 29 Mar 2018
this are all the equations and the pictuer contain the constant and KB= 8.617*10^-5; %eV/K
T= 0:300
the pictuer contain E2 there are 4 value of delta E which is (Ea-E0)
iam not wrong , you can make sure by your self and the prof say that he had the same picture by using matlab
Walter Roberson
on 29 Mar 2018
I am having difficulty reading some of the fine print. Which paper is this you are working on? Seems to be something about phonon scattering and microcavities. (The papers that I find with the appropriate keywords seem to be mostly behind paywalls unfortunately.)
Walter Roberson
on 29 Mar 2018
Fortunately the article is available for free via https://www.researchgate.net/publication/258276749_A_model_for_steady-state_luminescence_of_localized-state_ensemble
beso ss
on 30 Mar 2018
yesterday i was trying to use bisection but i didnt know how ??? and the prof told me he didn't want any soulution but in bisection
can any one help me to do it in matlab please
Walter Roberson
on 30 Mar 2018
"can any one help me to do it in matlab please"
No, the equations you have cannot be solved by bisection.
Suppose you were to take the Universe and magnify the smallest possible distance in it, the planck distance, to become the full size of the Universe. Like "Universe squared". Then the smallest possible distance in that magnified Universe^2 would still be larger than the largest possible result for E in the range T = 0 to T = 300. MATLAB is not able to represent numbers that small.
beso ss
on 30 Mar 2018
but the prof says that he use it and he had an accurate result .... is there any aother way to solve it rather than vpasolve..couse he says it dosent give us a physical meaning
beso ss
on 30 Mar 2018
what do you think the best equation to solve this ,,, and git the same as this figer ???? rather than vpasolve
Walter Roberson
on 1 Apr 2018
I will have to defer to your professor on this. Different specialties have different notations, that lead to different interpretations of formulas. For example in some places the [] brackets are just a different style to make it easier to find matches, but in other specialties they mean that the integer part is to be taken. I must be lacking knowledge of special notation for the topic that you are working on — because it is absolutely certain that if you use the regular interpretation of operations and if you have copied the equations and constants correctly then the solution is many orders of magnitude different from what would be produced by whatever specialized notation is being used for your question.
beso ss
on 1 Apr 2018
i understand what do you mean ,, and i know that but the problem is i don't know nothing about matlab ,, and the prof want me to do this equations without helping me ,, and he want it tomorrw to pass this semester ,,and he know that i dont know anything about matlab or programming.. every day im trying to read about matlab ,, but im still dont know how to do that equation ,, thats way iam asking for help ,,, i want to see the steps of biesction and i will try to do my equation with it ,, even if it didn't work ,, i want to show it to the prof ,,
Walter Roberson
on 1 Apr 2018
You can search MATLAB Answers for the key word bisection
I guarantee that bisection will fail for the combination of equations and constants that you have.
beso ss
on 2 May 2018
Edited: Walter Roberson
on 2 May 2018
hi it is me again hhhhh
i write the code by using biesction
and here it is ;
syms sigma deltaE E0 tautr taur KB E2
sigma = 13*10^-3; %eV
deltaE = -0.073; %eV
E0 = 1.185; %eV
tautr = 0.027;
taur =250;
KB= 8.617*10^-5; %eV/K
T=0:300;
X=zeros(1,numel(T));
E2 = zeros(size(T));
for i=1:numel(T)
syms x
min=0;
max=(sigma/(KB*T(i)))^2;
f=@(x) ((sigma/(KB*T(i)))^2-x)*(tautr/taur)*exp(deltaE/(KB*T(i)))-x*exp(x);
X(i)=bisection(f,f(max),f(min))
E2(i) = E0 -( alfa * T(i)^2 )/(theta + T(i)) - X(i) * KB * T(i);
clear x
end
subplot(1,3,1)
subplot(1,3,2)
plot(T,E2)
title('E2')
subplot(1,3,3)
plot(T,X)
title('X')
i have a problem with it ,, that x begins with minus but i wrote a condition that x bigger than zero and less than (sigma/(KB*T(i)))^2
beso ss
on 21 Mar 2018
2 Comments
Greg Heath
on 28 Mar 2018
Edited: Greg Heath
on 31 Mar 2018
CORRECTION:
God bless John and Walter!
Greg
Yi-Lin Tsai
on 4 May 2018
I have a nonlinear equation I see some references but i don't know to get each value about sigma tautr taur please help me
2 Comments
Walter Roberson
on 4 May 2018
You should start a new Question on that topic. Be sure to include a copy of the equations.
khouloud abiedh
on 3 Jul 2018
Edited: Walter Roberson
on 4 Jul 2018
hello ,i m debutant in matlab and i need to do a fit of my experimental data by LSE model in matlab iwrite the last comment i write the same code of the energetic position of a gaussian but i find a problem to do an exact fit of the data by matlab ,also i need to do a fit of the LMH of the gaussian by this model in matlab .when i read the article i understand that it must resolve another equation numerically to find the LMH described in the LSE model.i use the same code but i change the equation but it display an error.i need our help please help me and thank you in advance.
%%trouver x et tracer la position energetique
sigma = 13*10^-3;%eV
deltaE = -73e-3;%eV
E0 = 1.185;%eV
Ea = E0+75e-3;
tautr=0.027 ;
taur=250;
Kb = 8.617*10^-5;%eV/K
theta = 270; %K
alfa = 0.48e-3;%eV/K
T= 0:10:300;
X = zeros(size(T));
E2 = zeros(size(T));
E1= zeros(size(T));
a = zeros(size(T));
syms x
for i=1:numel(T)
E1(i) = E0 -( alfa*T(i)^2)/(theta+T(i));
X(i) = vpasolve(((sigma/(Kb*T(i)))^2 - x)*(taur/tautr)*exp(deltaE/(Kb*T(i)))-x*exp(x)==0,x);
E2(i) = E0 -(alfa*T(i)^2)/(theta + T(i))-X(i)*Kb*T(i);
a(i) = X(i)*Kb*T(i);
end
figure(1);plot(T,E1,'r.');
figure(2);plot(T,E2,'g')
hold on ;plot(T,E1,'r')
figure(3);plot(T,X,'k.','Markersize',5)
figure(4);plot(T,a,'k');
%%to find LMH i must resolve n(E0-X(T)*Kb*T,T)/2=n(E,T)
trouver n(E0-X(T)*Kb*T,T)/2
n(E,T)= exp(-(E(i)-E0)^2/(2*sigma^2)/(exp((E(i)-Ea)/(Kb*T(i)))+(tautr/taur)))
E = zeros(size(T));
n = zeros(size(T));
n1=zeros(size(T));
for i=1:numel(T)
E(i)= E0-a(i);
n(i)=exp(-(E(i)-E0).^2/(2*sigma^2))/(exp((E(i)-Ea)/(Kb*T(i)))+(tautr/taur));
n1(i)=n(i)/2;
end
figure (5)
plot(E1,n1)
X1=zeros(size(T));
syms x1
for i=1:numel(T)
X1(i) = vpasolve(exp(-(E(i)-E0)^2/(2*sigma^2))/(exp((E(i)-Ea)/(Kb*T(i)))+(tautr/taur))-n1(i)==0,x1);
end
figure(6);plot(T,X1,'r.')
3 Comments
khouloud abiedh
on 7 Jul 2018
hello dear Walter Roberson i start a new question in this topic i tag you but i don't know if i do a fault and i tag another person.
See Also
Categories
Find more on Nonlinear Analysis in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)