Tutorial 4 / CW
on Numerical Accuracy and Precision

Treat part of this as lecture materials. You should study and practice at home.

1. ListPlot/LineLinePlot practice

Plot Cos[x] for x from 0 to 50 by making a ~500-point array. Do with an array of point instead of built-in function.

Plot[Cos[x],{x,0,50}]

Answer

xarr1=Range[0,499]*50./500;
cosxarr1=Cos[xarr1];
ListLinePlot[Transpose[{xarr1,cosxarr1}]
,Filling-> {1-> {0,Hue[0.45,1,1]}}]

2. Problem 2 - EM wave long distance

Light is electromagnetic wave, and for some, such as a laser beam, we can describe the E field as:
                                 Ef=A (e^(i (k x-ω t)))              (3.1)
where k=((2 π)/λ)=(ω/c) where λ is the wavelength of light, ω is the frequency, and c is the speed of light which is:     c=299.792458((μm)/ps)
Note: μm is micrometer = (10^(-6)) meter
         ps is picosecond = (10^(-12)) second.

Consider a laser beam shot by NASA Lunar Atmosphere and Dust Environment Explorer (LADEE) from Moon to Earth. Although its wavelength is infrared (for minimal atmospheric scattering loss), for the sake of this HW, we will consider laser wavelength be 0.5 μm (green).

If we plot the E field at a distance d, with x over a range of 5 μm at 0.02 μm step at a moment in time, we can write as follow:
                       Ef=A (e^(i (k (d+x)-ω t))) -> A Cos[k (d+x) ]   for x =0 to 5 μm.
                       λ=0.5 μm;  k=((2π)/λ); and Let A=1.

2.1 Plot E-field at source

Let d=0, plot E field.

Answer

x range is from d+0 to d+50 (μm). Efield is the cos of that.
                   Ef=A (e^(i (k x-ω t)))      k=((2 π)/λ)

xarr2=Range[0,5,0.02];
λ=0.5 ; k=2 π/λ ;
d=0.; 
ef1=Cos[k (d+xarr2) ];
ListLinePlot[Transpose[{xarr2,ef1}],Filling-> {1-> {0,Hue[0.45,1,1]}}]

2.2 Plot E-field at d=1000 km

Plot the field after having traveled a distance d=1000 km away from its starting position.
Note: to convert km to μm:  1 km=(10^9) μm.

Answer

xarr2=Range[0,5,0.02];
λ=0.5 ; k=2 π/λ ;
d=1000*(10^9); 
ef2=Cos[k (d+xarr2) ];
ListLinePlot[Transpose[{xarr2,ef2}],Filling-> {1-> {0,Hue[0.45,1,1]}}]

2.3 Plot the E field when it reaches Earth at d = 384000 km.

Plot the E field when it reaches Earth at d = 384000 km. (384000*(10^9) μm)

Answer

xarr2=Range[0,5,0.02];
λ=0.5 ; k=2 π/λ ;
d=384000 (10^9); 
ef3=Cos[k (d+xarr2) ];
ListLinePlot[Transpose[{xarr2,ef3}],Filling-> {1-> {0,Hue[0.45,1,1]}}]

2.4 Discussion

Is the above correct result? i. e. is that the electromagnetic field when it reaches Earth? Discuss.

Answer

No. Electric field as expressed by      Ef=A (e^(i (k x-ω t)))  remains the same regardless of distance: periodic vs. λ=((2 π)/k). The graph is caused by a loss of accuracy and precision in the computation.

2.5 Is there any difference between E field after 1000 km vs E field at the starting position? They “look” the same

Is there any difference between ef2 (E field after 1000 km) with ef1 (E field at the starting position)? 
Hint: plot xarr2 vs. ef2-ef1

Answer

“Look” can be deceiving. We must be aware of the error long before it becomes obvious to the eye.

ListLinePlot[Transpose[{xarr2,ef1-ef2}]]

We see that the error is already there. quite big, ~ (10^(-3))

ListLinePlot[Transpose[{xarr2,ef1-ef3}]]

It is the same error, just that it becomes so big that it is now obvious.

Understanding (operational knowledge) of computer limits means to be aware of errors long before they become disastrous.

2.6 What is the “Numerical-Methods”-wise correct approach?

Answer

Don’t force a function out of its computationally valid range:
      (e^x), or Sin, Cos ((e^(i x)),  (e^(-i x))): all of them, as performed by the FPU of your CPU, have errors depending on the input variables.
Here, the correct solution is to use the periodicity property:

xarr2=Range[0,5,0.02];
λ=0.5 ; k=2 π/λ ;
d=384000 (10^9); d0=Mod[d,λ];
efcorrect=Cos[k (d0+xarr2) ];
ListLinePlot[Transpose[{xarr2,efcorrect}],Filling-> {1-> {0,Hue[0.45,1,1]}}]

Lecture break: mantissa, exponent, precision, accuracy

Goto Course Notes webpage for Lecture 2. Link here.

3. Problem 3 - numbers inside a computer

3.1 (Subscript[y, 1])= (Sqrt[(Subscript[x, 1])]), find (Subscript[x, 1])-((Subscript[y, 1])^2)

Let x1 be some number,  y1=(Sqrt[x1]), find z = x1-(y1^2).

Answer

3.2 Do 3.1 above with (Subscript[x, 1])=5.

Now, let x1=5., use computer to calculate z

x1=5. ; y1=(Sqrt[x1]);

Answer

x1=5. ; y1=(Sqrt[x1]);
z=x1-(y1^2)

3.3 Take (Subscript[Log, 2]) of z

Answer

Log[2,Abs[z]]

If (Subscript[log, 2])[z]= -50, it means z=(2^(-50)) =(1/(2^50))  . Let’s verified:

(1./(2^50))

why? we’ll see that in the lecture about binary computer architecture: Lecture Set 3. Link here.

3.4 Make a Table of many numbers

Make a table (array) of {x1 (integer), z=Abs[ x1-((((Sqrt[1.*x1])))^2)],  and Log[2,z]} for x1 from 0 to 50.

Answer

u1=Table[ {x1,z=Abs[x1-((((Sqrt[1.*x1])))^2) ],Log[2,z]}  ,{x1,0,50}] ;
Grid[u1, Frame-> All]

3.5 Obtain array {x1,z}, name it u2, and ListPlot x1 vs. z

u2=u1[[ ;; ,1 ;; 2]];
ListPlot[u2, PlotRange-> All]

3.6 Do the same as 3.5 but x1 from 0 to 1000000 in step of 2503 and use ListLogLogPlot

step=NextPrime[2500]
u3=Table[ {x1,z=Abs[x1-((((Sqrt[1.*x1])))^2) ]}  ,{x1,0,1000000, step}] ;
u3graph1=ListLogLogPlot[u3, Frame-> True, GridLines->Automatic, AspectRatio->1]

Notice something: the error grows proportionally with the number magnitude: the error is in the mantissa, and it grows proportionately with the exponent.

Extra: For those who want to see an example of linear regression

step=NextPrime[2500] ;
u3=Table[ {x1,z=Abs[x1-((((Sqrt[1.*x1])))^2) ]}  ,{x1,0.,1000000, step}] ;
(* just the same as above, but x1 is now floating point, and we take the log of it  *)
u3Log=Select[Log[u3],NumericQ[#[[1]]] &&  NumericQ[#[[2]]] &] ;
u3linearregression=LinearModelFit[u3Log,{1,x},x]
xlogMinMax=MinMax[Transpose[u3Log][[1]]];
u3LinRegFitplot=Plot[u3linearregression[x],{x,xlogMinMax[[1]],xlogMinMax[[2]]}
,PlotRange-> All, PlotStyle->{Thick,Red}]

Show[u3graph1,u3LinRegFitplot]

The power coefficient is ~1. It means the error grows linearly with respect to the number magnitude. The reason is that the number is stored as two parts: mantissa and component. Example:

Doesn’t matter how big a number is, the main change is in the 11-bit exponent, the mantissa precision error is almost the same (multiplied by the exponent).
Rule of thumb: Precision error is in the 52-bit mantissa. Magnitude (out-of-range) errors such as overflow and underflow are in the 11-bit exponent.

Lecture break: Binary, 64-bit float number - intro

Goto Course Notes webpage for Lecture 3: Binary Computer Architecture: Link here.

4. Revisit Problem 2 - EM wave

2.6 What is the “Numerical-Methods”-wise correct approach?

Answer

Don’t force a function out of its computationally correct range:
      (e^x), or Sin, Cos ((e^(i x)),  (e^(-i x)))
Use periodicity:

xarr2=Range[0,5,0.02];
λ=0.5 ; k=2 π/λ ;
d=384000 (10^9); d0=Mod[d,λ];

why can’t we use Mod for the sum of d+xarr2? What if we do?

efcorrect2=Cos[k Mod[(d+xarr2),λ] ];
ef3=Cos[k Mod[xarr2,λ] ];
ListLinePlot[Transpose[{xarr2,efcorrect2}],Filling-> {1-> {0,Hue[0.45,1,1]}}
,PlotLabel->Style["Error due to loss of precision: d is too big
to be added to a small-value array like xarr2",18,Red,Bold,FontFamily->"Calibri"]]
ListLinePlot[Transpose[{xarr2,ef3}],Filling-> {1-> {0,Hue[0.45,1,1]}}
,PlotLabel->Style["Everything is fine if we just use xarr2
(no big number added)",18,Blue,Bold,FontFamily->"Calibri"]]

This is what happens if we add a big number (d) to a much smaller (relatively) array:

xarr3=d+xarr2;
ListLinePlot[{xarr2,xarr3-d},PlotRange->{{0,20},{0,0.4}}
,PlotStyle->{{Thick,Blue},{Thick,Red}}]

Go to lecture and app for more details:

Lecture break: Floating-point - IEEE-754- detail

Use this app to get a feeling:

Home assignment (possibly) discussion

Consider this RC circuit:

Let the input be a square wave:

Plot[UnitBox[2u-0.5] /. u-> Mod[t,1],{t,0,4}
,PlotStyle-> {Thickness[0.005],Hue[0.1,1,1]}
,Filling-> {1-> {0.,Hue[0.16,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White
,FrameLabel->(Style[#,16,White,FontFamily->"Arial"] & /@ 
{"time","amplitude"})
,Background->Black,AspectRatio->0.5,ImageSize->400]

The output is: 
       (Subscript[v, out])=Π[2 t-0.5] -(e^(-(t/τ)))((1/(1+(e^(-(0.5/τ))))) -(e^(0.5/τ)) Θ[t-0.5])  for  0≤ t < 1  (2.1)
and is periodic with period of 1; τ is the RC time constant:
                        τ=((((Subscript[R, 1])+(Subscript[R, Th])) (Subscript[R, 2]))/(((Subscript[R, 1])+(Subscript[R, Th])) +  (Subscript[R, 2])))C 
For the notation of Eq. (2.1) above:
                  Π[u] is the UnitBox function
                  Θ[u] is the Heaviside UnitStep function

Below is an example of the output for τ=0.1

code

τ=0.1;Plot[(UnitBox[2u-0.5]- (E^(-(u/τ)))((1/(1+(E^(-(0.5/τ)))))-(E^(0.5/τ)) UnitStep[u-0.5]) 
/. u-> Mod[t,1]),{t,0,4} 
,PlotRange-> {0,1.1}
,PlotStyle-> {Thickness[0.005],Hue[0.35,1,1]}
,Filling-> {1-> {0.,Hue[0.3,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White
,FrameLabel->(Style[#,16,White,FontFamily->"Arial"] & /@ 
{"time","amplitude"})
,Background->Black,AspectRatio->0.5,ImageSize->400]

1 Plot the output for τ=0.01

Answer

code

We set τ=0.01;

τ=0.01;Plot[(UnitBox[2u-0.5]- (E^(-(u/τ)))((1/(1+(E^(-(0.5/τ)))))-(E^(0.5/τ)) UnitStep[u-0.5]) 
/. u-> Mod[t,1]),{t,0,4} 
,PlotRange-> {0,1.1}
,PlotStyle-> {Thickness[0.005],Hue[0.35,1,1]}
,Filling-> {1-> {0.,Hue[0.3,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White
,FrameLabel->(Style[#,16,White,FontFamily->"Arial"] & /@ 
{"time","amplitude"})
,Background->Black,AspectRatio->0.5,ImageSize->400]

2 Plot the output for τ=(10^(-4)), for t from 0.5 to 0.5005

You must explain your approach how to get the solution.

Answer

See mini homework 2