Problem with integration script

[eazye]

Hi,

I'm trying to code a program that performs an analysis of the integral I = 0.9sin(pi*x) in the interval 0.13<x<0.87 using crude Monte Carlo methods. I began with seeding the pseudo-random number generator and then use the equation I=(0.87-0.13)*(average f(x)) to calculate an estimate of I. The same calculation is repeated 100 times for instances where x is sampled N=10, N=100 and N=1000 times. The estimates for N=10, N=100 and N=1000 are then averaged and compared. The script also includes a calculation of the error Itheoretical - I estimate.

However, when you run the code you will notice that the estimate of I with N=10 is almost always as accurate as that of N=1000 and both are more always accurate N=100. This simply contradicts probability theory where the accuracy of the estimated value of I should increase with increasing sample number N.

I suspect the problem might be with the way I am using the random number generator but I cannot identify where exactly I went wrong. Using arrays was also not the best idea but I do not understand vectors. However, I know that arrays can sometimes pose problems - might they be the cause of my problem here?

I would really, really appreciate if someone could help me in identifying where I made the mistake.

#include <cmath>
#include <iostream>
#include <fstream>

#define PI 3.14159265

using namespace std;

int main()
{
    srand((unsigned)time(0));
    
    float sine10[10]={0.0}, sine100[100]={0.0}, sine1000[1000]={0.0}; ///sin
    float a = 0.87;/// upper limit
    float b = 0.13; /// lower limit
    float N10[10]={0.0}; ///N=10 array
    float N100[100]={0.0}; /// N=100 array
    float N1000[1000]={0.0}; /// N=1000 array
    float n10avg = 0, n100avg = 0, n1000avg = 0; /// averages
    float err10[100]={0.0};
    float err100[100]={0.0};  //// errors
    float err1000[100]={0.0};
    float sumerr10 =0.0, sumerr100 = 0.0, sumerr1000 = 0.0; ////sum of errors
    float avgerr10 = 0.0, avgerr100 = 0.0, avgerr1000 = 0.0; /// average error
    int Nten = 10, Nhun = 100, Nthou = 1000; /// N=10 N=100 N=1000
    float sumI10 = 0.0, sumI100 = 0.0, sumI1000 = 0.0;
    
    float Itru[100];  //// reference value for integral I
    for(int i = 0; i < 100; i++)
    {
        Itru[i]=0.525833;
    }
    
    for(int i = 0; i < 100; i++) /// loop required to obtain 100 estimates I
    {
        float sumn10i[10]={0.0};
        float avgfx10[100] = {0.0};
        float I10[100]={0.0}; /// array for 100 estimates of I for N = 10
        
        float sumn100i[100]={0.0};
        float avgfx100[100] = {0.0};
        float I100[100]={0.0}; /// N = 100
        
        float sumn1000i[1000]={0.0};
        float avgfx1000[100] = {0.0};
        float I1000[100]={0.0}; /// N = 1000
        
        for(int n = 0; n < 10; n++)
        {
            float array[n];
            float x;
            x = rand() % 740; /// obtain 10 random values
            array[n]=(x+13)/1000; ///limit the 10 random values between 0.13 and 0.87
            sine10[n] = sin (PI * array[n]); /// calculate the sin for every x
            N10[n] = 0.9*sine10[n];
            sumn10i[i]= sumn10i[i] + N10[n]; /// sum 10 N values
        }
        
        for(int n = 0; n < 100; n++)
        {
            float array[n];
            float x;
            x = rand() % 740;
            array[n]=(x+13)/1000;
            sine100[n] = sin (PI * array[n]);
            N100[n] = 0.9*sine100[n];
            sumn100i[i]= sumn100i[i] + N100[n];
        }

        for(int n = 0; n < 1000; n++)
        {
            float array[n];
            float x;
            x = rand() % 740;
            array[n]=(x+13)/1000;
            sine1000[n] = sin (PI * array[n]);
            N1000[n] = 0.9*sine1000[n];
            sumn1000i[i]= sumn1000i[i] + N1000[n];
        }
        
        avgfx10[i] = sumn10i[i]/Nten; //// average estimate of f(x) for every N (N=10)
        I10[i] = (a - b) * avgfx10[i]; /// estimate I based on I = (a-b)[average f(x)]
        sumI10 = sumI10 + I10[i];
        err10[i]=Itru[i]-I10[i];
        sumerr10 = sumerr10+err10[i];
        
        avgfx100[i] = sumn100i[i]/Nhun; // N = 100
        I100[i] = (a - b) * avgfx100[i];
        sumI100 = sumI100 + I100[i];
        err100[i]=Itru[i]-I100[i];
        sumerr100 = sumerr100+err100[i];
        
        avgfx1000[i] = sumn1000i[i]/Nthou; /// N = 1000
        I1000[i] = (a - b) * avgfx1000[i];
        sumI1000 = sumI10 + I1000[i];
        err1000[i]=Itru[i]-I1000[i];
        sumerr1000 = sumerr1000+err1000[i];
        
        //cout << sumn10i[i]<<endl;
        //cout << sumn100i[i]<< endl<<endl;
        
    }
    
    n10avg = sumI10/100; //// Average I value for N = 10
    avgerr10=abs(sumerr10/100); // Average |I-Iest|
    
    cout << "The average estimate of I for N = 10: " << n10avg << endl <<
    "The average |I-Iest| = " << avgerr10 << endl;
    
    n100avg = sumI100/100; /// N = 100
    avgerr100=abs(sumerr100/100);
    
    cout << "The average estimate of I for N = 100: " << n100avg << endl <<
    "The average |I-Iest| = " << avgerr100 << endl;
    
    n1000avg = sumI1000/100; /// N = 1000
    avgerr1000=abs(sumerr1000/100);
    
    cout << "The average estimate of I for N = 1000: " << n1000avg << endl <<
    "The average |I-Iest| = " << avgerr1000 << endl;
    
    
    return 0;