I have a nonlinear equation ,, which i didn't know how to put it in matlab

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

[John D'Errico]

First, you need to explain what you need out of it.

[beso ss]

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

is it clear now href = ""</a>>>/?????

Answer 1

[John D'Errico]

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.

[beso ss]

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 ?????????

[Torsten]

Do you have the symbolic toolbox installed ?

Test it by typing "ver" at the command line.

Best wishes

Torsten.

[beso ss]

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]

As the symbolic toolbox is not listed as licensed product, you'll have to buy and/or install it.

[beso ss]

how can i install it >>??? :(

[Torsten]

Download and read the installation instructions supplied from Mathworks.

[beso ss]

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]

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]

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]

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]

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]

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?

[beso ss]

yes i know that the value of E1 it will be so tiny ,, but not equal to 1.185 it will be smaller than it ,,, it will be around 1 and 1.185 , and also ,,E1 and E2 don't have to be the same result

[John D'Errico]

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.

Answer 2

[beso ss]

can anyone answer me please

[John D'Errico]

I did answer you. Did you think about what I said?

[beso ss]

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

Answer 3

[beso ss]

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

[beso ss]

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]

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]

but when i plot it ,, there is somthing wrong i don't know where is the wrong in my solution

[Walter Roberson]

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]

[beso ss]

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]

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]

when i put T=3000 it works with x correct , but it didn't work when i plot E2 ???

:(

[Walter Roberson]

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]

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]

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.

[beso ss]

but the T must be from 0 to 300

[Walter Roberson]

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]

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

[beso ss]

i don't know way here x=0

[Walter Roberson]

deltaE/KB is -928397.354067541 so exp(deltaE/(KB*T(i))) underflows to 0.

[beso ss]

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]

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]

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]

"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.

[Walter Roberson]

My suspicion is that you have the wrong unit on k_B

[beso ss]

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]

i dont'n know what to do to solve this equation :(( :(( :((

[Walter Roberson]

Would solutions on the order of 1E-133056 to 1E-1348 be acceptable ?

[beso ss]

it should be no minus in the answer of E1 or E2 as will

[beso ss]

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]

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]

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]

"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]

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]

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]

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]

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]

(the best solving method is "bisection" Try to find accurate results !!)

that what he write for me ..

[Walter Roberson]

One of the following must be true:

1. you made a mistake in copying one or more of the constants
2. you made a mistake in forming the equation implementation leading up to the vpasolve
3. the paper is wrong
4. your professor is wrong
5. both John and I have independently mis-analyzed it and two very different programs are making major major arithmetic errors in calculation.

[beso ss]

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

[beso ss]

:(

:(

[Walter Roberson]

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

[beso ss]

in the figer peak position means E2

the paper is ( A model for steady-state luminescence of localized-state ensemble)

[Walter Roberson]

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]

[beso ss]

:(( :(

[beso ss]

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]

"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]

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]

what do you think the best equation to solve this ,,, and git the same as this figer ???? rather than vpasolve

[beso ss]

this is condtion for x how can i write it in the code ????

[Walter Roberson]

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]

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]

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]

thank you :(

[beso ss]

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]

can anyone help me please

[beso ss]

:( anyone????

Answer 4

[beso ss]

[Greg Heath]

CORRECTION:

 God bless John and Walter!

 Greg

[beso ss]

:)

Answer 5

[Yi-Lin Tsai]

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

[Walter Roberson]

You should start a new Question on that topic. Be sure to include a copy of the equations.

[beso ss]

ok i will ,i am going to write a new question.

 thanks

Answer 6

[khouloud abiedh]

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.')

[Walter Roberson]

Please start a new question on this topic

[khouloud abiedh]

ok

[khouloud abiedh]

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.