Pregunta

I have a code for doing integrations and in the last for loop starting on line 67 I have a for loop which accumulates the values of the function at randomly generated points to get the integral(Monte Carlo integration). Unfortunately after the loop finishes I get NAN as the result for "monte2" variable. I have written a printf statement inside the for loop to pinpoint the mistake only to notice that after 235.494781 the sum turns into -nan. What may be the reason behind this problem? I am running Ubuntu 12.04.3 LTS 32-bit and compile the plain C code with GCC version 4.6.3. I appreciate your help, the code is as follows:

P.S: The code is originally written in Code Blocks on Windows 8 64-bit by a friend of mine if this makes a difference.

#include<stdio.h>
#include<math.h>

float I1(float x)
{
  return exp(-x)*cos(x*x)*cos(x*x);
}

float I2(float t)
{
  return cos(log(t)*log(t))*cos(log(t)*log(t));
}
float random()
{
  float a;
  a=rand()%1000;
  a=a/1000*20;
  //  printf("%.15f\t%f\n",I1(a),a);
  return a;
}
float random2()
{
  float a;
  a=rand()%1000;
  a/=1000;
  //  printf("%.15f\t%f\n",I2(a),a);
  return a;
}

int main()
{
  FILE *data=fopen("data.txt","w");
  FILE *data2=fopen("data2.txt","w");

  float trap=0,monte=0,sum=0, monte2=0;
  float a[1000],b[1000],dt=0.005;
  int i;

  /* Part 1 */

  for(i=0;i<1000;i++)
    a[i]=I1(i*dt);
  for(i=0;i<1000;i++)
    fprintf(data,"%f\t%f\n",i*dt,a[i]);

  for(i=1;i<1000;i++)
    trap+=(a[i]+a[i-1])/2*dt;
  printf("The integral value of I1 is = %f with trapezoid rule\n",trap);


  for(i=0;i<500;i++)
    monte+=I1(random());
  printf("The Monte Carlo Technique value for I1 is %f with 500 samples\n",monte/500*20);

  /* Part 2 */
  dt=0.001;
   printf(" \n");
  for(i=1;i<=1000;i++)
    b[i]=I2(i*dt);
  for(i=1;i<=1000;i++)
    fprintf(data2,"%f\t%f\n",i*dt,b[i]);

  for(i=2;i<=1000;i++)
    trap+=(b[i]+b[i-1])/2*dt;
  printf("The integral value of I2 is = %f with trapezoid rule\n",trap/2);

  for(i=0;i<500;i++)
  {
    monte2+=I2(random2());
    printf("%f \n", monte2);
  }
  printf("The Monte Carlo Technique value of I2 is %f with 500 samples\n",monte2/500);
  printf("\n");
  printf("Comment 1: Two values obtained with trapezoid rule is close to each other;however,they are not exactly same.\n");
  printf("\n");
  printf("Comment 2: The integral value and monte carlo value of I1 is closer than the integral value and monte carlo value of I2.This means that we have better expectation value of I1 with monte carlo technique with 500 samples.\n");
  fclose(data2);
  fclose(data);
  return 0;
}
¿Fue útil?

Solución

Your function call

monte2+=I2(random2());  

may produce NaN. This is because random2 may returns 0. log 0 is infinity. This will cause cos(log(t)*log(t))*cos(log(t)*log(t)) to produce NaN.

See the graph for log function:

enter image description here

Note that the graph gets arbitrarily close to the y axis, but does not meet or intersect it1.


1. Source Wikipedia

Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top