Scientific Computing with Matlab and Octave by Quarteroni, Saleri

I am trying to graph two bands which represent the upper and lower error bounds for the function e

^{x}. I am using a 64 bit operating system so the absolute roundoff error is bounded by

|x-x^hat|<= (1/2)eps

_{M}|x|, where eps_M =2

^{-52}~2.2

^{-16}

So, in a worst case scenario the error, p (since e has already been used) equals .5*eps

_{M}.

Then,

y=e

^{x+px}

The important section of the code is below. It lists my experimentation results in the comments. I will try to upload the entire code later if I can’t solve this with some tweaking in the next hour. It uses a lot of graphics header files that I downloaded from Stroustrup's website since I am using his book to learn C++. No sense posting a bunch of headers if not necessary. I will double check the numerics theory, I haven't used it in awhile. However, I think that my confusion towards the results is from a lack of practice not theoretical understanding.

At the moment, I am only considering the upper bound, band_u, so all the comment refer to that. Thanks for your time and comments!

int fac(int n) //factorial(n); n! { int r=1; while (n>1){ r*=n; --n; } return r; } double term(double x, int n){ return pow(x,n)/fac(n); } //nth term of series double expe(double x, int n) //sum of n terms for x { double sum = 0; for (int i=0; i<n; ++i) sum+=term(x,i); return sum; } int expN_number_of_terms=10; double expN(double x) { return expe(x,expN_number_of_terms); } const double eps_M=2^(-52); //machine error ~ 2.204*10^-16 //eps_M=10^(-9) yields a cusp graph which "cusps" at (0,1); //the behavior after the y-axis is unexpected. The graph should be exponentially growing. //eps_M=10^(-11) the graph looks similar to y=e^x as expected, except //acts as a lower bound instead of an upper bound //eps_M=10^(-12) yields a flatline on y=1 after y-axis //eps_M=10^(-14) yeilds a cusp again //eps_M=2^(-52) yields an extremely sharp cusp; almost flatline on x-axis double band_u(double x){ //upper errror bound double y; y=exp(x+.5*abs(x)*eps_M); return y; } double band_l(double x){ //lower error bound double y; y=exp(x-.5*abs(x)*eps_M); return y; }

Forgot to make clear that this function is the upper bound,

y=e

^{x+px}