skip to Main Content

I’m writing an integration script, and it uses a non-uniform grid. I’m using the pow function in math.h to calculate by what each value of x should be multiplied to get to the next value of x (this value is xmult in my code). However, for some reason, my pow function outputs zero, even though my calculator says that it should not be. Here are the relevant snippets from my code:

#include "stdio.h"
#include "math.h"
...
int N = 10240;
double recN = 1/N;
double upperlim = 4;
double lowerlim = 0.0001;
double limratio = upperlim/lowerlim;
double xmult;
xmult = pow(limratio,recN);
printf("%en",xmult);

The output of that printf function is 1.000000e+00 even though my calculated value does not round to that. I am compiling using gcc on CentOS with the -lm flag, it does not compile without that flag set. If it is any more help, the full code is below (the final output of the program being 0.000000e+00).

#include "stdio.h"
#include "math.h"

int main() {
  int N = 10240;
  double recN = 1/N;
  double pi = M_PI;
  double upperlim = 4;
  double lowerlim = 0.0001;
  double limratio = upperlim/lowerlim;
  double h;
  double f;
  double fnext;
  double ftot = 0;
  double ftot2 = 0;
  double x=lowerlim;
  double xmult;
  double xnext;
  int i;
  double func(double x){
    return sin(x) * exp(-x);
    }
  fnext = func(lowerlim);
  xmult = pow(limratio,recN);
  printf("%en",xmult);
  for (i=0;i<N;i++){
    f = fnext;
    fnext = func(x*xmult);
    xnext = x * xmult;
    h = xnext - x;
    ftot = 0.5*h*(f+fnext);
    ftot2 = ftot2  + ftot;
    x = xnext;
  }
  printf("%en",ftot2);
 }

2

Answers


  1. More than likely, your problem is a common error when dealing with numbers in C. Consider this:

    double recN = 1/N;

    Here, you define the variable as double, but you then perform integer division. If N>1 then this will result in zero.

    Instead, look for places you’ve done this type of calculation and convert it to:

    double recN = 1.0/N;

    Login or Signup to reply.
  2.   int N = 10240;
      double recN = 1/N;
    

    1 and N are both of type int, so the division is integer division, which rounds to an integer. Thus 1/N is 0.

    You probably want double recN = 1.0/N;.

    Login or Signup to reply.
Please signup or login to give your own answer.
Back To Top
Search