How to calculate exp() of a small number:

%% Absorbcja Prawo Lamberta i Transmisja
%format shortE
%format compact
%digits(100);
%% constance
um = 10^(-6); %micro

%% Variables
R = 0.5; %reflactance
alpha = 0.1;
L= linspace(0,10*um) % <- PROBLEM [m] lenght of the absorber
I0=1;

%% calculation
IT = I0 *((1-R)^2)*exp(-alpha.*L);
%IT = (exp(-L));
%% ploting
plot(L,IT);

The problem is how matlab rounded the samll value:

For L= linspace(0,10*um):

but the correct answer for L= linspace(0,10), is part of exp function:

So what should I do to get proper calculation?

Answer 1

[Walter Roberson]

For values in the range of 1e-7 you can use the Taylor series expansion .

1 + x + x^2/2! + x^3/3! ....

now 1e-7^3 is about 1e-21 which is less than eps(1) so you can stop before then .

Answer 2

[John D'Errico]

This is not at all surprising. For VERY small values of L, over a tiny interval, exp is extremely well represented by a straight line. That is generally true of any well behaved function. So exp actually did a very good job computing the result you wanted. And, yes, of course you see something happen when L varies from 0 to 10. Now the function shows curvature, completely as expected.

So exp did a VERY good job in evaluating the exponential. Even had you more accuracy in exp, you simply can do little more, in terms of double precision.

You can in fact get higher accuracy, if you are willing to work in symbolic or my own HPF. So, just to see if it would help, and if the computations are correct as they are seen, lets work in 50 digits of precision, using my HPF toolbox. First, using doubles as a reference:

um = 1e-6; %micro
%% Variables
R = 0.5; %reflactance
alpha = 0.1;
L = linspace(0,10*um,21)'; % 21 steps
I0=1;
%% calculation
ITdouble = I0 *((1-R)^2)*exp(-alpha.*L);

So what did doubles give?

ITdouble
ITdouble =
                      0.25
              0.2499999875
         0.249999975000001
         0.249999962500003
         0.249999950000005
         0.249999937500008
         0.249999925000011
         0.249999912500015
          0.24999990000002
         0.249999887500025
         0.249999875000031
         0.249999862500038
         0.249999850000045
         0.249999837500053
         0.249999825000061
          0.24999981250007
          0.24999980000008
          0.24999978750009
         0.249999775000101
         0.249999762500113
         0.249999750000125

That seems eminently reasonable. How about if we work in 50 digits of precision?

DefaultNumberOfDigits 50
um = hpf('1e-6'); %micro
%% Variables
R = hpf(0.5); %reflactance
alpha = hpf('0.1');
L = (0:20)*10*um/20; % Ensures the sequence is computed to full precision
I0=1;
%% calculation
IT = I0 *((1-R)^2)*exp(-alpha.*L);

So in high precision, we got what looks to be darn close.

IT
IT =
HPF array of size: 1  21
|1,1| 0.25
|1,2| 0.24999998750000031249999479166673177083268229167209
|1,3| 0.24999997500000124999995833333437499997916666701389
|1,4| 0.24999996250000281249985937500527343734179687895508
|1,5| 0.24999995000000499999966666668333333266666668888889
|1,6| 0.24999993750000781249934895837402343546549487643771
|1,7| 0.24999992500001124999887500008437499493750025312499
|1,8| 0.24999991250001531249821354182298175989127667995331
|1,9| 0.24999990000001999999733333359999997866666808888881
|1,10| 0.24999988750002531249620312542714839905664350825177
|1,11| 0.24999987500003124999479166731770826822917209201350
|1,12| 0.24999986250003781249306770928652333264909815299946
|1,13| 0.24999985000004499999100000134999983800001619999861
|1,14| 0.24999983750005281248855729352610652910614597877904
|1,15| 0.24999982500006124998570833583437464985420751700980
|1,16| 0.24999981250007031248242187829589794311529617308908
|1,17| 0.24999980000007999997866667093333265066675768887849
|1,18| 0.24999978750009031247441146377089751311406324634369
|1,19| 0.24999977500010124996962500683437376981268452810127
|1,20| 0.24999976250011281246427605015110515878996878187183
|1,21| 0.24999975000012499995833334374999791666701388883929

double(IT' - ITdouble)
ans =
                         0
     -1.32685728855204e-17
     -1.21598405648759e-17
       3.3261657119357e-18
      5.43383907928946e-18
     -5.83685171280916e-18
     -2.73036229872421e-18
      -1.3002299544076e-17
     -8.89711908322552e-18
      9.58514783383902e-18
      -1.3066681274127e-17
      6.41408918977814e-18
      1.25162767443129e-17
      5.23985013949525e-18
      1.23403537409737e-17
      6.06218068314036e-18
     -1.35947002839821e-17
      8.88083082088846e-18
     -9.77798409910904e-18
      1.36955505529395e-17
     -3.96532331982364e-18

I think for some reason, you were expecting somethign else. Perhaps, you thought that if you worked over a large range, that you see curvature, so if you narrow the range, you should still see something similar. exp did VERY well in fact, never needing to go to high precision. That curvature goesaway over a small interval.

If you really need more to understand what happens, look at that Taylor series.

exp(x) ~ 1 + x + x^2/2 + x^3/6 + x^4/24 + ...

Now, think about what happens when x is on the order of 1e-7? What is x^2/2? It is 0.5e-14. So for VERY small x, a very good approximation to exp(x) is just 1+x, accurate to nearly full double precision accuracy. While that second order term is still just barely significant at 1e-7, you could not see the difference from a perfectly straight line in any plot.