Why are there different results of 'int' function? (regarding 'int' and 'double' function

%%Version_1(without 'double' function)


syms e e_1 e_2 e_3 e_4                                                                 % e is energy


e_i = 1.2;                                  
r = 0.0036;
l = 0.78e-9;                                                               %[m]
p = 10.4*sqrt(e/e_i)*(((e-e_i)/e_i)^2);

Eox_raw_1 = [5.13411,5.36051,5.55457,5.81331,6.00737,6.20142,6.36314,6.55719,6.75125,6.91296,7.1717,7.36576,7.55982,7.72153,7.91559,8.10964,8.3037,8.46541,8.65947,8.88587,9.14461,9.46804,9.75912,9.95318,10.21192,10.40598,10.60003,10.92346,11.1822,11.44094,11.76437,12.02311,12.28186,12.5406,12.86403,13.18745,13.47854,13.70494,13.89899,14.12539,14.38413,14.61053,14.90162];
w_raw_1 = [3.57723,3.70289,3.79265,3.91832,4.02603,4.11579,4.20555,4.29531,4.38507,4.49279,4.6005,4.70821,4.81593,4.90569,4.99545,5.10316,5.19292,5.28268,5.35449,5.48016,5.60582,5.74944,5.89306,5.98282,6.09053,6.18029,6.2521,6.41367,6.52139,6.6291,6.75476,6.84453,6.93429,7.042,7.14971,7.25743,7.34719,7.419,7.4908,7.54466,7.61647,7.68828,7.76009];    

mu = 0.14;   %[m^2/vs]                                                     % mobility of electron in silicon
                                                                           
w_in = w_raw_1;

Unrecognized function or variable 'w_raw_1'.

E_ox1 = Eox_raw_1.*10^8;
E_si = E_ox1.*(3.9/11.7);


sz = size(w_in);
cnt = 0;

for i = 1:1:sz(2) 

cnt = cnt + 1;

v_g(e) = (r*p/(E_si(cnt)*l))*(e_i/e);   %
v_d(e) = (2*r*p*e_i)/((E_si(cnt)*l)^2);  %*sqrt(e*e_i)


P_bl_1(cnt) = int(exp(-int(v_g(e_1), e_1, 0, e))*exp(-int(v_g(e_2), e_2, e_i, e)) ...
    *(v_g(e)) ,e , e_i, w_in(cnt));

P_dr_1(cnt) = int(int(exp(-int(v_g(e_1), e_1, 0, e_2))*(v_g(e_2)) ...
    *exp(-int(v_d(e_1), e_1, e_2, e)) , e_2, 0, e)*exp(-int(v_d(e_4), e_4, e_i, e)) ...
    *(v_d(e)),e , e_i, w_in(cnt));

P_II_1(cnt) = P_bl_1(cnt) + P_dr_1(cnt);

end


%%Version_2(with 'double' function)

syms e e_1 e_2 e_3 e_4                                                                 % e is energy


e_i = 1.2;                                  
r = 0.0036;
l = 0.78e-9;                                                               %[m]
p = 10.4*sqrt(e/e_i)*(((e-e_i)/e_i)^2);


mu = 0.14;   %[m^2/vs]                                                     % mobility of electron in silicon
                                                                           
w_in = w_raw_1;
E_ox1 = Eox_raw_1.*10^8;
E_si = E_ox1.*(3.9/11.7);


sz = size(w_in);
cnt = 0;

for i = 1:1:sz(2) 

cnt = cnt + 1;

v_g(e) = (r*p/(E_si(cnt)*l))*(e_i/e);   %
v_d(e) = (2*r*p*e_i)/((E_si(cnt)*l)^2);  %*sqrt(e*e_i)


P_bl_1(cnt) = double(int(exp(-int(v_g(e_1), e_1, 0, e))*exp(-int(v_g(e_2), e_2, e_i, e)) ...
    *(v_g(e)) ,e , e_i, w_in(cnt)));

P_dr_1(cnt) = double(int(int(exp(-int(v_g(e_1), e_1, 0, e_2))*(v_g(e_2)) ...
    *exp(-int(v_d(e_1), e_1, e_2, e)) , e_2, 0, e)*exp(-int(v_d(e_4), e_4, e_i, e)) ...
    *(v_d(e)),e , e_i, w_in(cnt)));

P_II_1(cnt) = P_bl_1(cnt) + P_dr_1(cnt);

end

The difference between 'version_1' and 'Version_2' is only whether the double function is used in the integration result(At 'P_bl_1' and 'P_dr_1' lines). But Matlab returns different results.

This is result of version_1(without 'double'). Matlab returns an uncalculated integral.(above figure)

But, Matlab returns a calculated integral at Version_2(with 'double'). (below figure)

I don't know why Matlab returns different results. Do you know exactly why this difference in results occurs?

My ultimate goal is to obtain the symbolic integration result of below equation and put it to the 'origin' for fitting.

So, I want to understand why this difference occurs in the first place.

Please, Help me...

exp(-int(v_g(e_1), e_1, 0, e))*exp(-int(v_g(e_2), e_2, e_i, e))*(v_g(e))

int(exp(-int(v_g(e_1), e_1, 0, e_2))*(v_g(e_2)) ...
    *exp(-int(v_d(e_1), e_1, e_2, e)), e_2, 0, e)*exp(-int(v_d(e_4), e_4, e_i, e))*(v_d(e))

[Ayush Modi]

Can you provide the data ?

It will be easier to check and help you out.

[han-byeul Oh]

Sure

I think all that is needed for considering the above script are these.

w_in(1) = 3.577

E_si1(1) = 1.71e+8

are there any other additional data required?

[Walter Roberson]

In the current posted version, w_raw_1 is not defined.

[han-byeul Oh]

w_raw_1 is 1x47 double data. the range of w_raw_1 is between 3.577 and 7.76.

w_raw_1 = [3.57723,3.70289,3.79265,3.91832,4.02603,4.11579,4.20555,4.29531,4.38507,4.49279,4.6005,4.70821,4.81593,4.90569,4.99545,5.10316,5.19292,5.28268,5.35449,5.48016,5.60582,5.74944,5.89306,5.98282,6.09053,6.18029,6.2521,6.41367,6.52139,6.6291,6.75476,6.84453,6.93429,7.042,7.14971,7.25743,7.34719,7.419,7.4908,7.54466,7.61647,7.68828,7.76009];    

these are values extracted from the paper.

[Walter Roberson]

w_raw_1 = [3.57723,3.70289,3.79265,3.91832,4.02603,4.11579,4.20555,4.29531,4.38507,4.49279,4.6005,4.70821,4.81593,4.90569,4.99545,5.10316,5.19292,5.28268,5.35449,5.48016,5.60582,5.74944,5.89306,5.98282,6.09053,6.18029,6.2521,6.41367,6.52139,6.6291,6.75476,6.84453,6.93429,7.042,7.14971,7.25743,7.34719,7.419,7.4908,7.54466,7.61647,7.68828,7.76009];    

%%Version_1(without 'double' function)


syms e e_1 e_2 e_3 e_4                                                                 % e is energy


e_i = 1.2;                                  
r = 0.0036;
l = 0.78e-9;                                                               %[m]
p = 10.4*sqrt(e/e_i)*(((e-e_i)/e_i)^2);


mu = 0.14;   %[m^2/vs]                                                     % mobility of electron in silicon
                                                                           
w_in = w_raw_1;
E_ox1 = Eox_raw_1.*10^8;

Unrecognized function or variable 'Eox_raw_1'.

E_si = E_ox1.*(3.9/11.7);


sz = size(w_in);
cnt = 0;

for i = 1:1:sz(2) 

cnt = cnt + 1;

v_g(e) = (r*p/(E_si(cnt)*l))*(e_i/e);   %
v_d(e) = (2*r*p*e_i)/((E_si(cnt)*l)^2);  %*sqrt(e*e_i)


P_bl_1(cnt) = int(exp(-int(v_g(e_1), e_1, 0, e))*exp(-int(v_g(e_2), e_2, e_i, e)) ...
    *(v_g(e)) ,e , e_i, w_in(cnt));

P_dr_1(cnt) = int(int(exp(-int(v_g(e_1), e_1, 0, e_2))*(v_g(e_2)) ...
    *exp(-int(v_d(e_1), e_1, e_2, e)) , e_2, 0, e)*exp(-int(v_d(e_4), e_4, e_i, e)) ...
    *(v_d(e)),e , e_i, w_in(cnt));

P_II_1(cnt) = P_bl_1(cnt) + P_dr_1(cnt);

end


%%Version_2(with 'double' function)

syms e e_1 e_2 e_3 e_4                                                                 % e is energy


e_i = 1.2;                                  
r = 0.0036;
l = 0.78e-9;                                                               %[m]
p = 10.4*sqrt(e/e_i)*(((e-e_i)/e_i)^2);


mu = 0.14;   %[m^2/vs]                                                     % mobility of electron in silicon
                                                                           
w_in = w_raw_1;
E_ox1 = Eox_raw_1.*10^8;
E_si = E_ox1.*(3.9/11.7);


sz = size(w_in);
cnt = 0;

for i = 1:1:sz(2) 

cnt = cnt + 1;

v_g(e) = (r*p/(E_si(cnt)*l))*(e_i/e);   %
v_d(e) = (2*r*p*e_i)/((E_si(cnt)*l)^2);  %*sqrt(e*e_i)


P_bl_1(cnt) = double(int(exp(-int(v_g(e_1), e_1, 0, e))*exp(-int(v_g(e_2), e_2, e_i, e)) ...
    *(v_g(e)) ,e , e_i, w_in(cnt)));

P_dr_1(cnt) = double(int(int(exp(-int(v_g(e_1), e_1, 0, e_2))*(v_g(e_2)) ...
    *exp(-int(v_d(e_1), e_1, e_2, e)) , e_2, 0, e)*exp(-int(v_d(e_4), e_4, e_i, e)) ...
    *(v_d(e)),e , e_i, w_in(cnt)));

P_II_1(cnt) = P_bl_1(cnt) + P_dr_1(cnt);

end

Answer 1

[Walter Roberson]

It looks okay to me.

int() is not a numeric integral: it is a request to calculate the integral to full theoretical precision.

When you request to calculate an integral to full theoretical precision using int(), there are several possibilities:

1. The calculation involves symbolic floating point numbers. In this case, MATLAB attempts to calculate using vpaintegral()
2. MATLAB might know how to calculate it in full form, and does so, returning the exact theoretical form.a formula is returned
3. sometimes MATLAB is able to determine that the integral is undefined, and so returns NaN
4. The exact "closed form" integral might exist but it might be too complicated or beyond MATLAB's ability, and either way MATLAB returns the integral un-evaluated
5. There might be no known "closed form" integral; it might even be known that there is no closed form for the integral. This is in effect the same as the previous: no matter what the reason, MATLAB does not know a closed form solution for the integral and so returns the integral un-evaluated

The only case in which int() will return a floating point number is the one in which the formula involved a symbolic floating point number. Normal floating point numbers are automatically converted to rationals or square roots of rationals or to π when used in symbolic expressions, so in order to get a symbolic floating point number, the user either had to construct the number to be deliberately symbolic floating point, such as sym('0.1') -- or else vpa() had to have been used at some stage.

You did not use any symbolic floating point numbers, so you are left with (4) or (5) -- that MATLAB did not know how to calculate the exact integral and so returned the int() unevaluated.

Now, when you have an unevaluated symbolic formula with int() in it, there are a few different possibilities:

1. the symbolic formula might have a "free variable" in it -- a symbolic variable that is not being used as the dummy variable for an int() or symsum() or symprod() . In such a case, double() will not be able to process the expression and an error will be given. This applies even if you could theoretically zero out the coefficients of the variable, such as a*(b - 1)*(b + 1) - a*(b^2 - 1)
2. the int() might not have definite bounds; in this case the integral would be a formula in the variable of integration, and the variable of integration would be a free variable and double() would not be able to process that
3. the int() might have an exact symbolic numeric solution. However, if that were the case, the int() would already have created the exact symbolic numeric solution so this would violate our assumption that you have an unevaluated symbolic formula with int() in it.
4. a numeric integration will be attempted; if the integral does not converge to the default precision levels, then the numeric integral failed and the numeric integral would internally return an unevaluated numeric integration and double() would not be able to process that so an error would be given
5. numeric integration might find the integral is infinite or undefined; double() is able to process those results to return inf or nan.
6. numeric integration might succeed and will internally return a symbolic floating point number that double() will then be able to process, returning a double precision number.

You are encountering the last of those situations: int() was not able to create an exact solution, but symbolic numeric integration worked just fine, and double() converted that result to a double precision number.

There is no conflict here.

[han-byeul Oh]

Ultimately, MATLAB is not able to create an exact symbolic results, but numeric results for above equation. Right?. Then... there may be no solution for getting symbolic results..?

[Walter Roberson]

I need Eox_raw_1 and any other undefined numeric variable to test with.

[han-byeul Oh]

Sorry I'm late.

Eox_raw_1 = [5.13411,5.36051,5.55457,5.81331,6.00737,6.20142,6.36314,6.55719,6.75125,6.91296,7.1717,7.36576,7.55982,7.72153,7.91559,8.10964,8.3037,8.46541,8.65947,8.88587,9.14461,9.46804,9.75912,9.95318,10.21192,10.40598,10.60003,10.92346,11.1822,11.44094,11.76437,12.02311,12.28186,12.5406,12.86403,13.18745,13.47854,13.70494,13.89899,14.12539,14.38413,14.61053,14.90162];

these data are values extracted from the paper, too.

Thank you

[Walter Roberson]

In the below, if you use Maple to rewrite temp1 in terms of sinh cosh, then the result is something that Maple can integrate to closed form involving terms of exp()

Unfortunately, Maple is not able to find a closed form for the inner int() of the temp2

At the moment Maple is telling me that the second integral has a problem with division by 0 inside a natural log... I do not know at the moment what is causing that.

Q = @(v) sym(v);
Eox_raw_1 = [5.13411,5.36051,5.55457,5.81331,6.00737,6.20142,6.36314,6.55719,6.75125,6.91296,7.1717,7.36576,7.55982,7.72153,7.91559,8.10964,8.3037,8.46541,8.65947,8.88587,9.14461,9.46804,9.75912,9.95318,10.21192,10.40598,10.60003,10.92346,11.1822,11.44094,11.76437,12.02311,12.28186,12.5406,12.86403,13.18745,13.47854,13.70494,13.89899,14.12539,14.38413,14.61053,14.90162];
w_raw_1 = [3.57723,3.70289,3.79265,3.91832,4.02603,4.11579,4.20555,4.29531,4.38507,4.49279,4.6005,4.70821,4.81593,4.90569,4.99545,5.10316,5.19292,5.28268,5.35449,5.48016,5.60582,5.74944,5.89306,5.98282,6.09053,6.18029,6.2521,6.41367,6.52139,6.6291,6.75476,6.84453,6.93429,7.042,7.14971,7.25743,7.34719,7.419,7.4908,7.54466,7.61647,7.68828,7.76009];    

%%Version_1(without 'double' function)


syms e e_1 e_2 e_3 e_4                                                                 % e is energy


e_i = Q(1.2);                                  
r = Q(0.0036);
l = Q(0.78)*Q(10)^-9;                                                               %[m]
p = Q(10.4)*sqrt(e/e_i)*(((e-e_i)/e_i)^2);


mu = Q(0.14);   %[m^2/vs]                                                     % mobility of electron in silicon
                                                                           
w_in = Q(w_raw_1);
E_ox1 = Q(Eox_raw_1).*Q(10)^8;
E_si = E_ox1.*(Q(3.9)/Q(11.7));


sz = size(w_in);
cnt = 0;

for i = 1:1:1 %sz(2) 

cnt = cnt + 1;

v_g(e) = (r*p/(E_si(cnt)*l))*(e_i/e);   %
v_d(e) = (2*r*p*e_i)/((E_si(cnt)*l)^2);  %*sqrt(e*e_i)


temp1 = int(exp(-int(v_g(e_1), e_1, 0, e))*exp(-int(v_g(e_2), e_2, e_i, e)) ...
    *(v_g(e)) ,e , e_i, w_in(cnt))
P_bl_1(cnt) = rewrite(temp1, 'sinhcosh');

temp2 = int(int(exp(-int(v_g(e_1), e_1, 0, e_2))*(v_g(e_2)) ...
    *exp(-int(v_d(e_1), e_1, e_2, e)) , e_2, 0, e)*exp(-int(v_d(e_4), e_4, e_i, e)) ...
    *(v_d(e)),e , e_i, w_in(cnt));

P_dr_1(cnt) = temp2; %rewrite(temp2, 'sinhcosh');

P_II_1(cnt) = P_bl_1(cnt) + P_dr_1(cnt);

end

temp1 = 

P_bl_1(1)

ans = 

P_bl_1(1)

ans = 

char(P_bl_1(1))

ans = 'int((30399297484750848*((5*e)/6 - 1)^2*(cosh((281474976710656*30^(1/2)*e^(1/2)*(5*e^2 - 20*e + 36))/90320218292496625) - sinh((281474976710656*30^(1/2)*e^(1/2)*(5*e^2 - 20*e + 36))/90320218292496625))*(cosh((27021597764222976*5^(1/2)*6^(1/2)*30^(1/2))/2258005457312415625 - (281474976710656*30^(1/2)*e^(1/2)*(5*e^2 - 20*e + 36))/90320218292496625) + sinh((27021597764222976*5^(1/2)*6^(1/2)*30^(1/2))/2258005457312415625 - (281474976710656*30^(1/2)*e^(1/2)*(5*e^2 - 20*e + 36))/90320218292496625))*((5*e)/6)^(1/2))/(90320218292496625*e), e, 6/5, 251725182734665/70368744177664)'

char(P_dr_1(1))

ans = 'int((34226566206162193840410986545152*exp((182541686432865033815525261574144*5^(1/2)*6^(1/2)*30^(1/2))/5938836053990288142572260663292375 - (316912650057057350374175801344*30^(1/2)*e^(3/2)*(25*e^2 - 84*e + 84))/47510688431922305140578085306339)*((5*e)/6 - 1)^2*int((30399297484750848*exp((316912650057057350374175801344*30^(1/2)*e_2^(3/2)*(25*e_2^2 - 84*e_2 + 84))/47510688431922305140578085306339 - (316912650057057350374175801344*30^(1/2)*e^(3/2)*(25*e^2 - 84*e + 84))/47510688431922305140578085306339)*exp(-(281474976710656*30^(1/2)*e_2^(1/2)*(5*e_2^2 - 20*e_2 + 36))/90320218292496625)*((5*e_2)/6 - 1)^2*((5*e_2)/6)^(1/2))/(90320218292496625*e_2), e_2, 0, e)*((5*e)/6)^(1/2))/6787241204560329305796869329477, e, 6/5, 251725182734665/70368744177664)'