6 Replies - 288 Views - Last Post: 22 September 2012 - 07:02 PM Rate Topic: -----

#1 InterzoneAgent  Icon User is offline

  • New D.I.C Head

Reputation: 1
  • View blog
  • Posts: 20
  • Joined: 19-May 12

Scientific Computing: possible underflow

Posted 22 September 2012 - 11:49 AM

Reference texts: Programming:Principles and Practice by Stroustrup
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 ex. I am using a 64 bit operating system so the absolute roundoff error is bounded by

|x-x^hat|<= (1/2)epsM|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*epsM.

Then,

y=ex+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=ex+px

Is This A Good Question/Topic? 0
  • +

Replies To: Scientific Computing: possible underflow

#2 Skydiver  Icon User is offline

  • Code herder
  • member icon

Reputation: 3548
  • View blog
  • Posts: 10,989
  • Joined: 05-May 12

Re: Scientific Computing: possible underflow

Posted 22 September 2012 - 12:37 PM

Have you tried using a long double instead of just a double?

If you need more precision than that, there is also GMP: http://gmplib.org/
Was This Post Helpful? 0
  • +
  • -

#3 InterzoneAgent  Icon User is offline

  • New D.I.C Head

Reputation: 1
  • View blog
  • Posts: 20
  • Joined: 19-May 12

Re: Scientific Computing: possible underflow

Posted 22 September 2012 - 03:16 PM

Thanks. I am trying to get a function in the headers Stroustrup supplied to take either a double or a long double. Hopefully, when I try it out with the
long double
type that will fix things.
Was This Post Helpful? 0
  • +
  • -

#4 InterzoneAgent  Icon User is offline

  • New D.I.C Head

Reputation: 1
  • View blog
  • Posts: 20
  • Joined: 19-May 12

Re: Scientific Computing: possible underflow

Posted 22 September 2012 - 06:08 PM

I tried long double. Still the same results, so it must not be a precision problem.
Was This Post Helpful? 0
  • +
  • -

#5 Skydiver  Icon User is offline

  • Code herder
  • member icon

Reputation: 3548
  • View blog
  • Posts: 10,989
  • Joined: 05-May 12

Re: Scientific Computing: possible underflow

Posted 22 September 2012 - 06:41 PM

Is the version of pow() you are using taking and returning long doubles? Or does it simply truncate the long double down to a double, and returning a double? Same question for your abs().
Was This Post Helpful? 0
  • +
  • -

#6 InterzoneAgent  Icon User is offline

  • New D.I.C Head

Reputation: 1
  • View blog
  • Posts: 20
  • Joined: 19-May 12

Re: Scientific Computing: possible underflow

Posted 22 September 2012 - 06:52 PM

I fixed the problem, apparently multiplying by a small number causes much more error than I had expected. I knew that dividing by small numbers did. Here is my new code. I divide by machine epsilon which means that I can't use the proper representation of machine epsilon, but I avoid massive roundoff error. I would appreciate help on ways of maintaining code which represents the idea I am trying to illustrate while avoiding massive roundoff error. Thanks!

Here are my musings. Perhaps if I simply change the order of multiplication when dealing with small numbers? 10^-16 should be able to fit in long double, so I see no reason for underflow.

const long double eps_M=10^(16);		//machine error ~ 2.204*10^-16

long double band_u(long double x){	//upper errror bound
	long double y;
	y=exp(x+0.5*abs(2)/eps_M);
	//y=(abs(x)/2)/eps_M;
	
	return y;
}


Success! At the mere cost of a pseudo-magic constant.
Was This Post Helpful? 0
  • +
  • -

#7 InterzoneAgent  Icon User is offline

  • New D.I.C Head

Reputation: 1
  • View blog
  • Posts: 20
  • Joined: 19-May 12

Re: Scientific Computing: possible underflow

Posted 22 September 2012 - 07:02 PM

I missed your post earlier. Yeah, I was using ^ instead of pow. No doubt that was my problem. Thanks!!
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1